/[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.6 by stephd, Thu May 31 21:08:25 2012 UTC revision 1.13 by stephd, Tue Oct 23 16:39:32 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 90  c ANNA define variables for wavebands Line 92  c ANNA define variables for wavebands
92         _RL PARw_k(tlam,Nr)         _RL PARw_k(tlam,Nr)
93         _RL PARwup(tlam)         _RL PARwup(tlam)
94         _RL acdom_k(Nr,tlam)         _RL acdom_k(Nr,tlam)
95           _RL Ek_nll(npmax,tlam)
96           _RL EkoverE_nll(npmax,tlam)
97  #ifdef DAR_RADTRANS  #ifdef DAR_RADTRANS
98         integer iday,iyr,imon,isec,lp,wd,mydate(4)         integer iday,iyr,imon,isec,lp,wd,mydate(4)
99         _RL Edwsf(tlam),Eswsf(tlam)         _RL Edwsf(tlam),Eswsf(tlam)
100         _RL Edz(tlam,Nr),Esz(tlam,Nr),Euz(tlam,Nr),Eutop(tlam,Nr)         _RL Edz(tlam,Nr),Esz(tlam,Nr),Euz(tlam,Nr)
101           _RL Estop(tlam,Nr),Eutop(tlam,Nr)
102         _RL tirrq(nr)         _RL tirrq(nr)
103         _RL tirrwq(tlam,nr)         _RL tirrwq(tlam,nr)
104           _RL amp1(tlam,nr), amp2(tlam,nr)
105         _RL solz         _RL solz
106         _RL rmud         _RL rmud
107         _RL actot,bctot,bbctot         _RL actot,bctot,bbctot
108         _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)
109         _RL bt_k(Nr,tlam), bb_k(Nr,tlam)         _RL bt_k(Nr,tlam), bb_k(Nr,tlam)
110           _RL discEs, discEu
111           INTEGER idiscEs,jdiscEs,kdiscEs,ldiscEs
112           INTEGER idiscEu,jdiscEu,kdiscEu,ldiscEu
113  #else  #else
114         _RL PARwdn(tlam)         _RL PARwdn(tlam)
115  #endif  #endif
# Line 124  C      always need for diagnostics Line 133  C      always need for diagnostics
133  #ifdef GEIDER  #ifdef GEIDER
134        _RL phychl(npmax)        _RL phychl(npmax)
135        _RL phychl_k(npmax,Nr)        _RL phychl_k(npmax,Nr)
136          _RL Ekl(npmax)
137          _RL EkoverEl(npmax)
138          _RL chl2cl(npmax)
139  #ifdef DYNAMIC_CHL  #ifdef DYNAMIC_CHL
140        _RL dphychl(npmax)        _RL dphychl(npmax)
141        _RL chlup(npmax)        _RL chlup(npmax)
142          _RL accliml(npmax)
143  #endif  #endif
144  #endif  #endif
145  #ifdef ALLOW_CDOM  #ifdef ALLOW_CDOM
# Line 194  c some local variables Line 207  c some local variables
207         _RL PSiupl         _RL PSiupl
208         _RL Tlocal         _RL Tlocal
209         _RL Slocal         _RL Slocal
210           _RL pCO2local
211         _RL Qswlocal         _RL Qswlocal
212         _RL NH4l         _RL NH4l
213         _RL NO2l         _RL NO2l
# Line 326  COJ Line 340  COJ
340          enddo          enddo
341         ENDDO         ENDDO
342         ENDDO         ENDDO
343    
344    #ifdef DAR_RADTRANS
345           idiscEs = 0
346           jdiscEs = 0
347           kdiscEs = 0
348           ldiscEs = 0
349           idiscEu = 0
350           jdiscEu = 0
351           kdiscEu = 0
352           ldiscEu = 0
353           discEs = 0.
354           discEu = 0.
355    #endif
356  c  c
357  c bio-chemical time loop  c bio-chemical time loop
358  c--------------------------------------------------  c--------------------------------------------------
# Line 621  c           add water and CDOM Line 648  c           add water and CDOM
648              bt_k(k,ilam) = bw(ilam) + bctot + bpart_k(k,ilam)              bt_k(k,ilam) = bw(ilam) + bctot + bpart_k(k,ilam)
649              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)
650              bb_k(k,ilam) = MAX(darwin_bbmin, bb_k(k,ilam))              bb_k(k,ilam) = MAX(darwin_bbmin, bb_k(k,ilam))
651    c           initialize output variables
652                Edz(ilam,k) = 0.0
653                Esz(ilam,k) = 0.0
654                Euz(ilam,k) = 0.0
655                Estop(ilam,k) = 0.0
656                Eutop(ilam,k) = 0.0
657                amp1(ilam,k) = 0.0
658                amp2(ilam,k) = 0.0
659            ENDDO            ENDDO
660           ENDDO           ENDDO
661    
662  #ifdef DAR_RADTRANS_ITERATIVE           IF (darwin_radtrans_niter.GE.0) THEN
663           call MONOD_RADTRANS_ITER(             call MONOD_RADTRANS_ITER(
664       I                dz_k,rmud,Edwsf,Eswsf,a_k,bt_k,bb_k,       I                dz_k,rmud,Edwsf,Eswsf,a_k,bt_k,bb_k,
665       I                darwin_radtrans_kmax,darwin_radtrans_niter,       I                darwin_radtrans_kmax,darwin_radtrans_niter,
666       O                Edz,Esz,Euz,Eutop,       O                Edz,Esz,Euz,Eutop,
667       O                tirrq,tirrwq,       O                tirrq,tirrwq,
668         O                amp1,amp2,
669       I                myThid)       I                myThid)
670  #else           ELSEIF (darwin_radtrans_niter.EQ.-1) THEN
671  c dzlocal ?????  c dzlocal ?????
672           call MONOD_RADTRANS(             call MONOD_RADTRANS(
673       I                drF,rmud,Edwsf,Eswsf,a_k,bt_k,bb_k,       I                drF,rmud,Edwsf,Eswsf,a_k,bt_k,bb_k,
674       O                Edz,Esz,Euz,Eutop,       O                Edz,Esz,Euz,Eutop,
675       O                tirrq,tirrwq,       O                tirrq,tirrwq,
676       I                myThid)       I                myThid)
677             ELSE
678               call MONOD_RADTRANS_DIRECT(
679         I                dz_k,rmud,Edwsf,Eswsf,a_k,bt_k,bb_k,
680         I                darwin_radtrans_kmax,
681         O                Edz,Esz,Euz,Estop,Eutop,
682         O                tirrq,tirrwq,
683         O                amp1,amp2,
684         I                myThid)
685    #ifdef DAR_CHECK_IRR_CONT
686               IF( dz_k(1) .GT. 0.0 )THEN
687               DO ilam = 1,tlam
688               IF(Eswsf(ilam).GE.darwin_radmodThresh .OR.
689         &        Edwsf(ilam).GE.darwin_radmodThresh ) THEN
690                IF(ABS(Estop(ilam,1)-Eswsf(ilam)) .GT. discEs )THEN
691                  discEs = ABS(Estop(ilam,1)-Eswsf(ilam))
692                  idiscEs = i
693                  jdiscEs = j
694                  kdiscEs = 1
695                  ldiscEs = ilam
696                ENDIF
697                DO k=1,darwin_radtrans_kmax-1
698                 IF(ABS(Estop(ilam,k+1)-Esz(ilam,k)) .GT. discEs)THEN
699                  discEs = ABS(Estop(ilam,k+1)-Esz(ilam,k))
700                  idiscEs = i
701                  jdiscEs = j
702                  kdiscEs = k+1
703                  ldiscEs = ilam
704                 ENDIF
705                 IF(ABS(Eutop(ilam,k+1)-Euz(ilam,k)) .GT. discEu)THEN
706                  discEu = ABS(Eutop(ilam,k+1)-Euz(ilam,k))
707                  idiscEu = i
708                  jdiscEu = j
709                  kdiscEu = k+1
710                  ldiscEu = ilam
711                 ENDIF
712                ENDDO
713               ENDIF
714               ENDDO
715               ENDIF
716  #endif  #endif
717             ENDIF
718  c  c
719  c uses chl from prev timestep (as wavebands does)  c uses chl from prev timestep (as wavebands does)
720  c keep like this in case need to consider upwelling irradiance as affecting the grid box above  c keep like this in case need to consider upwelling irradiance as affecting the grid box above
# Line 810  c for explicit sinking of particulate ma Line 886  c for explicit sinking of particulate ma
886               Slocal = salt(i,j,k,bi,bj)               Slocal = salt(i,j,k,bi,bj)
887  #endif  #endif
888    
889    c choice where to get pCO2 from
890    c taking from igsm dic run - fed through Tflux array
891    c               pCO2local=surfaceForcingT(i,j,bi,bj)
892    c or from darwin carbon module
893    #ifdef ALLOW_CARBON
894                   pCO2local=pCO2(i,j,bi,bj)
895    #else
896                   pCO2local=280. _d -6
897    #endif
898    
899               freefu = max(freefe(i,j,k),0. _d 0)               freefu = max(freefe(i,j,k),0. _d 0)
900               if (k.eq.1) then               if (k.eq.1) then
901                 inputFel = inputFe(i,j,bi,bj)                 inputFel = inputFe(i,j,bi,bj)
# Line 884  c set other arguments to zero Line 970  c set other arguments to zero
970                  NfixPl(np)=0. _d 0                  NfixPl(np)=0. _d 0
971  #endif  #endif
972  #endif  #endif
973    #ifdef DAR_DIAG_PARW
974                    chl2cl(np)=0. _d 0
975    #endif
976    #ifdef DAR_DIAG_EK
977                    Ekl(np)=0. _d 0
978                    EkoverEl(np)=0. _d 0
979                    do ilam=1,tlam
980                      Ek_nll(np,ilam)=0. _d 0
981                      EkoverE_nll(np,ilam)=0. _d 0
982                    enddo
983    #endif
984               enddo               enddo
985    
986    
# Line 926  c ANNA pass extra variables if WAVEBANDS Line 1023  c ANNA pass extra variables if WAVEBANDS
1023       I                       pofeupl, psiupl,       I                       pofeupl, psiupl,
1024       I                       PARl,       I                       PARl,
1025       I                       Tlocal, Slocal,       I                       Tlocal, Slocal,
1026         I                       pCO2local,
1027       I                       freefu, inputFel,       I                       freefu, inputFel,
1028       I                       bottom, dzlocal,       I                       bottom, dzlocal,
1029       O                       Rstarl, RNstarl,       O                       Rstarl, RNstarl,
# Line 952  c ANNA pass extra variables if WAVEBANDS Line 1050  c ANNA pass extra variables if WAVEBANDS
1050  #endif  #endif
1051  #ifdef GEIDER  #ifdef GEIDER
1052       O                       phychl,       O                       phychl,
1053    #ifdef DAR_DIAG_EK
1054         I                       Ekl, EkoverEl,
1055    #endif
1056    #ifdef DAR_DIAG_PARW
1057         I                       chl2cl,
1058    #endif
1059  #ifdef DYNAMIC_CHL  #ifdef DYNAMIC_CHL
1060       I                       dphychl,       I                       dphychl,
1061       I                       chlup,       I                       chlup,
1062    #ifdef DAR_DIAG_EK
1063         O                       accliml,
1064    #endif
1065  #endif  #endif
1066  #ifdef ALLOW_CDOM  #ifdef ALLOW_CDOM
1067       O                       dcdoml,       O                       dcdoml,
# Line 962  c ANNA pass extra variables if WAVEBANDS Line 1069  c ANNA pass extra variables if WAVEBANDS
1069  #endif  #endif
1070  #ifdef WAVEBANDS  #ifdef WAVEBANDS
1071       I                       PARw_k(1,k),       I                       PARw_k(1,k),
1072    #ifdef DAR_DIAG_EK
1073         I                       Ek_nll, EkoverE_nll,
1074    #endif
1075  #endif  #endif
1076  #endif  #endif
1077  #ifdef ALLOW_PAR_DAY  #ifdef ALLOW_PAR_DAY
# Line 1300  c    &                       deltaTclock Line 1410  c    &                       deltaTclock
1410       &                           phychl(np)*dtplankton       &                           phychl(np)*dtplankton
1411               enddo               enddo
1412  #endif  #endif
1413    #ifdef DAR_DIAG_PARW
1414                do ilam=1,tlam
1415                   PARwave(i,j,k,bi,bj,ilam)=PARwave(i,j,k,bi,bj,ilam)+
1416         &                           PARw_k(ilam,k)*dtplankton
1417                enddo
1418                do np=1,npmax
1419                  chl2cave(i,j,k,bi,bj,np)=chl2cave(i,j,k,bi,bj,np)+
1420         &                          chl2cl(np)*dtplankton
1421                enddo
1422    #endif
1423  #ifdef DAR_DIAG_ACDOM  #ifdef DAR_DIAG_ACDOM
1424  c            print*,'acdom',k,acdom_k(k,darwin_diag_acdom_ilam)  c            print*,'acdom',k,acdom_k(k,darwin_diag_acdom_ilam)
1425               aCDOMave(i,j,k,bi,bj)=aCDOMave(i,j,k,bi,bj)+               aCDOMave(i,j,k,bi,bj)=aCDOMave(i,j,k,bi,bj)+
# Line 1321  Coj            no Eu at surface (yet) Line 1441  Coj            no Eu at surface (yet)
1441                 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)+
1442       &                                 Euz(ilam,k-1)*dtplankton       &                                 Euz(ilam,k-1)*dtplankton
1443                endif                endif
1444                  Estave(i,j,k,bi,bj,ilam)=Estave(i,j,k,bi,bj,ilam)+
1445         &                                 Estop(ilam,k)*dtplankton
1446                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)+
1447       &                                 Eutop(ilam,k)*dtplankton       &                                 Eutop(ilam,k)*dtplankton
1448               enddo               enddo
1449  #endif  #endif
1450    #ifdef DAR_DIAG_IRR_AMPS
1451                 do ilam = 1,tlam
1452                   amp1ave(i,j,k,bi,bj,ilam)=amp1ave(i,j,k,bi,bj,ilam)+
1453         &                                 amp1(ilam,k)*dtplankton
1454                   amp2ave(i,j,k,bi,bj,ilam)=amp2ave(i,j,k,bi,bj,ilam)+
1455         &                                 amp2(ilam,k)*dtplankton
1456                 enddo
1457    #endif
1458  #ifdef DAR_DIAG_ABSORP  #ifdef DAR_DIAG_ABSORP
1459               do ilam = 1,tlam               do ilam = 1,tlam
1460                 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 1349  Coj            no Eu at surface (yet) Line 1479  Coj            no Eu at surface (yet)
1479       &                                 bbpart_k(k,ilam)*dtplankton       &                                 bbpart_k(k,ilam)*dtplankton
1480               enddo               enddo
1481  #endif  #endif
1482    #ifdef DAR_RADTRANS
1483                 if (k.eq.1) then
1484                   rmudave(i,j,bi,bj)=rmudave(i,j,bi,bj)+
1485         &                                 rmud*dtplankton
1486                 endif
1487    #endif
1488    #ifdef DAR_DIAG_EK
1489                do np=1,npmax
1490                 Ekave(i,j,k,bi,bj,np)=Ekave(i,j,k,bi,bj,np)+
1491         &                        Ekl(np)*dtplankton
1492                 EkoverEave(i,j,k,bi,bj,np)=EkoverEave(i,j,k,bi,bj,np)+
1493         &                        EkoverEl(np)*dtplankton
1494                 acclimave(i,j,k,bi,bj,np)=acclimave(i,j,k,bi,bj,np)+
1495         &                        accliml(np)*dtplankton
1496                 do ilam=1,tlam
1497                    Ek_nlave(i,j,k,bi,bj,np,ilam)=
1498         &                        Ek_nlave(i,j,k,bi,bj,np,ilam)+
1499         &                        Ek_nll(np,ilam)*dtplankton
1500                    EkoverE_nlave(i,j,k,bi,bj,np,ilam)=
1501         &                        EkoverE_nlave(i,j,k,bi,bj,np,ilam)+
1502         &                        EkoverE_nll(np,ilam)*dtplankton
1503                 enddo
1504                enddo
1505    #endif
1506  #ifdef DAR_DIAG_RSTAR  #ifdef DAR_DIAG_RSTAR
1507               do np=1,npmax               do np=1,npmax
1508                 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 1454  C       reset the other slot for averagi Line 1608  C       reset the other slot for averagi
1608  C itistime  C itistime
1609  #endif  #endif
1610    
1611    #ifdef DAR_CHECK_IRR_CONT
1612           i = myXGlobalLo-1+(bi-1)*sNx+idiscEs
1613           j = myYGlobalLo-1+(bj-1)*sNy+jdiscEs
1614           write(6,'(I4.4,X,A,4(X,I4),1PE24.16)')myProcId,'max Es disc',
1615         &                                   i,j,kdiscEs,ldiscEs,discEs
1616           i = myXGlobalLo-1+(bi-1)*sNx+idiscEu
1617           j = myYGlobalLo-1+(bj-1)*sNy+jdiscEu
1618           write(6,'(I4.4,X,A,4(X,I4),1PE24.16)')myProcId,'max Eu disc',
1619         &                                   i,j,kdiscEu,ldiscEu,discEu
1620    #endif
1621    
1622  COJ fill diagnostics  COJ fill diagnostics
1623  #ifdef ALLOW_DIAGNOSTICS  #ifdef ALLOW_DIAGNOSTICS
1624         IF ( useDiagnostics ) THEN         IF ( useDiagnostics ) THEN
Legend:
Removed from v.1.6  
changed lines
  Added in v.1.13
  ViewVC Help
Powered by ViewVC 1.1.22