pkg/seaice/seaice_tracer_phys.F

f50f58ec54 Gael*0001 #include "SEAICE_OPTIONS.h"
                0002
                0003 CStartOfInterface
                0004       SUBROUTINE SEAICE_TRACER_PHYS( myTime, myIter, myThid )
dc0a9dc58e Jean*0005 C     *=======================================================*
                0006 C     | SUBROUTINE seaice_tracer_phys
                0007 C     | o Time step SItr/SItrEFF as a result of
                0008 C     |   seaice thermodynamics and specific tracer physics
                0009 C     *=======================================================*
f50f58ec54 Gael*0010       IMPLICIT NONE
                0011
                0012 C     === Global variables ===
                0013 #include "SIZE.h"
                0014 #include "EEPARAMS.h"
f1c91a595d Jean*0015 #include "PARAMS.h"
f50f58ec54 Gael*0016 #include "FFIELDS.h"
                0017 #include "DYNVARS.h"
                0018 #include "SEAICE_SIZE.h"
                0019 #include "SEAICE.h"
                0020 #include "SEAICE_PARAMS.h"
                0021 #include "SEAICE_TRACER.h"
8bc8bee483 Gael*0022 #ifdef ALLOW_SALT_PLUME
                0023 # include "SALT_PLUME.h"
                0024 #endif
f50f58ec54 Gael*0025
                0026 C     === Routine arguments ===
                0027 C     INPUT:
                0028 C     myTime  :: Simulation time
                0029 C     myIter  :: Simulation timestep number
                0030 C     myThid  :: Thread no. that called this routine.
                0031 C     OUTPUT:
                0032       _RL myTime
                0033       INTEGER myIter, myThid
                0034 CEndOfInterface
                0035
                0036 C     === Local variables ===
                0037 #ifdef ALLOW_SITRACER
                0038
                0039       INTEGER iTr, jTh, I, J, bi, bj, ks
                0040       _RL SItrFromOcean  (1:sNx,1:sNy)
bb24b8a3e6 Gael*0041       _RL SItrFromFlood   (1:sNx,1:sNy)
f50f58ec54 Gael*0042       _RL HEFFprev, HEFFpost, growFact, meltPart, tmpscal1
bb24b8a3e6 Gael*0043       _RL SItrExpand  (1:sNx,1:sNy)
                0044       _RL AREAprev, AREApost, expandFact
3721cfe5e4 Gael*0045       CHARACTER*8   diagName
f50f58ec54 Gael*0046
3721cfe5e4 Gael*0047 #ifdef ALLOW_SITRACER_DEBUG_DIAG
f50f58ec54 Gael*0048       _RL DIAGarray     (1:sNx,1:sNy,Nr)
                0049 #endif
                0050
                0051 cgf for now I do not fully account for ocean-ice fluxes of tracer
                0052 cgf -> I just prescribe it consistent with age tracer
                0053 cgf eventually I will need to handle them as function params
                0054
                0055       ks=1
                0056
                0057       DO bj=myByLo(myThid),myByHi(myThid)
                0058       DO bi=myBxLo(myThid),myBxHi(myThid)
8bc8bee483 Gael*0059       DO iTr=1,SItrNumInUse
f50f58ec54 Gael*0060
fcdcf3fa1d Jean*0061 c 0) set ice-ocean and ice-snow exchange values
f50f58ec54 Gael*0062 c =============================================
                0063       DO J=1,sNy
                0064        DO I=1,sNx
8bc8bee483 Gael*0065         SItrFromOcean(i,j)=SItrFromOcean0(iTr)
                0066         SItrFromFlood(i,j)=SItrFromFlood0(iTr)
                0067         SItrExpand(i,j)=SItrExpand0(iTr)
f50f58ec54 Gael*0068        ENDDO
                0069       ENDDO
fcdcf3fa1d Jean*0070 c salinity tracer:
dc0a9dc58e Jean*0071       if ( (SItrName(iTr).EQ.'salinity').AND.
8bc8bee483 Gael*0072      &      (SItrFromOceanFrac(iTr).GT.ZERO) ) then
f50f58ec54 Gael*0073        DO J=1,sNy
                0074         DO I=1,sNx
8bc8bee483 Gael*0075          SItrFromOcean(i,j)=SItrFromOceanFrac(iTr)*salt(I,j,ks,bi,bj)
                0076          SItrFromFlood(i,j)=SItrFromFloodFrac(iTr)*salt(I,j,ks,bi,bj)
f50f58ec54 Gael*0077         ENDDO
                0078        ENDDO
                0079       endif
                0080 c 1) seaice thermodynamics processes
                0081 c ==================================
bb24b8a3e6 Gael*0082       if (SItrMate(iTr).EQ.'HEFF') then
f50f58ec54 Gael*0083       DO J=1,sNy
                0084        DO I=1,sNx
                0085         HEFFprev=SItrHEFF(i,j,bi,bj,1)
3721cfe5e4 Gael*0086 #ifdef ALLOW_SITRACER_DEBUG_DIAG
fcdcf3fa1d Jean*0087         DIAGarray(I,J,5+(iTr-1)*5) =
f50f58ec54 Gael*0088      &    HEFFprev*SItracer(i,j,bi,bj,iTr) + SItrBucket(i,j,bi,bj,iTr)
                0089 #endif
f61838dfc1 Torg*0090 c apply the sequence of thermodynamics increments to actual tracer
f50f58ec54 Gael*0091 c (see seaice_growth.F)
                0092 c (jTh=1 tendency due to ice-ocean interaction)
                0093 c (jTh=2 tendency due to the atmosphere, over ice covered part)
                0094 c (jTh=3 tendency due to the atmosphere, over open water part)
                0095 c (jTh=4 tendency due to flooding)
                0096         DO jTh=1,3
                0097          HEFFprev=SItrHEFF(i,j,bi,bj,jTh)
                0098          HEFFpost=SItrHEFF(i,j,bi,bj,jTh+1)
                0099 c compute ratio in [0. 1.] range for either growth or melt
                0100          growFact=1. _d 0
                0101          meltPart=0. _d 0
                0102          if (HEFFpost.GT.HEFFprev) growFact=HEFFprev/HEFFpost
                0103          if (HEFFpost.LT.HEFFprev) meltPart=HEFFprev-HEFFpost
                0104 c update SItr accordingly
                0105          SItracer(i,j,bi,bj,iTr)=SItracer(i,j,bi,bj,iTr)*growFact
                0106      &                      +SItrFromOcean(i,j)*(1. _d 0 - growFact)
                0107          SItrBucket(i,j,bi,bj,iTr)=SItrBucket(i,j,bi,bj,iTr)
                0108      &             -HEFFpost*SItrFromOcean(i,j)*(1. _d 0 - growFact)
                0109          SItrBucket(i,j,bi,bj,iTr)=SItrBucket(i,j,bi,bj,iTr)
                0110      &                            +meltPart*SItracer(i,j,bi,bj,iTr)
                0111         ENDDO
                0112 c apply flooding term
                0113         growFact=1. _d 0
                0114         HEFFprev=SItrHEFF(i,j,bi,bj,4)
                0115         HEFFpost=SItrHEFF(i,j,bi,bj,5)
                0116         if (HEFFpost.GT.HEFFprev) growFact=HEFFprev/HEFFpost
                0117         SItracer(i,j,bi,bj,iTr)=SItracer(i,j,bi,bj,iTr)*growFact
bb24b8a3e6 Gael*0118      &                     +SItrFromFlood(i,j) *(1. _d 0 - growFact)
                0119 c rk: flooding can only imply an ocean-ice tracer exchange, as long
                0120 c as we dont have snow tracers, so it goes through SItrBucket.
f50f58ec54 Gael*0121           SItrBucket(i,j,bi,bj,iTr)=SItrBucket(i,j,bi,bj,iTr)
bb24b8a3e6 Gael*0122      &             -HEFFpost*SItrFromFlood(i,j)*(1. _d 0 - growFact)
3721cfe5e4 Gael*0123 #ifdef ALLOW_SITRACER_DEBUG_DIAG
fcdcf3fa1d Jean*0124         DIAGarray(I,J,5+(iTr-1)*5) = HEFFpost*SItracer(i,j,bi,bj,iTr)
3721cfe5e4 Gael*0125      &  +SItrBucket(i,j,bi,bj,iTr)-DIAGarray(I,J,5+(iTr-1)*5)
f50f58ec54 Gael*0126 #endif
                0127        ENDDO
                0128       ENDDO
bb24b8a3e6 Gael*0129 c TAF?      if (SItrMate(iTr).EQ.'AREA') then
                0130       else
                0131 c 1) or seaice cover expansion
                0132 c ============================
fcdcf3fa1d Jean*0133 c this is much simpler than for ice volume/mass tracers, because
                0134 c properties of the ice surface are not be conserved across the
bb24b8a3e6 Gael*0135 c ocean-ice system, the contraction/expansion terms are all
fcdcf3fa1d Jean*0136 c simultaneous (which is sane), and the only generic effect
bb24b8a3e6 Gael*0137 c is due to expansion (new cover).
                0138       DO J=1,sNy
                0139        DO I=1,sNx
                0140 c apply expansion
                0141         AREAprev=SItrAREA(i,j,bi,bj,2)
                0142         AREApost=SItrAREA(i,j,bi,bj,3)
                0143 c compute ratio in [0. 1.] range for expansion/contraction
                0144         expandFact=1. _d 0
                0145         if (AREApost.GT.AREAprev) expandFact=AREAprev/AREApost
                0146 c update SItr accordingly
                0147          SItracer(i,j,bi,bj,iTr)=SItracer(i,j,bi,bj,iTr)*expandFact
                0148      &                      +SItrExpand(i,j)*(1. _d 0 - expandFact)
                0149        ENDDO
                0150       ENDDO
                0151       endif
f50f58ec54 Gael*0152 c 2) very ice tracer processes
                0153 c ============================
                0154       if (SItrName(iTr).EQ.'age') then
                0155 c age tracer: grow old as time passes by
                0156        DO J=1,sNy
                0157         DO I=1,sNx
bb24b8a3e6 Gael*0158           if (( (SItrHEFF(i,j,bi,bj,5).GT.0. _d 0).AND.(SItrMate(iTr)
                0159      &     .EQ.'HEFF') ).OR.( (SItrAREA(i,j,bi,bj,3).GT.0. _d 0).AND.
                0160      &     (SItrMate(iTr).EQ.'AREA') )) then
f50f58ec54 Gael*0161             SItracer(i,j,bi,bj,iTr)=
                0162      &      SItracer(i,j,bi,bj,iTr)+SEAICE_deltaTtherm
                0163           else
                0164             SItracer(i,j,bi,bj,iTr)=0. _d 0
                0165           endif
                0166         ENDDO
                0167        ENDDO
                0168       elseif (SItrName(iTr).EQ.'salinity') then
                0169 c salinity tracer: no specific process
                0170       elseif (SItrName(iTr).EQ.'one') then
                0171 c "ice concentration" tracer: no specific process
bb24b8a3e6 Gael*0172       elseif (SItrName(iTr).EQ.'ridge') then
                0173 c simple, made up, ice surface roughness index prototype
                0174        DO J=1,sNy
                0175         DO I=1,sNx
                0176 c ridging increases roughness
                0177           SItracer(i,j,bi,bj,iTr)=SItracer(i,j,bi,bj,iTr)+
                0178      &    MAX(0. _d 0, SItrAREA(i,j,bi,bj,1)-SItrAREA(i,j,bi,bj,2))
                0179 c ice melt reduces ridges/roughness
                0180           HEFFprev=SItrHEFF(i,j,bi,bj,1)
                0181           HEFFpost=SItrHEFF(i,j,bi,bj,4)
                0182           tmpscal1=1. _d 0
                0183           if (HEFFprev.GT.HEFFpost) tmpscal1=HEFFpost/HEFFprev
                0184           SItracer(i,j,bi,bj,iTr)=SItracer(i,j,bi,bj,iTr)*tmpscal1
                0185         ENDDO
                0186        ENDDO
f50f58ec54 Gael*0187       endif
bb24b8a3e6 Gael*0188 c 3) ice-ocean tracer exchange/mapping to external variables
                0189 c ==========================================================
3721cfe5e4 Gael*0190 #ifdef ALLOW_DIAGNOSTICS
f1c91a595d Jean*0191       IF ( useDiagnostics .AND. SItrMate(iTr).EQ.'HEFF') THEN
2eab512e2a Gael*0192         WRITE(diagName,'(A4,I2.2,A2)') 'SItr',iTr,'Fx'
3721cfe5e4 Gael*0193         tmpscal1=-ONE/SEAICE_deltaTtherm*SEAICE_rhoIce
dc0a9dc58e Jean*0194         CALL DIAGNOSTICS_SCALE_FILL(SItrBucket(1-OLx,1-OLy,bi,bj,iTr),
3721cfe5e4 Gael*0195      &   tmpscal1, 1, diagName,0,1,2,bi,bj,myThid)
f1c91a595d Jean*0196       ENDIF
3721cfe5e4 Gael*0197 #endif
8bc8bee483 Gael*0198
                0199       if ( (SItrName(iTr).EQ.'salinity').AND.
                0200      &     (SEAICE_salinityTracer) ) then
f50f58ec54 Gael*0201 c salinity tracer: salt flux
8bc8bee483 Gael*0202         DO J=1,sNy
                0203          DO I=1,sNx
                0204           saltFlux(I,J,bi,bj) = - SItrBucket(i,j,bi,bj,iTr)
                0205      &      *HEFFM(I,J,bi,bj)/SEAICE_deltaTtherm*SEAICE_rhoIce
f50f58ec54 Gael*0206 c note: at this point of the time step, that is the correct sign
8bc8bee483 Gael*0207 #ifdef ALLOW_SALT_PLUME
                0208 c should work for both constant and variable ice salinity -- to be tested
dc0a9dc58e Jean*0209           saltPlumeFlux(I,J,bi,bj) = MAX(ZERO,saltFlux(I,J,bi,bj))
8bc8bee483 Gael*0210      &      *SPsalFRAC*(salt(I,j,ks,bi,bj)-SItrFromOcean(i,j))
                0211 #endif
                0212          ENDDO
                0213         ENDDO
f50f58ec54 Gael*0214       endif
dc0a9dc58e Jean*0215
f50f58ec54 Gael*0216       DO J=1,sNy
                0217        DO I=1,sNx
3721cfe5e4 Gael*0218 #ifdef ALLOW_SITRACER_DEBUG_DIAG
                0219         DIAGarray(I,J,4+(iTr-1)*5) = - SItrBucket(i,j,bi,bj,iTr)
f50f58ec54 Gael*0220      &  *HEFFM(I,J,bi,bj)/SEAICE_deltaTtherm*SEAICE_rhoIce
                0221 #endif
                0222 c empty bucket
f61838dfc1 Torg*0223 c  but not for 'grease' (see seaice_growth.F)
0320e25227 Mart*0224         if (SItrName(iTr).NE.'grease')
f61838dfc1 Torg*0225      &      SItrBucket(i,j,bi,bj,iTr)=0. _d 0
f50f58ec54 Gael*0226        ENDDO
                0227       ENDDO
dc0a9dc58e Jean*0228
bb24b8a3e6 Gael*0229 c TAF? elseif (SItrMate(iTr).EQ.'AREA') then
8bc8bee483 Gael*0230
f50f58ec54 Gael*0231 c 4) diagnostics
                0232 c ==============
3721cfe5e4 Gael*0233 #ifdef ALLOW_SITRACER_DEBUG_DIAG
bb24b8a3e6 Gael*0234       if (SItrMate(iTr).EQ.'HEFF') then
f50f58ec54 Gael*0235       DO J=1,sNy
                0236        DO I=1,sNx
                0237         HEFFpost=SItrHEFF(i,j,bi,bj,5)
                0238         DIAGarray(I,J,1+(iTr-1)*5) = SItracer(i,j,bi,bj,iTr)
                0239         DIAGarray(I,J,2+(iTr-1)*5) = SItracer(i,j,bi,bj,iTr)*HEFFpost
3721cfe5e4 Gael*0240 c DIAGarray(:,:,3) is the term of comparison for DIAGarray(:,:,2)
e54fe3e1f9 Gael*0241         if (SItrName(iTr).EQ.'salinity') then
f50f58ec54 Gael*0242           DIAGarray(I,J,3+(iTr-1)*5) = HSALT(i,j,bi,bj)/SEAICE_rhoIce
                0243         elseif (SItrName(iTr).EQ.'one') then
                0244           DIAGarray(I,J,3+(iTr-1)*5) = HEFFpost
                0245         endif
3721cfe5e4 Gael*0246 c DIAGarray(:,:,4) allows check of conservation : del(SItrBucket)+del(SItr*HEFF)=0. over do_phys
                0247 c DIAGarray(:,:,5) is the tracer flux from the ocean (<0 incr. ocean tracer)
f50f58ec54 Gael*0248        ENDDO
                0249       ENDDO
bb24b8a3e6 Gael*0250       else
                0251       DO J=1,sNy
                0252        DO I=1,sNx
                0253         AREApost=SItrAREA(i,j,bi,bj,3)
                0254         DIAGarray(I,J,1+(iTr-1)*5) = SItracer(i,j,bi,bj,iTr)
                0255         DIAGarray(I,J,2+(iTr-1)*5) = SItracer(i,j,bi,bj,iTr)*AREApost
                0256        ENDDO
                0257       ENDDO
                0258       endif
f50f58ec54 Gael*0259 #endif
                0260       ENDDO
3721cfe5e4 Gael*0261 #ifdef ALLOW_SITRACER_DEBUG_DIAG
fcdcf3fa1d Jean*0262 c     CALL DIAGNOSTICS_FILL(DIAGarray,'UDIAG1  ',0,Nr,3,bi,bj,myThid)
f50f58ec54 Gael*0263 #endif
                0264       ENDDO
                0265       ENDDO
                0266
                0267 #endif /* ALLOW_SITRACER */
                0268
                0269       RETURN
                0270       END