C $Header: /u/gcmpack/MITgcm/pkg/atm2d/relax_add.F,v 1.5 2012/01/20 20:32:57 jmc Exp $
C $Name: $
#include "ctrparam.h"
#include "ATM2D_OPTIONS.h"
C !INTERFACE:
SUBROUTINE RELAX_ADD( wght0, wght1,
& intime0, intime1, iftime, myIter, myThid)
C *==========================================================*
C | Adds restoring terms to surface forcing. Note that: |
C | - restoring is phased out as restor (or act.) SST <2C |
C | - if nsTypeRelax NE 0, salt rest. phased out nr poles |
C | - if ntTypeRelax NE 0, temp rest. phased out nr poles |
C *==========================================================*
IMPLICIT NONE
#include "ATMSIZE.h"
#include "SIZE.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "GRID.h"
#include "THSICE_VARS.h"
#include "ATM2D_VARS.h"
c include ocean and seaice vars
C !INPUT/OUTPUT PARAMETERS:
C === Routine arguments ===
C wght0, wght1 - weight of first and second month, respectively
C intime0,intime1- month id # for first and second months
C iftime - true -> prompts a reloading of data from disk
C myIter - Ocean iteration number
C myThid - Thread no. that called this routine.
_RL wght0
_RL wght1
INTEGER intime0
INTEGER intime1
LOGICAL iftime
INTEGER myIter
INTEGER myThid
C LOCAL VARIABLES:
C Save below so that continual file reloads aren't necessary
COMMON /OCEANRELAX/
& sst0, sst1, sss0, sss1
_RS sst0(1-Olx:sNx+Olx,1-Oly:sNy+Oly,1,1)
_RS sst1(1-Olx:sNx+Olx,1-Oly:sNy+Oly,1,1)
_RS sss0(1-Olx:sNx+Olx,1-Oly:sNy+Oly,1,1)
_RS sss1(1-Olx:sNx+Olx,1-Oly:sNy+Oly,1,1)
_RL lambdaTheta,lambdaSalt
_RS nearIce ! constant used to phase out rest near frz point
_RL qrelflux, frelflux
_RL sstRelax(1:sNx,1:sNy) ! relaxation sst as computed from file
_RL sssRelax(1:sNx,1:sNy) ! relaxation sss as computed from file
INTEGER i,j
IF (ifTime) THEN
C If the above condition is met then we need to read in
C data for the period ahead and the period behind current time.
WRITE(*,*) 'S/R RELAX_ADD: Reading new data'
IF ( thetaRelaxFile .NE. ' ' ) THEN
CALL READ_REC_XY_RS( thetaRelaxFile,sst0,intime0,
& myIter,myThid )
CALL READ_REC_XY_RS( thetaRelaxFile,sst1,intime1,
& myIter,myThid )
ENDIF
IF ( saltRelaxFile .NE. ' ' ) THEN
CALL READ_REC_XY_RS( saltRelaxFile,sss0,intime0,
& myIter,myThid )
CALL READ_REC_XY_RS( saltRelaxFile,sss1,intime1,
& myIter,myThid )
ENDIF
ENDIF
IF ((thetaRelaxFile.NE.' ').OR.(saltRelaxFile.NE.' ')) THEN
C-- Interpolate and add to anomaly
DO j=1,sNy
IF (ntTypeRelax .EQ. 0) THEN
lambdaTheta = r_tauThetaRelax
ELSE
lambdaTheta = r_tauThetaRelax *
& max(cos(1.5 _d 0*yC(1,j,1,1)*deg2rad),0. _d 0)
ENDIF
IF (nsTypeRelax .EQ. 0) THEN
lambdaSalt = r_tauSaltRelax
ELSE
lambdaSalt = r_tauSaltRelax *
& max(cos(1.5 _d 0*yC(1,j,1,1)*deg2rad),0. _d 0)
ENDIF
DO i=1,sNx
IF (maskC(i,j,1,1,1) .EQ. 1.) THEN
IF (thetaRelaxFile.NE.' ') THEN
sstRelax(i,j)= (wght0*sst0(i,j,1,1) + wght1*sst1(i,j,1,1))
ELSE !no T restoring; use actual SST to determine if nr freezing
sstRelax(i,j)= sstFromOcn(i,j)
ENDIF
IF (saltRelaxFile.NE.' ') THEN
sssRelax(i,j)= (wght0*sss0(i,j,1,1) + wght1*sss1(i,j,1,1))
ELSE ! no S restoring; this ensures frelflux=0
sssRelax(i,j)= sssFromOcn(i,j)
ENDIF
C Next lines: linearly phase out SST restoring between 2C and -1C
C ONLY if seaice is present
IF ((sstRelax(i,j).GT.2. _d 0).OR.
& (iceMask(i,j,1,1) .EQ. 0. _d 0)) THEN
nearIce=1.0
ELSEIF (sstRelax(i,j) .LE. -1. _d 0) THEN
nearIce=0.0
ELSE
nearIce=(sstRelax(i,j)+1.0)/3.0
endif
qrelflux= lambdaTheta*(sstFromOcn(i,j)-sstRelax(i,j))
& * (HeatCapacity_Cp*rhoNil*drF(1))*nearIce
C- should use rhoConst instead of rhoNil:
c & * (HeatCapacity_Cp*rhoConst*drF(1))*nearIce
C no longer restore on top of ice, but effectively full gridpoint UNDER ice
C (unless gridpoint is fully ice covered)
IF (iceMask(i,j,1,1) .LT. 0.999 _d 0) THEN
qneto_2D(i,j)= qneto_2D(i,j) + qrelflux
& / (1. _d 0 - iceMask(i,j,1,1))
ENDIF
C qneto_2D(i,j)= qneto_2D(i,j) + qrelflux
C qneti_2D(i,j)= qneti_2D(i,j) + qrelflux
frelflux= -lambdaSalt*(sssFromOcn(i,j)-sssRelax(i,j))/
& (convertFW2Salt *recip_drF(1))*nearIce
C or use actual salt instead of convertFW2salt above?
IF (frelflux .GT. 0. _d 0) THEN
evapo_2D(i,j)= evapo_2D(i,j) - frelflux
C note most of the time, evapi=0 when iceMask>0 anyway
C (i.e., only when relaxing SST >2 but ocn still ice-covered)
IF (iceMask(i,j,1,1).GT.0. _d 0)
& evapi_2D(i,j)= evapi_2D(i,j) - frelflux
ELSE
precipo_2D(i,j)= precipo_2D(i,j) + frelflux
IF (iceMask(i,j,1,1).GT.0. _d 0)
& precipi_2D(i,j)= precipi_2D(i,j) + frelflux
ENDIF
C IF (iceMask(i,j,1,1) .GT. 0. _d 0) THEN
C PRINT *,'Frelflux',frelflux,precipi_2D(i,j),atm_precip(j+1)
C ENDIF
C Diagnostics
sum_qrel(i,j)= sum_qrel(i,j) + qrelflux*dtatmo
sum_frel(i,j)= sum_frel(i,j) + frelflux*dtatmo
ENDIF
ENDDO
ENDDO
ENDIF
C PRINT *,'***bottom of relaxadd',wght0,wght1,intime0,intime1
C PRINT *,'evapo_2d: ',evapo_2D(JBUGI,JBUGJ)
C PRINT *,'precipo_2d: ',precipo_2D(JBUGI,JBUGJ)
C PRINT *,'qneto_2d: ',qneto_2D(JBUGI,JBUGJ)
C PRINT *,'SStfrom Ocn: ',sstfromocn(JBUGI,JBUGJ)
C PRINT *,'SSSfrom Ocn: ',sssfromocn(JBUGI,JBUGJ)
RETURN
END