/[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.4 by jahn, Mon Nov 7 21:55:11 2011 UTC revision 1.11 by jahn, Thu Aug 23 21:49:33 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 90  c ANNA define variables for wavebands Line 95  c ANNA define variables for wavebands
95  #ifdef DAR_RADTRANS  #ifdef DAR_RADTRANS
96         integer iday,iyr,imon,isec,lp,wd,mydate(4)         integer iday,iyr,imon,isec,lp,wd,mydate(4)
97         _RL Edwsf(tlam),Eswsf(tlam)         _RL Edwsf(tlam),Eswsf(tlam)
98         _RL Edz(tlam,Nr),Esz(tlam,Nr),Euz(tlam,Nr),Eutop(tlam,Nr)         _RL Edz(tlam,Nr),Esz(tlam,Nr),Euz(tlam,Nr)
99           _RL Estop(tlam,Nr),Eutop(tlam,Nr)
100         _RL tirrq(nr)         _RL tirrq(nr)
101         _RL tirrwq(tlam,nr)         _RL tirrwq(tlam,nr)
102           _RL amp1(tlam,nr), amp2(tlam,nr)
103         _RL solz         _RL solz
104         _RL rmud         _RL rmud
105         _RL actot,bctot,bbctot         _RL actot,bctot,bbctot
106         _RL apart_k(Nr,tlam),bpart_k(Nr,tlam),bbpart_k(Nr,tlam)         _RL apart_k(Nr,tlam),bpart_k(Nr,tlam),bbpart_k(Nr,tlam)
107         _RL bt_k(Nr,tlam), bb_k(Nr,tlam)         _RL bt_k(Nr,tlam), bb_k(Nr,tlam)
108           _RL discEs, discEu
109  #else  #else
110         _RL PARwdn(tlam)         _RL PARwdn(tlam)
111  #endif  #endif
# Line 111  C      always need for diagnostics Line 119  C      always need for diagnostics
119        _RL  Diver2(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL  Diver2(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
120        _RL  Diver3(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL  Diver3(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
121        _RL  Diver4(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL  Diver4(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
122          _RL Shannon(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
123          _RL Simpson(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
124    
125        _RL  tmpphy(npmax)        _RL  tmpphy(npmax)
126        _RL  totphy, biotot, maxphy, phymax        _RL  totphy, biotot, maxphy, phymax
# Line 124  C      always need for diagnostics Line 134  C      always need for diagnostics
134        _RL chlup(npmax)        _RL chlup(npmax)
135  #endif  #endif
136  #endif  #endif
137    #ifdef ALLOW_CDOM
138          _RL cdoml  
139          _RL dcdoml
140    #endif
141    
142  #ifdef ALLOW_DIAGNOSTICS  #ifdef ALLOW_DIAGNOSTICS
143  COJ for diagnostics  COJ for diagnostics
# Line 185  c some local variables Line 199  c some local variables
199         _RL PSiupl         _RL PSiupl
200         _RL Tlocal         _RL Tlocal
201         _RL Slocal         _RL Slocal
202           _RL pCO2local
203         _RL Qswlocal         _RL Qswlocal
204         _RL NH4l         _RL NH4l
205         _RL NO2l         _RL NO2l
# Line 277  c Line 292  c
292             Diver2(i,j,k)=0. _d 0             Diver2(i,j,k)=0. _d 0
293             Diver3(i,j,k)=0. _d 0             Diver3(i,j,k)=0. _d 0
294             Diver4(i,j,k)=0. _d 0             Diver4(i,j,k)=0. _d 0
295               Shannon(i,j,k)=0. _d 0
296               Simpson(i,j,k)=1. _d 0
297  #endif  #endif
298    
299  #ifdef ALLOW_DIAGNOSTICS  #ifdef ALLOW_DIAGNOSTICS
# Line 315  COJ Line 332  COJ
332          enddo          enddo
333         ENDDO         ENDDO
334         ENDDO         ENDDO
335    
336    #ifdef DAR_RADTRANS
337           discEs = 0.
338           discEu = 0.
339    #endif
340  c  c
341  c bio-chemical time loop  c bio-chemical time loop
342  c--------------------------------------------------  c--------------------------------------------------
# Line 413  C ========================== i,j loops = Line 435  C ========================== i,j loops =
435  c ------------ these are convenient ------------------------------------  c ------------ these are convenient ------------------------------------
436           DO k=1,Nr           DO k=1,Nr
437            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)
438    #ifdef ALLOW_CDOM
439              cdom_k(k)  = max(Ptr(i,j,k,bi,bj,iCDOM),0. _d 0)
440    #endif
441            DO np = 1,npmax            DO np = 1,npmax
442              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)
443  #ifdef GEIDER  #ifdef GEIDER
# Line 428  c ------------ these are convenient ---- Line 453  c ------------ these are convenient ----
453  c ------------ GET CDOM_k FOR WAVEBANDS_3D and RADTRANS ----------------  c ------------ GET CDOM_k FOR WAVEBANDS_3D and RADTRANS ----------------
454  #ifdef WAVEBANDS  #ifdef WAVEBANDS
455  #if defined(DAR_CALC_ACDOM) || defined(DAR_RADTRANS)  #if defined(DAR_CALC_ACDOM) || defined(DAR_RADTRANS)
456    #ifdef ALLOW_CDOM
457            call MONOD_ACDOM(cdom_k,
458         O                     acdom_k,
459         I                     myThid)
460    #else
461           call MONOD_ACDOM(phychl_k,aphy_chl,aw,           call MONOD_ACDOM(phychl_k,aphy_chl,aw,
462       O                     acdom_k,       O                     acdom_k,
463       I                     myThid)       I                     myThid)
464    #endif
465  #else  #else
466           DO k=1,Nr           DO k=1,Nr
467            DO ilam = 1,tlam            DO ilam = 1,tlam
# Line 601  c           add water and CDOM Line 632  c           add water and CDOM
632              bt_k(k,ilam) = bw(ilam) + bctot + bpart_k(k,ilam)              bt_k(k,ilam) = bw(ilam) + bctot + bpart_k(k,ilam)
633              bb_k(k,ilam) = darwin_bbw*bw(ilam)+bbctot+bbpart_k(k,ilam)              bb_k(k,ilam) = darwin_bbw*bw(ilam)+bbctot+bbpart_k(k,ilam)
634              bb_k(k,ilam) = MAX(darwin_bbmin, bb_k(k,ilam))              bb_k(k,ilam) = MAX(darwin_bbmin, bb_k(k,ilam))
635    c           initialize output variables
636                Edz(ilam,k) = 0.0
637                Esz(ilam,k) = 0.0
638                Euz(ilam,k) = 0.0
639                Estop(ilam,k) = 0.0
640                Eutop(ilam,k) = 0.0
641                amp1(ilam,k) = 0.0
642                amp2(ilam,k) = 0.0
643            ENDDO            ENDDO
644           ENDDO           ENDDO
645    
646  #ifdef DAR_RADTRANS_ITERATIVE           IF (darwin_radtrans_niter.GE.0) THEN
647           call MONOD_RADTRANS_ITER(             call MONOD_RADTRANS_ITER(
648       I                dz_k,rmud,Edwsf,Eswsf,a_k,bt_k,bb_k,       I                dz_k,rmud,Edwsf,Eswsf,a_k,bt_k,bb_k,
649       I                darwin_radtrans_kmax,darwin_radtrans_niter,       I                darwin_radtrans_kmax,darwin_radtrans_niter,
650       O                Edz,Esz,Euz,Eutop,       O                Edz,Esz,Euz,Eutop,
651       O                tirrq,tirrwq,       O                tirrq,tirrwq,
652         O                amp1,amp2,
653       I                myThid)       I                myThid)
654  #else           ELSEIF (darwin_radtrans_niter.EQ.-1) THEN
655  c dzlocal ?????  c dzlocal ?????
656           call MONOD_RADTRANS(             call MONOD_RADTRANS(
657       I                drF,rmud,Edwsf,Eswsf,a_k,bt_k,bb_k,       I                drF,rmud,Edwsf,Eswsf,a_k,bt_k,bb_k,
658       O                Edz,Esz,Euz,Eutop,       O                Edz,Esz,Euz,Eutop,
659       O                tirrq,tirrwq,       O                tirrq,tirrwq,
660       I                myThid)       I                myThid)
661             ELSE
662               call MONOD_RADTRANS_DIRECT(
663         I                dz_k,rmud,Edwsf,Eswsf,a_k,bt_k,bb_k,
664         I                darwin_radtrans_kmax,
665         O                Edz,Esz,Euz,Estop,Eutop,
666         O                tirrq,tirrwq,
667         O                amp1,amp2,
668         I                myThid)
669    #ifdef DAR_CHECK_IRR_CONT
670               DO ilam = 1,tlam
671                 discEs=MAX(discEs,ABS(Estop(ilam,1)-Eswsf(ilam)))
672               ENDDO
673               DO k=1,darwin_radtrans_kmax-1
674                DO ilam = 1,tlam
675                 discEs=MAX(discEs,ABS(Estop(ilam,k+1)-Esz(ilam,k)))
676                 discEu=MAX(discEu,ABS(Eutop(ilam,k+1)-Euz(ilam,k)))
677                ENDDO
678               ENDDO
679             ENDIF
680  #endif  #endif
681  c  c
682  c uses chl from prev timestep (as wavebands does)  c uses chl from prev timestep (as wavebands does)
# Line 665  c brute force... Line 724  c brute force...
724               psil  = max(Ptr(i,j,k,bi,bj,iPOSi ),0. _d 0)               psil  = max(Ptr(i,j,k,bi,bj,iPOSi ),0. _d 0)
725               NH4l  = max(Ptr(i,j,k,bi,bj,iNH4  ),0. _d 0)               NH4l  = max(Ptr(i,j,k,bi,bj,iNH4  ),0. _d 0)
726               NO2l  = max(Ptr(i,j,k,bi,bj,iNO2  ),0. _d 0)               NO2l  = max(Ptr(i,j,k,bi,bj,iNO2  ),0. _d 0)
727    #ifdef ALLOW_CDOM
728                 cdoml = max(Ptr(i,j,k,bi,bj,iCDOM  ),0. _d 0)
729    #endif
730  #ifdef ALLOW_CARBON  #ifdef ALLOW_CARBON
731               dicl  = max(Ptr(i,j,k,bi,bj,iDIC  ),0. _d 0)               dicl  = max(Ptr(i,j,k,bi,bj,iDIC  ),0. _d 0)
732               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 735  c ratio of maximum species Line 797  c ratio of maximum species
797                      Diver4(i,j,k)=Diver4(i,j,k)+1. _d 0                      Diver4(i,j,k)=Diver4(i,j,k)+1. _d 0
798                   endif                   endif
799                 enddo                 enddo
800    c           totphy > thresh0
801                endif
802    c Shannon and Simpson indices
803                Shannon(i,j,k) = 0. _d 0
804    c note: minimal valid value is 1, but we set to zero below threshold
805                Simpson(i,j,k) = 0. _d 0
806                if (totphy.gt.shannon_thresh) then
807                  do np=1,npmax
808                   if (Phy(np) .gt. 0. _d 0) then
809                    tmpphy(np) = Phy(np)/totphy
810                    Shannon(i,j,k)=Shannon(i,j,k)+tmpphy(np)*LOG(tmpphy(np))
811                    Simpson(i,j,k)=Simpson(i,j,k)+tmpphy(np)*tmpphy(np)
812                   endif
813                  enddo
814                  Shannon(i,j,k) = -Shannon(i,j,k)
815                  Simpson(i,j,k) = 1./Simpson(i,j,k)
816              endif              endif
817  #endif  #endif
818    
# Line 771  c for explicit sinking of particulate ma Line 849  c for explicit sinking of particulate ma
849               Slocal = salt(i,j,k,bi,bj)               Slocal = salt(i,j,k,bi,bj)
850  #endif  #endif
851    
852    c choice where to get pCO2 from
853    c taking from igsm dic run - fed through Tflux array
854    c               pCO2local=surfaceForcingT(i,j,bi,bj)
855    c or from darwin carbon module
856    #ifdef ALLOW_CARBON
857                   pCO2local=pCO2(i,j,bi,bj)
858    #else
859                   pCO2local=280. _d -6
860    #endif
861    
862               freefu = max(freefe(i,j,k),0. _d 0)               freefu = max(freefe(i,j,k),0. _d 0)
863               if (k.eq.1) then               if (k.eq.1) then
864                 inputFel = inputFe(i,j,bi,bj)                 inputFel = inputFe(i,j,bi,bj)
# Line 815  c set tendencies to 0 Line 903  c set tendencies to 0
903                 dphychl(np)=0. _d 0                 dphychl(np)=0. _d 0
904               enddo               enddo
905  #endif  #endif
906    #ifdef ALLOW_CDOM
907                dcdoml=0. _d 0
908    #endif
909  #ifdef ALLOW_CARBON  #ifdef ALLOW_CARBON
910               ddicl=0. _d 0               ddicl=0. _d 0
911               ddocl=0. _d 0               ddocl=0. _d 0
# Line 884  c ANNA pass extra variables if WAVEBANDS Line 975  c ANNA pass extra variables if WAVEBANDS
975       I                       pofeupl, psiupl,       I                       pofeupl, psiupl,
976       I                       PARl,       I                       PARl,
977       I                       Tlocal, Slocal,       I                       Tlocal, Slocal,
978         I                       pCO2local,
979       I                       freefu, inputFel,       I                       freefu, inputFel,
980       I                       bottom, dzlocal,       I                       bottom, dzlocal,
981       O                       Rstarl, RNstarl,       O                       Rstarl, RNstarl,
# Line 914  c ANNA pass extra variables if WAVEBANDS Line 1006  c ANNA pass extra variables if WAVEBANDS
1006       I                       dphychl,       I                       dphychl,
1007       I                       chlup,       I                       chlup,
1008  #endif  #endif
1009    #ifdef ALLOW_CDOM
1010         O                       dcdoml,
1011         I                       cdoml,
1012    #endif
1013  #ifdef WAVEBANDS  #ifdef WAVEBANDS
1014       I                       PARw_k(1,k),       I                       PARw_k(1,k),
1015  #endif  #endif
# Line 1120  c in darwin_plankton this is stored for Line 1216  c in darwin_plankton this is stored for
1216  #endif  #endif
1217  #endif  #endif
1218           ENDDO           ENDDO
1219    #ifdef ALLOW_CDOM
1220              Ptr(i,j,k,bi,bj,iCDOM ) = Ptr(i,j,k,bi,bj,iCDOM ) +
1221         &                                  dtplankton*dcdoml
1222    #endif
1223  #ifdef ALLOW_CARBON  #ifdef ALLOW_CARBON
1224            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 ) +
1225       &                                  dtplankton*ddicl       &                                  dtplankton*ddicl
# Line 1271  Coj            no Eu at surface (yet) Line 1371  Coj            no Eu at surface (yet)
1371                 Euave(i,j,k,bi,bj,ilam)=Euave(i,j,k,bi,bj,ilam)+                 Euave(i,j,k,bi,bj,ilam)=Euave(i,j,k,bi,bj,ilam)+
1372       &                                 Euz(ilam,k-1)*dtplankton       &                                 Euz(ilam,k-1)*dtplankton
1373                endif                endif
1374                  Estave(i,j,k,bi,bj,ilam)=Estave(i,j,k,bi,bj,ilam)+
1375         &                                 Estop(ilam,k)*dtplankton
1376                Eutave(i,j,k,bi,bj,ilam)=Eutave(i,j,k,bi,bj,ilam)+                Eutave(i,j,k,bi,bj,ilam)=Eutave(i,j,k,bi,bj,ilam)+
1377       &                                 Eutop(ilam,k)*dtplankton       &                                 Eutop(ilam,k)*dtplankton
1378               enddo               enddo
1379  #endif  #endif
1380    #ifdef DAR_DIAG_IRR_AMPS
1381                 do ilam = 1,tlam
1382                   amp1ave(i,j,k,bi,bj,ilam)=amp1ave(i,j,k,bi,bj,ilam)+
1383         &                                 amp1(ilam,k)*dtplankton
1384                   amp2ave(i,j,k,bi,bj,ilam)=amp2ave(i,j,k,bi,bj,ilam)+
1385         &                                 amp2(ilam,k)*dtplankton
1386                 enddo
1387    #endif
1388  #ifdef DAR_DIAG_ABSORP  #ifdef DAR_DIAG_ABSORP
1389               do ilam = 1,tlam               do ilam = 1,tlam
1390                 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 1299  Coj            no Eu at surface (yet) Line 1409  Coj            no Eu at surface (yet)
1409       &                                 bbpart_k(k,ilam)*dtplankton       &                                 bbpart_k(k,ilam)*dtplankton
1410               enddo               enddo
1411  #endif  #endif
1412    #ifdef DAR_RADTRANS
1413                 if (k.eq.1) then
1414                   rmudave(i,j,bi,bj)=rmudave(i,j,bi,bj)+
1415         &                                 rmud*dtplankton
1416                 endif
1417    #endif
1418  #ifdef DAR_DIAG_RSTAR  #ifdef DAR_DIAG_RSTAR
1419               do np=1,npmax               do np=1,npmax
1420                 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 1404  C       reset the other slot for averagi Line 1520  C       reset the other slot for averagi
1520  C itistime  C itistime
1521  #endif  #endif
1522    
1523    #ifdef DAR_CHECK_IRR_CONT
1524           print*,'max Es, Eu discontinuties',discEs,discEu
1525    #endif
1526    
1527  COJ fill diagnostics  COJ fill diagnostics
1528  #ifdef ALLOW_DIAGNOSTICS  #ifdef ALLOW_DIAGNOSTICS
1529         IF ( useDiagnostics ) THEN         IF ( useDiagnostics ) THEN
# Line 1446  c ANNA end TAVE Line 1566  c ANNA end TAVE
1566          WRITE(diagname,'(A8)') 'Diver4  '          WRITE(diagname,'(A8)') 'Diver4  '
1567          CALL DIAGNOSTICS_FILL( Diver4(1-Olx,1-Oly,1), diagname,          CALL DIAGNOSTICS_FILL( Diver4(1-Olx,1-Oly,1), diagname,
1568       &                         0,Nr,2,bi,bj,myThid )       &                         0,Nr,2,bi,bj,myThid )
1569            WRITE(diagname,'(A8)') 'Shannon '
1570            CALL DIAGNOSTICS_FILL( Shannon(1-Olx,1-Oly,1), diagname,
1571         &                         0,Nr,2,bi,bj,myThid )
1572            WRITE(diagname,'(A8)') 'Simpson '
1573            CALL DIAGNOSTICS_FILL( Simpson(1-Olx,1-Oly,1), diagname,
1574         &                         0,Nr,2,bi,bj,myThid )
1575  #endif  #endif
1576  #ifdef ALLOW_DIAZ  #ifdef ALLOW_DIAZ
1577  #ifdef DAR_DIAG_NFIXP  #ifdef DAR_DIAG_NFIXP
Legend:
Removed from v.1.4  
changed lines
  Added in v.1.11
  ViewVC Help
Powered by ViewVC 1.1.22