pkg/seaice/seaice_tracer_phys.F

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