/[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.3 by jahn, Tue Oct 11 18:10:55 2011 UTC revision 1.10 by jahn, Thu Aug 23 21:48:24 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
# Line 111  C      always need for diagnostics Line 118  C      always need for diagnostics
118        _RL  Diver2(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL  Diver2(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
119        _RL  Diver3(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL  Diver3(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
120        _RL  Diver4(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL  Diver4(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
121          _RL Shannon(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
122          _RL Simpson(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
123    
124        _RL  tmpphy(npmax)        _RL  tmpphy(npmax)
125        _RL  totphy, biotot, maxphy, phymax        _RL  totphy, biotot, maxphy, phymax
# Line 124  C      always need for diagnostics Line 133  C      always need for diagnostics
133        _RL chlup(npmax)        _RL chlup(npmax)
134  #endif  #endif
135  #endif  #endif
136    #ifdef ALLOW_CDOM
137          _RL cdoml  
138          _RL dcdoml
139    #endif
140    
141  #ifdef ALLOW_DIAGNOSTICS  #ifdef ALLOW_DIAGNOSTICS
142  COJ for diagnostics  COJ for diagnostics
# Line 185  c some local variables Line 198  c some local variables
198         _RL PSiupl         _RL PSiupl
199         _RL Tlocal         _RL Tlocal
200         _RL Slocal         _RL Slocal
201           _RL pCO2local
202         _RL Qswlocal         _RL Qswlocal
203         _RL NH4l         _RL NH4l
204         _RL NO2l         _RL NO2l
# Line 277  c Line 291  c
291             Diver2(i,j,k)=0. _d 0             Diver2(i,j,k)=0. _d 0
292             Diver3(i,j,k)=0. _d 0             Diver3(i,j,k)=0. _d 0
293             Diver4(i,j,k)=0. _d 0             Diver4(i,j,k)=0. _d 0
294               Shannon(i,j,k)=0. _d 0
295               Simpson(i,j,k)=1. _d 0
296  #endif  #endif
297    
298  #ifdef ALLOW_DIAGNOSTICS  #ifdef ALLOW_DIAGNOSTICS
299  COJ for diagnostics  COJ for diagnostics
300             PParr(i,j,k) = 0. _d 0             PParr(i,j,k) = 0. _d 0
301             Nfixarr(i,j,k) = 0. _d 0             Nfixarr(i,j,k) = 0. _d 0
302    #ifdef DAR_DIAG_CHL
303               GeiderChlarr(i,j,k) = 0. _d 0
304               GeiderChl2Carr(i,j,k) = 0. _d 0
305               DoneyChlarr(i,j,k) = 0. _d 0
306               DoneyChl2Carr(i,j,k) = 0. _d 0
307               CloernChlarr(i,j,k) = 0. _d 0
308               CloernChl2Carr(i,j,k) = 0. _d 0
309    #endif
310  c ANNA_TAVE  c ANNA_TAVE
311  #ifdef WAVES_DIAG_PCHL  #ifdef WAVES_DIAG_PCHL
312             DO np=1,npmax             DO np=1,npmax
# Line 405  C ========================== i,j loops = Line 429  C ========================== i,j loops =
429  c ------------ these are convenient ------------------------------------  c ------------ these are convenient ------------------------------------
430           DO k=1,Nr           DO k=1,Nr
431            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)
432    #ifdef ALLOW_CDOM
433              cdom_k(k)  = max(Ptr(i,j,k,bi,bj,iCDOM),0. _d 0)
434    #endif
435            DO np = 1,npmax            DO np = 1,npmax
436              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)
437  #ifdef GEIDER  #ifdef GEIDER
# Line 420  c ------------ these are convenient ---- Line 447  c ------------ these are convenient ----
447  c ------------ GET CDOM_k FOR WAVEBANDS_3D and RADTRANS ----------------  c ------------ GET CDOM_k FOR WAVEBANDS_3D and RADTRANS ----------------
448  #ifdef WAVEBANDS  #ifdef WAVEBANDS
449  #if defined(DAR_CALC_ACDOM) || defined(DAR_RADTRANS)  #if defined(DAR_CALC_ACDOM) || defined(DAR_RADTRANS)
450    #ifdef ALLOW_CDOM
451            call MONOD_ACDOM(cdom_k,
452         O                     acdom_k,
453         I                     myThid)
454    #else
455           call MONOD_ACDOM(phychl_k,aphy_chl,aw,           call MONOD_ACDOM(phychl_k,aphy_chl,aw,
456       O                     acdom_k,       O                     acdom_k,
457       I                     myThid)       I                     myThid)
458    #endif
459  #else  #else
460           DO k=1,Nr           DO k=1,Nr
461            DO ilam = 1,tlam            DO ilam = 1,tlam
# Line 593  c           add water and CDOM Line 626  c           add water and CDOM
626              bt_k(k,ilam) = bw(ilam) + bctot + bpart_k(k,ilam)              bt_k(k,ilam) = bw(ilam) + bctot + bpart_k(k,ilam)
627              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)
628              bb_k(k,ilam) = MAX(darwin_bbmin, bb_k(k,ilam))              bb_k(k,ilam) = MAX(darwin_bbmin, bb_k(k,ilam))
629    c           initialize output variables
630                Edz(ilam,k) = 0.0
631                Esz(ilam,k) = 0.0
632                Euz(ilam,k) = 0.0
633                Estop(ilam,k) = 0.0
634                Eutop(ilam,k) = 0.0
635                amp1(ilam,k) = 0.0
636                amp2(ilam,k) = 0.0
637            ENDDO            ENDDO
638           ENDDO           ENDDO
639    
# Line 602  c           add water and CDOM Line 643  c           add water and CDOM
643       I                darwin_radtrans_kmax,darwin_radtrans_niter,       I                darwin_radtrans_kmax,darwin_radtrans_niter,
644       O                Edz,Esz,Euz,Eutop,       O                Edz,Esz,Euz,Eutop,
645       O                tirrq,tirrwq,       O                tirrq,tirrwq,
646         O                amp1,amp2,
647       I                myThid)       I                myThid)
648  #else  #else
649  c dzlocal ?????  c dzlocal ?????
# Line 657  c brute force... Line 699  c brute force...
699               psil  = max(Ptr(i,j,k,bi,bj,iPOSi ),0. _d 0)               psil  = max(Ptr(i,j,k,bi,bj,iPOSi ),0. _d 0)
700               NH4l  = max(Ptr(i,j,k,bi,bj,iNH4  ),0. _d 0)               NH4l  = max(Ptr(i,j,k,bi,bj,iNH4  ),0. _d 0)
701               NO2l  = max(Ptr(i,j,k,bi,bj,iNO2  ),0. _d 0)               NO2l  = max(Ptr(i,j,k,bi,bj,iNO2  ),0. _d 0)
702    #ifdef ALLOW_CDOM
703                 cdoml = max(Ptr(i,j,k,bi,bj,iCDOM  ),0. _d 0)
704    #endif
705  #ifdef ALLOW_CARBON  #ifdef ALLOW_CARBON
706               dicl  = max(Ptr(i,j,k,bi,bj,iDIC  ),0. _d 0)               dicl  = max(Ptr(i,j,k,bi,bj,iDIC  ),0. _d 0)
707               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 772  c ratio of maximum species
772                      Diver4(i,j,k)=Diver4(i,j,k)+1. _d 0                      Diver4(i,j,k)=Diver4(i,j,k)+1. _d 0
773                   endif                   endif
774                 enddo                 enddo
775    c           totphy > thresh0
776                endif
777    c Shannon and Simpson indices
778                Shannon(i,j,k) = 0. _d 0
779    c note: minimal valid value is 1, but we set to zero below threshold
780                Simpson(i,j,k) = 0. _d 0
781                if (totphy.gt.shannon_thresh) then
782                  do np=1,npmax
783                   if (Phy(np) .gt. 0. _d 0) then
784                    tmpphy(np) = Phy(np)/totphy
785                    Shannon(i,j,k)=Shannon(i,j,k)+tmpphy(np)*LOG(tmpphy(np))
786                    Simpson(i,j,k)=Simpson(i,j,k)+tmpphy(np)*tmpphy(np)
787                   endif
788                  enddo
789                  Shannon(i,j,k) = -Shannon(i,j,k)
790                  Simpson(i,j,k) = 1./Simpson(i,j,k)
791              endif              endif
792  #endif  #endif
793    
# Line 763  c for explicit sinking of particulate ma Line 824  c for explicit sinking of particulate ma
824               Slocal = salt(i,j,k,bi,bj)               Slocal = salt(i,j,k,bi,bj)
825  #endif  #endif
826    
827    c choice where to get pCO2 from
828    c taking from igsm dic run - fed through Tflux array
829    c               pCO2local=surfaceForcingT(i,j,bi,bj)
830    c or from darwin carbon module
831    #ifdef ALLOW_CARBON
832                   pCO2local=pCO2(i,j,bi,bj)
833    #else
834                   pCO2local=280. _d -6
835    #endif
836    
837               freefu = max(freefe(i,j,k),0. _d 0)               freefu = max(freefe(i,j,k),0. _d 0)
838               if (k.eq.1) then               if (k.eq.1) then
839                 inputFel = inputFe(i,j,bi,bj)                 inputFel = inputFe(i,j,bi,bj)
# Line 807  c set tendencies to 0 Line 878  c set tendencies to 0
878                 dphychl(np)=0. _d 0                 dphychl(np)=0. _d 0
879               enddo               enddo
880  #endif  #endif
881    #ifdef ALLOW_CDOM
882                dcdoml=0. _d 0
883    #endif
884  #ifdef ALLOW_CARBON  #ifdef ALLOW_CARBON
885               ddicl=0. _d 0               ddicl=0. _d 0
886               ddocl=0. _d 0               ddocl=0. _d 0
# Line 876  c ANNA pass extra variables if WAVEBANDS Line 950  c ANNA pass extra variables if WAVEBANDS
950       I                       pofeupl, psiupl,       I                       pofeupl, psiupl,
951       I                       PARl,       I                       PARl,
952       I                       Tlocal, Slocal,       I                       Tlocal, Slocal,
953         I                       pCO2local,
954       I                       freefu, inputFel,       I                       freefu, inputFel,
955       I                       bottom, dzlocal,       I                       bottom, dzlocal,
956       O                       Rstarl, RNstarl,       O                       Rstarl, RNstarl,
# Line 906  c ANNA pass extra variables if WAVEBANDS Line 981  c ANNA pass extra variables if WAVEBANDS
981       I                       dphychl,       I                       dphychl,
982       I                       chlup,       I                       chlup,
983  #endif  #endif
984    #ifdef ALLOW_CDOM
985         O                       dcdoml,
986         I                       cdoml,
987    #endif
988  #ifdef WAVEBANDS  #ifdef WAVEBANDS
989       I                       PARw_k(1,k),       I                       PARw_k(1,k),
990  #endif  #endif
# Line 1112  c in darwin_plankton this is stored for Line 1191  c in darwin_plankton this is stored for
1191  #endif  #endif
1192  #endif  #endif
1193           ENDDO           ENDDO
1194    #ifdef ALLOW_CDOM
1195              Ptr(i,j,k,bi,bj,iCDOM ) = Ptr(i,j,k,bi,bj,iCDOM ) +
1196         &                                  dtplankton*dcdoml
1197    #endif
1198  #ifdef ALLOW_CARBON  #ifdef ALLOW_CARBON
1199            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 ) +
1200       &                                  dtplankton*ddicl       &                                  dtplankton*ddicl
# Line 1263  Coj            no Eu at surface (yet) Line 1346  Coj            no Eu at surface (yet)
1346                 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)+
1347       &                                 Euz(ilam,k-1)*dtplankton       &                                 Euz(ilam,k-1)*dtplankton
1348                endif                endif
1349                  Estave(i,j,k,bi,bj,ilam)=Estave(i,j,k,bi,bj,ilam)+
1350         &                                 Estop(ilam,k)*dtplankton
1351                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)+
1352       &                                 Eutop(ilam,k)*dtplankton       &                                 Eutop(ilam,k)*dtplankton
1353               enddo               enddo
1354  #endif  #endif
1355    #ifdef DAR_DIAG_IRR_AMPS
1356                 do ilam = 1,tlam
1357                   amp1ave(i,j,k,bi,bj,ilam)=amp1ave(i,j,k,bi,bj,ilam)+
1358         &                                 amp1(ilam,k)*dtplankton
1359                   amp2ave(i,j,k,bi,bj,ilam)=amp2ave(i,j,k,bi,bj,ilam)+
1360         &                                 amp2(ilam,k)*dtplankton
1361                 enddo
1362    #endif
1363  #ifdef DAR_DIAG_ABSORP  #ifdef DAR_DIAG_ABSORP
1364               do ilam = 1,tlam               do ilam = 1,tlam
1365                 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 1384  Coj            no Eu at surface (yet)
1384       &                                 bbpart_k(k,ilam)*dtplankton       &                                 bbpart_k(k,ilam)*dtplankton
1385               enddo               enddo
1386  #endif  #endif
1387    #ifdef DAR_RADTRANS
1388                 if (k.eq.1) then
1389                   rmudave(i,j,bi,bj)=rmudave(i,j,bi,bj)+
1390         &                                 rmud*dtplankton
1391                 endif
1392    #endif
1393  #ifdef DAR_DIAG_RSTAR  #ifdef DAR_DIAG_RSTAR
1394               do np=1,npmax               do np=1,npmax
1395                 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 1438  c ANNA end TAVE Line 1537  c ANNA end TAVE
1537          WRITE(diagname,'(A8)') 'Diver4  '          WRITE(diagname,'(A8)') 'Diver4  '
1538          CALL DIAGNOSTICS_FILL( Diver4(1-Olx,1-Oly,1), diagname,          CALL DIAGNOSTICS_FILL( Diver4(1-Olx,1-Oly,1), diagname,
1539       &                         0,Nr,2,bi,bj,myThid )       &                         0,Nr,2,bi,bj,myThid )
1540            WRITE(diagname,'(A8)') 'Shannon '
1541            CALL DIAGNOSTICS_FILL( Shannon(1-Olx,1-Oly,1), diagname,
1542         &                         0,Nr,2,bi,bj,myThid )
1543            WRITE(diagname,'(A8)') 'Simpson '
1544            CALL DIAGNOSTICS_FILL( Simpson(1-Olx,1-Oly,1), diagname,
1545         &                         0,Nr,2,bi,bj,myThid )
1546  #endif  #endif
1547  #ifdef ALLOW_DIAZ  #ifdef ALLOW_DIAZ
1548  #ifdef DAR_DIAG_NFIXP  #ifdef DAR_DIAG_NFIXP
Legend:
Removed from v.1.3  
changed lines
  Added in v.1.10
  ViewVC Help
Powered by ViewVC 1.1.22