Back to home page

MITgcm

 
 

    


File indexing completed on 2018-03-02 18:37:32 UTC

view on githubraw file Latest commit add29e06 on 2018-01-31 20:35:05 UTC
09a6f3668a Jeff*0001 #include "ctrparam.h"
                0002 #include "ATM2D_OPTIONS.h"
                0003 
                0004 C     !INTERFACE:
                0005       SUBROUTINE RELAX_ADD( wght0, wght1,
                0006      &               intime0, intime1, iftime, myIter, myThid)
                0007 C     *==========================================================*
                0008 C     | Adds restoring terms to surface forcing. Note that:      |
                0009 C     |    - restoring is phased out as restor (or act.) SST <2C |
                0010 C     |    - if nsTypeRelax NE 0, salt rest. phased out nr poles |
                0011 C     |    - if ntTypeRelax NE 0, temp rest. phased out nr poles |
                0012 C     *==========================================================*
                0013         IMPLICIT NONE
                0014 
                0015 #include "ATMSIZE.h"
                0016 #include "SIZE.h"
                0017 #include "EEPARAMS.h"
                0018 #include "PARAMS.h"
                0019 #include "GRID.h"
                0020 #include "THSICE_VARS.h"
                0021 #include "ATM2D_VARS.h"
                0022 
                0023 c include ocean and seaice vars
                0024 
                0025 C     !INPUT/OUTPUT PARAMETERS:
                0026 C     === Routine arguments ===
                0027 C     wght0, wght1   - weight of first and second month, respectively
                0028 C     intime0,intime1- month id # for first and second months
                0029 C     iftime - true -> prompts a reloading of data from disk
                0030 C     myIter - Ocean iteration number
                0031 C     myThid - Thread no. that called this routine.
                0032       _RL  wght0
                0033       _RL  wght1
                0034       INTEGER intime0
                0035       INTEGER intime1
                0036       LOGICAL iftime
                0037       INTEGER myIter
                0038       INTEGER myThid
                0039 
                0040 C     LOCAL VARIABLES:
                0041 C     Save below so that continual file reloads aren't necessary
                0042       COMMON /OCEANRELAX/
                0043      &                 sst0, sst1, sss0, sss1
                0044 
                0045       _RS  sst0(1-Olx:sNx+Olx,1-Oly:sNy+Oly,1,1)
9274434acc Jean*0046       _RS  sst1(1-Olx:sNx+Olx,1-Oly:sNy+Oly,1,1)
                0047       _RS  sss0(1-Olx:sNx+Olx,1-Oly:sNy+Oly,1,1)
                0048       _RS  sss1(1-Olx:sNx+Olx,1-Oly:sNy+Oly,1,1)
09a6f3668a Jeff*0049       _RL lambdaTheta,lambdaSalt
                0050       _RS nearIce    ! constant used to phase out rest near frz point
                0051       _RL qrelflux, frelflux
                0052       _RL sstRelax(1:sNx,1:sNy) ! relaxation sst as computed from file
                0053       _RL sssRelax(1:sNx,1:sNy) ! relaxation sss as computed from file
                0054       INTEGER i,j
                0055 
                0056       IF (ifTime) THEN
                0057 
                0058 C      If the above condition is met then we need to read in
                0059 C      data for the period ahead and the period behind current time.
                0060 
                0061         WRITE(*,*) 'S/R RELAX_ADD: Reading new data'
                0062         IF ( thetaRelaxFile .NE. ' '  ) THEN
                0063           CALL READ_REC_XY_RS( thetaRelaxFile,sst0,intime0,
                0064      &                      myIter,myThid )
                0065           CALL READ_REC_XY_RS( thetaRelaxFile,sst1,intime1,
                0066      &                      myIter,myThid )
                0067         ENDIF
                0068         IF ( saltRelaxFile .NE. ' '  ) THEN
                0069           CALL READ_REC_XY_RS( saltRelaxFile,sss0,intime0,
                0070      &                      myIter,myThid )
                0071           CALL READ_REC_XY_RS( saltRelaxFile,sss1,intime1,
                0072      &                      myIter,myThid )
                0073         ENDIF
                0074 
                0075       ENDIF
                0076 
                0077       IF ((thetaRelaxFile.NE.' ').OR.(saltRelaxFile.NE.' ')) THEN
                0078 
                0079 C--   Interpolate and add to anomaly
                0080       DO j=1,sNy
                0081 
                0082         IF (ntTypeRelax .EQ. 0) THEN
                0083           lambdaTheta =  r_tauThetaRelax
                0084         ELSE
b7373c84c8 Jeff*0085           lambdaTheta = r_tauThetaRelax *
09a6f3668a Jeff*0086      &                 max(cos(1.5 _d 0*yC(1,j,1,1)*deg2rad),0. _d 0)
                0087         ENDIF
                0088         IF (nsTypeRelax .EQ. 0) THEN
                0089           lambdaSalt = r_tauSaltRelax
                0090         ELSE
b7373c84c8 Jeff*0091           lambdaSalt = r_tauSaltRelax *
09a6f3668a Jeff*0092      &                max(cos(1.5 _d 0*yC(1,j,1,1)*deg2rad),0. _d 0)
                0093         ENDIF
                0094 
                0095         DO i=1,sNx
                0096 
                0097           IF (maskC(i,j,1,1,1) .EQ. 1.) THEN
                0098 
                0099           IF (thetaRelaxFile.NE.' ') THEN
                0100             sstRelax(i,j)= (wght0*sst0(i,j,1,1) + wght1*sst1(i,j,1,1))
                0101           ELSE  !no T restoring; use actual SST to determine if nr freezing
                0102             sstRelax(i,j)= sstFromOcn(i,j)
                0103           ENDIF
                0104 
                0105           IF (saltRelaxFile.NE.' ') THEN
                0106             sssRelax(i,j)= (wght0*sss0(i,j,1,1) + wght1*sss1(i,j,1,1))
                0107           ELSE  ! no S restoring; this ensures frelflux=0
                0108             sssRelax(i,j)= sssFromOcn(i,j)
                0109           ENDIF
                0110 
                0111 
                0112 C         Next lines: linearly phase out SST restoring between 2C and -1C
                0113 C         ONLY if seaice is present
                0114           IF ((sstRelax(i,j).GT.2. _d 0).OR.
                0115      &        (iceMask(i,j,1,1) .EQ. 0. _d 0)) THEN
                0116               nearIce=1.0
                0117           ELSEIF (sstRelax(i,j) .LE. -1. _d 0) THEN
                0118               nearIce=0.0
                0119           ELSE
                0120               nearIce=(sstRelax(i,j)+1.0)/3.0
                0121           endif
9274434acc Jean*0122 
50ec6123fe Jean*0123           qrelflux= lambdaTheta*(sstFromOcn(i,j)-sstRelax(i,j))
                0124      &            * (HeatCapacity_Cp*rhoNil*drF(1))*nearIce
                0125 C-    should use rhoConst instead of rhoNil:
                0126 c    &            * (HeatCapacity_Cp*rhoConst*drF(1))*nearIce
09a6f3668a Jeff*0127 
709ddc5d16 Jeff*0128 C         no longer restore on top of ice, but effectively full gridpoint UNDER ice
                0129 C         (unless gridpoint is fully ice covered)
                0130           IF (iceMask(i,j,1,1) .LT. 0.999 _d 0) THEN
                0131                qneto_2D(i,j)= qneto_2D(i,j) + qrelflux
                0132      &                / (1. _d 0 - iceMask(i,j,1,1))
                0133           ENDIF
                0134 C          qneto_2D(i,j)= qneto_2D(i,j) + qrelflux
                0135 C          qneti_2D(i,j)= qneti_2D(i,j) + qrelflux
09a6f3668a Jeff*0136 
                0137           frelflux= -lambdaSalt*(sssFromOcn(i,j)-sssRelax(i,j))/
                0138      &                  (convertFW2Salt *recip_drF(1))*nearIce
                0139 
                0140 C         or use actual salt instead of convertFW2salt above?
                0141 
                0142           IF (frelflux .GT. 0. _d 0) THEN
                0143             evapo_2D(i,j)= evapo_2D(i,j) - frelflux
                0144 C           note most of the time, evapi=0 when iceMask>0 anyway
                0145 C           (i.e., only when relaxing SST >2 but ocn still ice-covered)
                0146             IF (iceMask(i,j,1,1).GT.0. _d 0)
                0147      &            evapi_2D(i,j)= evapi_2D(i,j) - frelflux
                0148           ELSE
                0149             precipo_2D(i,j)= precipo_2D(i,j) + frelflux
                0150             IF (iceMask(i,j,1,1).GT.0. _d 0)
                0151      &            precipi_2D(i,j)= precipi_2D(i,j) + frelflux
                0152           ENDIF
                0153 
                0154 C          IF (iceMask(i,j,1,1) .GT. 0. _d 0) THEN
                0155 C          PRINT *,'Frelflux',frelflux,precipi_2D(i,j),atm_precip(j+1)
                0156 C          ENDIF
                0157 
                0158 C         Diagnostics
                0159           sum_qrel(i,j)= sum_qrel(i,j) + qrelflux*dtatmo
                0160           sum_frel(i,j)= sum_frel(i,j) + frelflux*dtatmo
                0161 
                0162           ENDIF
                0163         ENDDO
                0164       ENDDO
                0165       ENDIF
                0166 
                0167 C      PRINT *,'***bottom of relaxadd',wght0,wght1,intime0,intime1
                0168 C      PRINT *,'evapo_2d: ',evapo_2D(JBUGI,JBUGJ)
                0169 C      PRINT *,'precipo_2d: ',precipo_2D(JBUGI,JBUGJ)
                0170 C      PRINT *,'qneto_2d: ',qneto_2D(JBUGI,JBUGJ)
                0171 C      PRINT *,'SStfrom Ocn: ',sstfromocn(JBUGI,JBUGJ)
                0172 C      PRINT *,'SSSfrom Ocn: ',sssfromocn(JBUGI,JBUGJ)
                0173 
                0174       RETURN
                0175       END