/[MITgcm]/MITgcm_contrib/darwin2/pkg/monod/monod_forcing.F
ViewVC logotype

Diff of /MITgcm_contrib/darwin2/pkg/monod/monod_forcing.F

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph | View Patch Patch

revision 1.1 by jahn, Wed Apr 13 18:56:25 2011 UTC revision 1.9 by jahn, Mon Jul 30 15:21:51 2012 UTC
# Line 21  C======================================= Line 21  C=======================================
21  #include "PARAMS.h"  #include "PARAMS.h"
22  #include "GRID.h"  #include "GRID.h"
23  #include "DYNVARS.h"  #include "DYNVARS.h"
24  #ifdef USE_QSW  c for Qsw and/or surfaceForcingT
25    c choice which field to take pCO2 from for pCO2limit
26    c this assumes we use Ttendency from offline
27  #include "FFIELDS.h"  #include "FFIELDS.h"
 #endif  
28  #ifdef ALLOW_LONGSTEP  #ifdef ALLOW_LONGSTEP
29  #include "LONGSTEP.h"  #include "LONGSTEP.h"
30  #endif  #endif
# Line 43  c ANNA include wavebands_params.h Line 44  c ANNA include wavebands_params.h
44  #include "WAVEBANDS_PARAMS.h"  #include "WAVEBANDS_PARAMS.h"
45  #endif  #endif
46    
47    
48  C     === Global variables ===  C     === Global variables ===
49  c tracers  c tracers
50         _RL Ptr(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy,nDarwin)         _RL Ptr(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy,nDarwin)
# Line 68  c plankton arrays Line 70  c plankton arrays
70        _RL  Phy_k(npmax,Nr)        _RL  Phy_k(npmax,Nr)
71        _RL  Phyup(npmax)        _RL  Phyup(npmax)
72        _RL  part_k(Nr)        _RL  part_k(Nr)
73    #ifdef ALLOW_CDOM
74          _RL cdom_k(Nr)
75    #endif
76  c iron partitioning  c iron partitioning
77        _RL  freefe(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL  freefe(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
78  c some working variables  c some working variables
# Line 93  c ANNA define variables for wavebands Line 98  c ANNA define variables for wavebands
98         _RL Edz(tlam,Nr),Esz(tlam,Nr),Euz(tlam,Nr),Eutop(tlam,Nr)         _RL Edz(tlam,Nr),Esz(tlam,Nr),Euz(tlam,Nr),Eutop(tlam,Nr)
99         _RL tirrq(nr)         _RL tirrq(nr)
100         _RL tirrwq(tlam,nr)         _RL tirrwq(tlam,nr)
101           _RL c1(tlam,nr), c2(tlam,nr)
102         _RL solz         _RL solz
103         _RL rmud         _RL rmud
104         _RL actot,bctot,bbctot         _RL actot,bctot,bbctot
# Line 111  C      always need for diagnostics Line 117  C      always need for diagnostics
117        _RL  Diver2(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL  Diver2(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
118        _RL  Diver3(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL  Diver3(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
119        _RL  Diver4(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL  Diver4(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
120          _RL Shannon(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
121          _RL Simpson(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
122    
123        _RL  tmpphy(npmax)        _RL  tmpphy(npmax)
124        _RL  totphy, biotot, maxphy, phymax        _RL  totphy, biotot, maxphy, phymax
# Line 124  C      always need for diagnostics Line 132  C      always need for diagnostics
132        _RL chlup(npmax)        _RL chlup(npmax)
133  #endif  #endif
134  #endif  #endif
135    #ifdef ALLOW_CDOM
136          _RL cdoml  
137          _RL dcdoml
138    #endif
139    
140  #ifdef ALLOW_DIAGNOSTICS  #ifdef ALLOW_DIAGNOSTICS
141  COJ for diagnostics  COJ for diagnostics
# Line 185  c some local variables Line 197  c some local variables
197         _RL PSiupl         _RL PSiupl
198         _RL Tlocal         _RL Tlocal
199         _RL Slocal         _RL Slocal
200           _RL pCO2local
201         _RL Qswlocal         _RL Qswlocal
202         _RL NH4l         _RL NH4l
203         _RL NO2l         _RL NO2l
# Line 248  c tendencies Line 261  c tendencies
261         _RL do2l         _RL do2l
262         _RL dZooCl(nzmax)         _RL dZooCl(nzmax)
263  c air-sea fluxes  c air-sea fluxes
264         _RL flxCO2(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)         _RL flxCO2(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
265         _RL flxALK(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)         _RL flxALK(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
266         _RL flxO2(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)         _RL flxO2(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
267  #endif  #endif
268    
269         _RL tot_Nfix         _RL tot_Nfix
# Line 277  c Line 290  c
290             Diver2(i,j,k)=0. _d 0             Diver2(i,j,k)=0. _d 0
291             Diver3(i,j,k)=0. _d 0             Diver3(i,j,k)=0. _d 0
292             Diver4(i,j,k)=0. _d 0             Diver4(i,j,k)=0. _d 0
293               Shannon(i,j,k)=0. _d 0
294               Simpson(i,j,k)=1. _d 0
295  #endif  #endif
296    
297  #ifdef ALLOW_DIAGNOSTICS  #ifdef ALLOW_DIAGNOSTICS
298  COJ for diagnostics  COJ for diagnostics
299             PParr(i,j,k) = 0. _d 0             PParr(i,j,k) = 0. _d 0
300             Nfixarr(i,j,k) = 0. _d 0             Nfixarr(i,j,k) = 0. _d 0
301    #ifdef DAR_DIAG_CHL
302               GeiderChlarr(i,j,k) = 0. _d 0
303               GeiderChl2Carr(i,j,k) = 0. _d 0
304               DoneyChlarr(i,j,k) = 0. _d 0
305               DoneyChl2Carr(i,j,k) = 0. _d 0
306               CloernChlarr(i,j,k) = 0. _d 0
307               CloernChl2Carr(i,j,k) = 0. _d 0
308    #endif
309  c ANNA_TAVE  c ANNA_TAVE
310  #ifdef WAVES_DIAG_PCHL  #ifdef WAVES_DIAG_PCHL
311             DO np=1,npmax             DO np=1,npmax
# Line 405  C ========================== i,j loops = Line 428  C ========================== i,j loops =
428  c ------------ these are convenient ------------------------------------  c ------------ these are convenient ------------------------------------
429           DO k=1,Nr           DO k=1,Nr
430            part_k(k)  = max(Ptr(i,j,k,bi,bj,iPOP),0. _d 0)            part_k(k)  = max(Ptr(i,j,k,bi,bj,iPOP),0. _d 0)
431    #ifdef ALLOW_CDOM
432              cdom_k(k)  = max(Ptr(i,j,k,bi,bj,iCDOM),0. _d 0)
433    #endif
434            DO np = 1,npmax            DO np = 1,npmax
435              Phy_k(np,k) = max(Ptr(i,j,k,bi,bj,iPhy+np-1),0. _d 0)              Phy_k(np,k) = max(Ptr(i,j,k,bi,bj,iPhy+np-1),0. _d 0)
436  #ifdef GEIDER  #ifdef GEIDER
# Line 420  c ------------ these are convenient ---- Line 446  c ------------ these are convenient ----
446  c ------------ GET CDOM_k FOR WAVEBANDS_3D and RADTRANS ----------------  c ------------ GET CDOM_k FOR WAVEBANDS_3D and RADTRANS ----------------
447  #ifdef WAVEBANDS  #ifdef WAVEBANDS
448  #if defined(DAR_CALC_ACDOM) || defined(DAR_RADTRANS)  #if defined(DAR_CALC_ACDOM) || defined(DAR_RADTRANS)
449    #ifdef ALLOW_CDOM
450            call MONOD_ACDOM(cdom_k,
451         O                     acdom_k,
452         I                     myThid)
453    #else
454           call MONOD_ACDOM(phychl_k,aphy_chl,aw,           call MONOD_ACDOM(phychl_k,aphy_chl,aw,
455       O                     acdom_k,       O                     acdom_k,
456       I                     myThid)       I                     myThid)
457    #endif
458  #else  #else
459           DO k=1,Nr           DO k=1,Nr
460            DO ilam = 1,tlam            DO ilam = 1,tlam
# Line 602  c           add water and CDOM Line 634  c           add water and CDOM
634       I                darwin_radtrans_kmax,darwin_radtrans_niter,       I                darwin_radtrans_kmax,darwin_radtrans_niter,
635       O                Edz,Esz,Euz,Eutop,       O                Edz,Esz,Euz,Eutop,
636       O                tirrq,tirrwq,       O                tirrq,tirrwq,
637         O                c1,c2,
638       I                myThid)       I                myThid)
639  #else  #else
640  c dzlocal ?????  c dzlocal ?????
# Line 657  c brute force... Line 690  c brute force...
690               psil  = max(Ptr(i,j,k,bi,bj,iPOSi ),0. _d 0)               psil  = max(Ptr(i,j,k,bi,bj,iPOSi ),0. _d 0)
691               NH4l  = max(Ptr(i,j,k,bi,bj,iNH4  ),0. _d 0)               NH4l  = max(Ptr(i,j,k,bi,bj,iNH4  ),0. _d 0)
692               NO2l  = max(Ptr(i,j,k,bi,bj,iNO2  ),0. _d 0)               NO2l  = max(Ptr(i,j,k,bi,bj,iNO2  ),0. _d 0)
693    #ifdef ALLOW_CDOM
694                 cdoml = max(Ptr(i,j,k,bi,bj,iCDOM  ),0. _d 0)
695    #endif
696  #ifdef ALLOW_CARBON  #ifdef ALLOW_CARBON
697               dicl  = max(Ptr(i,j,k,bi,bj,iDIC  ),0. _d 0)               dicl  = max(Ptr(i,j,k,bi,bj,iDIC  ),0. _d 0)
698               docl  = max(Ptr(i,j,k,bi,bj,iDOC  ),0. _d 0)               docl  = max(Ptr(i,j,k,bi,bj,iDOC  ),0. _d 0)
# Line 727  c ratio of maximum species Line 763  c ratio of maximum species
763                      Diver4(i,j,k)=Diver4(i,j,k)+1. _d 0                      Diver4(i,j,k)=Diver4(i,j,k)+1. _d 0
764                   endif                   endif
765                 enddo                 enddo
766    c           totphy > thresh0
767                endif
768    c Shannon and Simpson indices
769                Shannon(i,j,k) = 0. _d 0
770    c note: minimal valid value is 1, but we set to zero below threshold
771                Simpson(i,j,k) = 0. _d 0
772                if (totphy.gt.shannon_thresh) then
773                  do np=1,npmax
774                   if (Phy(np) .gt. 0. _d 0) then
775                    tmpphy(np) = Phy(np)/totphy
776                    Shannon(i,j,k)=Shannon(i,j,k)+tmpphy(np)*LOG(tmpphy(np))
777                    Simpson(i,j,k)=Simpson(i,j,k)+tmpphy(np)*tmpphy(np)
778                   endif
779                  enddo
780                  Shannon(i,j,k) = -Shannon(i,j,k)
781                  Simpson(i,j,k) = 1./Simpson(i,j,k)
782              endif              endif
783  #endif  #endif
784    
# Line 763  c for explicit sinking of particulate ma Line 815  c for explicit sinking of particulate ma
815               Slocal = salt(i,j,k,bi,bj)               Slocal = salt(i,j,k,bi,bj)
816  #endif  #endif
817    
818    c choice where to get pCO2 from
819    c taking from igsm dic run - fed through Tflux array
820    c               pCO2local=surfaceForcingT(i,j,bi,bj)
821    c or from darwin carbon module
822    #ifdef ALLOW_CARBON
823                   pCO2local=pCO2(i,j,bi,bj)
824    #else
825                   pCO2local=280. _d -6
826    #endif
827    
828               freefu = max(freefe(i,j,k),0. _d 0)               freefu = max(freefe(i,j,k),0. _d 0)
829               if (k.eq.1) then               if (k.eq.1) then
830                 inputFel = inputFe(i,j,bi,bj)                 inputFel = inputFe(i,j,bi,bj)
# Line 807  c set tendencies to 0 Line 869  c set tendencies to 0
869                 dphychl(np)=0. _d 0                 dphychl(np)=0. _d 0
870               enddo               enddo
871  #endif  #endif
872    #ifdef ALLOW_CDOM
873                dcdoml=0. _d 0
874    #endif
875  #ifdef ALLOW_CARBON  #ifdef ALLOW_CARBON
876               ddicl=0. _d 0               ddicl=0. _d 0
877               ddocl=0. _d 0               ddocl=0. _d 0
# Line 876  c ANNA pass extra variables if WAVEBANDS Line 941  c ANNA pass extra variables if WAVEBANDS
941       I                       pofeupl, psiupl,       I                       pofeupl, psiupl,
942       I                       PARl,       I                       PARl,
943       I                       Tlocal, Slocal,       I                       Tlocal, Slocal,
944         I                       pCO2local,
945       I                       freefu, inputFel,       I                       freefu, inputFel,
946       I                       bottom, dzlocal,       I                       bottom, dzlocal,
947       O                       Rstarl, RNstarl,       O                       Rstarl, RNstarl,
# Line 906  c ANNA pass extra variables if WAVEBANDS Line 972  c ANNA pass extra variables if WAVEBANDS
972       I                       dphychl,       I                       dphychl,
973       I                       chlup,       I                       chlup,
974  #endif  #endif
975    #ifdef ALLOW_CDOM
976         O                       dcdoml,
977         I                       cdoml,
978    #endif
979  #ifdef WAVEBANDS  #ifdef WAVEBANDS
980       I                       PARw_k(1,k),       I                       PARw_k(1,k),
981  #endif  #endif
# Line 931  c            endif Line 1001  c            endif
1001  c  c
1002  #ifdef IRON_SED_SOURCE  #ifdef IRON_SED_SOURCE
1003  c only above minimum depth (continental shelf)  c only above minimum depth (continental shelf)
1004               if (rF(k).lt.depthfesed) then               if (rF(k).gt.-depthfesed) then
1005  c only if bottom layer  c only if bottom layer
1006                 if (bottom.eq.1.0 _d 0) then                 if (bottom.eq.1.0 _d 0) then
1007  #ifdef IRON_SED_SOURCE_VARIABLE  #ifdef IRON_SED_SOURCE_VARIABLE
# Line 966  c Line 1036  c
1036               picupl = PICl               picupl = PICl
1037  c include surface forcing  c include surface forcing
1038               if (k.eq.1) then               if (k.eq.1) then
1039                ddicl  =  ddicl  + flxCO2(i,j,bi,bj)                ddicl  =  ddicl  + flxCO2(i,j)
1040                dalkl  =  dalkl  + flxALK(i,j,bi,bj)                dalkl  =  dalkl  + flxALK(i,j)
1041                do2l   =  do2l   + flxO2(i,j,bi,bj)                do2l   =  do2l   + flxO2(i,j)
1042               endif               endif
1043  #endif  #endif
1044  c  c
# Line 1112  c in darwin_plankton this is stored for Line 1182  c in darwin_plankton this is stored for
1182  #endif  #endif
1183  #endif  #endif
1184           ENDDO           ENDDO
1185    #ifdef ALLOW_CDOM
1186              Ptr(i,j,k,bi,bj,iCDOM ) = Ptr(i,j,k,bi,bj,iCDOM ) +
1187         &                                  dtplankton*dcdoml
1188    #endif
1189  #ifdef ALLOW_CARBON  #ifdef ALLOW_CARBON
1190            Ptr(i,j,k,bi,bj,iDIC ) = Ptr(i,j,k,bi,bj,iDIC ) +            Ptr(i,j,k,bi,bj,iDIC ) = Ptr(i,j,k,bi,bj,iDIC ) +
1191       &                                  dtplankton*ddicl       &                                  dtplankton*ddicl
# Line 1267  Coj            no Eu at surface (yet) Line 1341  Coj            no Eu at surface (yet)
1341       &                                 Eutop(ilam,k)*dtplankton       &                                 Eutop(ilam,k)*dtplankton
1342               enddo               enddo
1343  #endif  #endif
1344    #ifdef DAR_DIAG_IRR_AMPS
1345                 do ilam = 1,tlam
1346                   c1ave(i,j,k,bi,bj,ilam)=c1ave(i,j,k,bi,bj,ilam)+
1347         &                                 c1(ilam,k)*dtplankton
1348                   c2ave(i,j,k,bi,bj,ilam)=c2ave(i,j,k,bi,bj,ilam)+
1349         &                                 c2(ilam,k)*dtplankton
1350                 enddo
1351    #endif
1352  #ifdef DAR_DIAG_ABSORP  #ifdef DAR_DIAG_ABSORP
1353               do ilam = 1,tlam               do ilam = 1,tlam
1354                 aave(i,j,k,bi,bj,ilam)=aave(i,j,k,bi,bj,ilam)+                 aave(i,j,k,bi,bj,ilam)=aave(i,j,k,bi,bj,ilam)+
# Line 1291  Coj            no Eu at surface (yet) Line 1373  Coj            no Eu at surface (yet)
1373       &                                 bbpart_k(k,ilam)*dtplankton       &                                 bbpart_k(k,ilam)*dtplankton
1374               enddo               enddo
1375  #endif  #endif
1376    #ifdef DAR_RADTRANS
1377                 if (k.eq.1) then
1378                   rmudave(i,j,bi,bj)=rmudave(i,j,bi,bj)+
1379         &                                 rmud*dtplankton
1380                 endif
1381    #endif
1382  #ifdef DAR_DIAG_RSTAR  #ifdef DAR_DIAG_RSTAR
1383               do np=1,npmax               do np=1,npmax
1384                 Rstarave(i,j,k,bi,bj,np)=Rstarave(i,j,k,bi,bj,np)+                 Rstarave(i,j,k,bi,bj,np)=Rstarave(i,j,k,bi,bj,np)+
# Line 1331  Coj            no Eu at surface (yet) Line 1419  Coj            no Eu at surface (yet)
1419  #ifdef ALLOW_CARBON  #ifdef ALLOW_CARBON
1420               if (k.eq.1) then               if (k.eq.1) then
1421                SURave(i,j,bi,bj)    =SURave(i,j,bi,bj)+                SURave(i,j,bi,bj)    =SURave(i,j,bi,bj)+
1422       &                              flxCO2(i,j,bi,bj)*dtplankton       &                              flxCO2(i,j)*dtplankton
1423                SURCave(i,j,bi,bj)    =SURCave(i,j,bi,bj)+                SURCave(i,j,bi,bj)    =SURCave(i,j,bi,bj)+
1424       &                              FluxCO2(i,j,bi,bj)*dtplankton       &                              FluxCO2(i,j,bi,bj)*dtplankton
1425                SUROave(i,j,bi,bj)   =SUROave(i,j,bi,bj)+                SUROave(i,j,bi,bj)   =SUROave(i,j,bi,bj)+
1426       &                              flxO2(i,j,bi,bj)*dtplankton       &                              flxO2(i,j)*dtplankton
1427                pCO2ave(i,j,bi,bj)   =pCO2ave(i,j,bi,bj)+                pCO2ave(i,j,bi,bj)   =pCO2ave(i,j,bi,bj)+
1428       &                              pCO2(i,j,bi,bj)*dtplankton       &                              pCO2(i,j,bi,bj)*dtplankton
1429                pHave(i,j,bi,bj)     =pHave(i,j,bi,bj)+                pHave(i,j,bi,bj)     =pHave(i,j,bi,bj)+
# Line 1438  c ANNA end TAVE Line 1526  c ANNA end TAVE
1526          WRITE(diagname,'(A8)') 'Diver4  '          WRITE(diagname,'(A8)') 'Diver4  '
1527          CALL DIAGNOSTICS_FILL( Diver4(1-Olx,1-Oly,1), diagname,          CALL DIAGNOSTICS_FILL( Diver4(1-Olx,1-Oly,1), diagname,
1528       &                         0,Nr,2,bi,bj,myThid )       &                         0,Nr,2,bi,bj,myThid )
1529            WRITE(diagname,'(A8)') 'Shannon '
1530            CALL DIAGNOSTICS_FILL( Shannon(1-Olx,1-Oly,1), diagname,
1531         &                         0,Nr,2,bi,bj,myThid )
1532            WRITE(diagname,'(A8)') 'Simpson '
1533            CALL DIAGNOSTICS_FILL( Simpson(1-Olx,1-Oly,1), diagname,
1534         &                         0,Nr,2,bi,bj,myThid )
1535  #endif  #endif
1536  #ifdef ALLOW_DIAZ  #ifdef ALLOW_DIAZ
1537  #ifdef DAR_DIAG_NFIXP  #ifdef DAR_DIAG_NFIXP
# Line 1463  c ANNA end TAVE Line 1557  c ANNA end TAVE
1557       &                         0,Nr,2,bi,bj,myThid )       &                         0,Nr,2,bi,bj,myThid )
1558  #endif  #endif
1559  #ifdef ALLOW_CARBON  #ifdef ALLOW_CARBON
1560          CALL DIAGNOSTICS_FILL( flxCO2(1-Olx,1-Oly,bi,bj), 'DICTFLX ',          CALL DIAGNOSTICS_FILL( flxCO2(1-Olx,1-Oly), 'DICTFLX ',
1561       &                         0,1,2,bi,bj,myThid )       &                         0,1,2,bi,bj,myThid )
1562          CALL DIAGNOSTICS_FILL( FluxCO2(1-Olx,1-Oly,bi,bj), 'DICCFLX ',          CALL DIAGNOSTICS_FILL( FluxCO2(1-Olx,1-Oly,bi,bj), 'DICCFLX ',
1563       &                         0,1,2,bi,bj,myThid )       &                         0,1,2,bi,bj,myThid )
1564          CALL DIAGNOSTICS_FILL( flxO2(1-Olx,1-Oly,bi,bj), 'DICOFLX ',          CALL DIAGNOSTICS_FILL( flxO2(1-Olx,1-Oly), 'DICOFLX ',
1565       &                         0,1,2,bi,bj,myThid )       &                         0,1,2,bi,bj,myThid )
1566          CALL DIAGNOSTICS_FILL( pCO2(1-Olx,1-Oly,bi,bj), 'DICPCO2 ',          CALL DIAGNOSTICS_FILL( pCO2(1-Olx,1-Oly,bi,bj), 'DICPCO2 ',
1567       &                         0,1,2,bi,bj,myThid )       &                         0,1,2,bi,bj,myThid )
Legend:
Removed from v.1.1  
changed lines
  Added in v.1.9
  ViewVC Help
Powered by ViewVC 1.1.22