tutorial_global_oce_latlon/code/ptracers_forcing_surf.F

4307c62f79 Jean*0001 #include "PTRACERS_OPTIONS.h"
                0002
                0003 CBOP
                0004 C !ROUTINE: PTRACERS_FORCING_SURF
                0005
                0006 C !INTERFACE: ==========================================================
                0007       SUBROUTINE PTRACERS_FORCING_SURF(
3f4989e08f Jean*0008      I                            relaxForcingS,
4307c62f79 Jean*0009      I                            bi, bj, iMin, iMax, jMin, jMax,
                0010      I                            myTime,myIter,myThid )
                0011
                0012 C !DESCRIPTION:
                0013 C     Precomputes surface forcing term for pkg/ptracers.
                0014 C     Precomputation is needed because of non-local KPP transport term,
                0015 C     routine KPP_TRANSPORT_PTR.
                0016
                0017 C !USES: ===============================================================
                0018       IMPLICIT NONE
                0019 #include "SIZE.h"
                0020 #include "EEPARAMS.h"
                0021 #include "PARAMS.h"
                0022 #include "GRID.h"
                0023 #include "SURFACE.h"
3f4989e08f Jean*0024 #include "DYNVARS.h"
4307c62f79 Jean*0025 #include "FFIELDS.h"
                0026 #include "PTRACERS_SIZE.h"
22e0cff85e Jean*0027 #include "PTRACERS_PARAMS.h"
d6215f7b79 Jean*0028 #include "PTRACERS_START.h"
22e0cff85e Jean*0029 #include "PTRACERS_FIELDS.h"
4307c62f79 Jean*0030
                0031 C !INPUT PARAMETERS: ===================================================
3f4989e08f Jean*0032 C  relaxForcingS        :: Salt forcing due to surface relaxation
4307c62f79 Jean*0033 C  bi,bj                :: tile indices
                0034 C  myTime               :: model time
                0035 C  myIter               :: time-step number
                0036 C  myThid               :: thread number
3f4989e08f Jean*0037       _RL relaxForcingS(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
4307c62f79 Jean*0038       INTEGER bi, bj, iMin, iMax, jMin, jMax
                0039       _RL myTime
                0040       INTEGER myIter
                0041       INTEGER myThid
                0042
                0043 #ifdef ALLOW_PTRACERS
                0044
                0045 C !LOCAL VARIABLES: ====================================================
                0046 C  i,j                  :: loop indices
                0047 C  iTrc                 :: tracer index
                0048 C  ks                   :: surface level index
                0049       INTEGER i, j
                0050       INTEGER iTrc, ks
3f4989e08f Jean*0051       _RL add2EmP(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0052       _RL epsil, cutoff, tmpVar
4307c62f79 Jean*0053 CEOP
                0054
                0055       IF ( usingPCoords ) THEN
                0056         ks = Nr
                0057       ELSE
                0058         ks = 1
                0059       ENDIF
                0060
                0061 C Example of how to add forcing at the surface
                0062       DO iTrc=1,PTRACERS_numInUse
d6215f7b79 Jean*0063 c       IF ( PTRACERS_StepFwd(iTrc) ) THEN
4307c62f79 Jean*0064           DO j = jMin, jMax
                0065            DO i = iMin, iMax
d6215f7b79 Jean*0066 C-    Case of age-tracer: at the surface, add 10.days relaxation towards zero:
22e0cff85e Jean*0067              surfaceForcingPTr(i,j,bi,bj,iTrc) =
62fd6ae4e5 Jean*0068      &        + 1. _d 0 / (10. _d 0 * 86400. _d 0)
                0069      &                  * ( 0. _d 0 - pTracer(i,j,ks,bi,bj,iTrc) )
97f6d163ab Davi*0070      &                  * drF(ks) * _hFacC(i,j,ks,bi,bj)
4307c62f79 Jean*0071            ENDDO
                0072           ENDDO
d6215f7b79 Jean*0073 c       ENDIF
4307c62f79 Jean*0074       ENDDO
                0075
3f4989e08f Jean*0076 C--   Option to convert Salt-relaxation into additional EmP contribution
                0077       IF ( PTRACERS_addSrelax2EmP ) THEN
                0078 C-    here we assume that salt_EvPrRn = 0
                0079 C     set cutoff value to prevent too large additional EmP:
                0080 C       current limit is set to 0.1 CFL
                0081         epsil = 1. _d -10
                0082         cutoff = 0.1 _d 0 *drF(ks)/PTRACERS_dTLev(ks)
                0083         IF ( ( (nonlinFreeSurf.GT.0 .OR. usingPCoords)
                0084      &         .AND. useRealFreshWaterFlux )
                0085      &     .OR.convertFW2Salt .EQ. -1. ) THEN
                0086          DO j = jMin, jMax
                0087           DO i = iMin, iMax
                0088             tmpVar = MAX( salt(i,j,ks,bi,bj), epsil )
                0089             add2EmP(i,j) = relaxForcingS(i,j)/tmpVar
                0090             add2EmP(i,j) = rUnit2mass
                0091      &                  *MAX( -cutoff, MIN( add2EmP(i,j), cutoff ) )
                0092           ENDDO
                0093          ENDDO
                0094         ELSE
                0095          DO j = jMin, jMax
                0096           DO i = iMin, iMax
                0097             add2EmP(i,j) = relaxForcingS(i,j)/convertFW2Salt
                0098             add2EmP(i,j) = rUnit2mass
                0099      &                  *MAX( -cutoff, MIN( add2EmP(i,j), cutoff ) )
                0100           ENDDO
                0101          ENDDO
                0102         ENDIF
                0103 #ifdef ALLOW_DIAGNOSTICS
                0104         IF ( useDiagnostics ) THEN
                0105          CALL DIAGNOSTICS_FILL(add2EmP,'Add2EmP ',0,1,2,bi,bj,myThid)
                0106         ENDIF
                0107 #endif /* ALLOW_DIAGNOSTICS */
                0108       ELSE
                0109         DO j = jMin, jMax
                0110           DO i = iMin, iMax
                0111             add2EmP(i,j) = 0. _d 0
                0112           ENDDO
                0113         ENDDO
                0114       ENDIF
                0115 C-- end of "addEmP" setting
                0116
4307c62f79 Jean*0117       IF ( (nonlinFreeSurf.GT.0 .OR. usingPCoords)
                0118      &     .AND. useRealFreshWaterFlux ) THEN
                0119
                0120        DO iTrc=1,PTRACERS_numInUse
                0121
                0122 c-  NonLin_FrSurf and RealFreshWaterFlux : PmEpR effectively changes
                0123 c   the water column height ; temp., salt, (tracer) flux associated
                0124 c   with this input/output of water is added here to the surface tendency.
                0125 c
d6215f7b79 Jean*0126          IF ( PTRACERS_StepFwd(iTrc) .AND.
                0127      &        PTRACERS_EvPrRn(iTrc).NE.UNSET_RL ) THEN
4307c62f79 Jean*0128           DO j = jMin, jMax
                0129            DO i = iMin, iMax
22e0cff85e Jean*0130              surfaceForcingPTr(i,j,bi,bj,iTrc) =
                0131      &          surfaceForcingPTr(i,j,bi,bj,iTrc)
3f4989e08f Jean*0132      &        + ( PmEpR(i,j,bi,bj) - add2EmP(i,j) )
4307c62f79 Jean*0133      &          *( PTRACERS_EvPrRn(iTrc) - pTracer(i,j,ks,bi,bj,iTrc) )
62fd6ae4e5 Jean*0134      &          *mass2rUnit
4307c62f79 Jean*0135            ENDDO
                0136           ENDDO
                0137          ENDIF
                0138
                0139        ENDDO
a09a74749d Davi*0140
                0141 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0142       ELSE
4307c62f79 Jean*0143
a09a74749d Davi*0144 C--   EmPmR does not really affect the water column height (for tracer budget)
                0145 C     and is converted to a salt tendency.
                0146
                0147        IF (convertFW2Salt .EQ. -1.) THEN
                0148 C-    use local surface tracer field to calculate forcing term:
                0149
                0150         DO iTrc=1,PTRACERS_numInUse
                0151
d6215f7b79 Jean*0152          IF ( PTRACERS_StepFwd(iTrc) .AND.
                0153      &        PTRACERS_EvPrRn(iTrc).NE.UNSET_RL ) THEN
22e0cff85e Jean*0154 C        account for Rain/Evap tracer content (PTRACERS_EvPrRn) using
a09a74749d Davi*0155 C        local surface tracer
                0156           DO j = jMin, jMax
                0157            DO i = iMin, iMax
22e0cff85e Jean*0158             surfaceForcingPTr(i,j,bi,bj,iTrc) =
                0159      &          surfaceForcingPTr(i,j,bi,bj,iTrc)
3f4989e08f Jean*0160      &        + ( EmPmR(i,j,bi,bj) + add2EmP(i,j) )
22e0cff85e Jean*0161      &          *( pTracer(i,j,ks,bi,bj,iTrc) - PTRACERS_EvPrRn(iTrc) )
a09a74749d Davi*0162      &          *mass2rUnit
                0163            ENDDO
                0164           ENDDO
                0165          ENDIF
                0166
                0167         ENDDO
                0168
                0169        ELSE
                0170 C-    use uniform tracer value to calculate forcing term:
                0171
                0172         DO iTrc=1,PTRACERS_numInUse
                0173
d6215f7b79 Jean*0174          IF ( PTRACERS_StepFwd(iTrc) .AND.
                0175      &        PTRACERS_EvPrRn(iTrc).NE.UNSET_RL ) THEN
22e0cff85e Jean*0176 C     account for Rain/Evap tracer content (PTRACERS_EvPrRn) assuming uniform
a09a74749d Davi*0177 C     surface tracer (=PTRACERS_ref)
                0178           DO j = jMin, jMax
                0179            DO i = iMin, iMax
22e0cff85e Jean*0180             surfaceForcingPTr(i,j,bi,bj,iTrc) =
                0181      &          surfaceForcingPTr(i,j,bi,bj,iTrc)
3f4989e08f Jean*0182      &        + ( EmPmR(i,j,bi,bj) + add2EmP(i,j) )
22e0cff85e Jean*0183      &            *( PTRACERS_ref(ks,iTrc) - PTRACERS_EvPrRn(iTrc) )
a09a74749d Davi*0184      &            *mass2rUnit
                0185            ENDDO
                0186           ENDDO
                0187          ENDIF
                0188
                0189         ENDDO
                0190
                0191 C-    end local-surface-tracer / uniform-value distinction
                0192        ENDIF
                0193
                0194       ENDIF
                0195
                0196 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0197
4307c62f79 Jean*0198 #endif /* ALLOW_PTRACERS */
                0199
                0200       RETURN
                0201       END