/[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.18 by jahn, Fri Jul 24 21:47:26 2015 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 87  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 111  C      always need for diagnostics Line 123  C      always need for diagnostics
123        _RL  Diver2(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL  Diver2(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
124        _RL  Diver3(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL  Diver3(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
125        _RL  Diver4(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL  Diver4(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
126          _RL Shannon(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
127          _RL Simpson(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
128    
129        _RL  tmpphy(npmax)        _RL  tmpphy(npmax)
130        _RL  totphy, biotot, maxphy, phymax        _RL  totphy, biotot, maxphy, phymax
# Line 119  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
144  #endif  #endif
145    #ifdef ALLOW_CDOM
146          _RL cdoml  
147          _RL dcdoml
148  #endif  #endif
149    
150  #ifdef ALLOW_DIAGNOSTICS  #ifdef ALLOW_DIAGNOSTICS
# Line 185  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 277  c Line 300  c
300             Diver2(i,j,k)=0. _d 0             Diver2(i,j,k)=0. _d 0
301             Diver3(i,j,k)=0. _d 0             Diver3(i,j,k)=0. _d 0
302             Diver4(i,j,k)=0. _d 0             Diver4(i,j,k)=0. _d 0
303               Shannon(i,j,k)=0. _d 0
304               Simpson(i,j,k)=1. _d 0
305  #endif  #endif
306    
307  #ifdef ALLOW_DIAGNOSTICS  #ifdef ALLOW_DIAGNOSTICS
# Line 315  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 384  C        always fill; this will be the s Line 422  C        always fill; this will be the s
422  C        way it won't blow up for weird diagnostics periods.  C        way it won't blow up for weird diagnostics periods.
423  C        we fill before updating, so the diag is the one used in this time  C        we fill before updating, so the diag is the one used in this time
424  C        step  C        step
425           CALL DIAGNOSTICS_FILL(           IF ( useDiagnostics ) THEN
426              CALL DIAGNOSTICS_FILL(
427       &         PARday(1-Olx,1-Oly,1,bi,bj,PARiprev),'PARday  ',       &         PARday(1-Olx,1-Oly,1,bi,bj,PARiprev),'PARday  ',
428       &         0,Nr,2,bi,bj,myThid )       &         0,Nr,2,bi,bj,myThid )
429             ENDIF
430  #endif  #endif
431  #endif /* ALLOW_PAR_DAY */  #endif /* ALLOW_PAR_DAY */
432    
# Line 413  C ========================== i,j loops = Line 453  C ========================== i,j loops =
453  c ------------ these are convenient ------------------------------------  c ------------ these are convenient ------------------------------------
454           DO k=1,Nr           DO k=1,Nr
455            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)
456    #ifdef ALLOW_CDOM
457              cdom_k(k)  = max(Ptr(i,j,k,bi,bj,iCDOM),0. _d 0)
458    #endif
459            DO np = 1,npmax            DO np = 1,npmax
460              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)
461  #ifdef GEIDER  #ifdef GEIDER
# Line 428  c ------------ these are convenient ---- Line 471  c ------------ these are convenient ----
471  c ------------ GET CDOM_k FOR WAVEBANDS_3D and RADTRANS ----------------  c ------------ GET CDOM_k FOR WAVEBANDS_3D and RADTRANS ----------------
472  #ifdef WAVEBANDS  #ifdef WAVEBANDS
473  #if defined(DAR_CALC_ACDOM) || defined(DAR_RADTRANS)  #if defined(DAR_CALC_ACDOM) || defined(DAR_RADTRANS)
474    #ifdef ALLOW_CDOM
475            call MONOD_ACDOM(cdom_k,
476         O                     acdom_k,
477         I                     myThid)
478    #else
479           call MONOD_ACDOM(phychl_k,aphy_chl,aw,           call MONOD_ACDOM(phychl_k,aphy_chl,aw,
480       O                     acdom_k,       O                     acdom_k,
481       I                     myThid)       I                     myThid)
482    #endif
483  #else  #else
484           DO k=1,Nr           DO k=1,Nr
485            DO ilam = 1,tlam            DO ilam = 1,tlam
# Line 458  c ------------ GET INCIDENT NON-SPECTRAL Line 507  c ------------ GET INCIDENT NON-SPECTRAL
507    
508  #else /* not USE_QSW */  #else /* not USE_QSW */
509    
510           lite = sfac(j)*maskC(i,j,1,bi,bj)/86400*1 _d 6  C        convert W/m2 to uEin/s/m2
511             lite = sfac(j)*parconv*maskC(i,j,1,bi,bj)
512    
513  #endif /* not USE_QSW */  #endif /* not USE_QSW */
514  #endif /* not READ_PAR */  #endif /* not READ_PAR */
# Line 601  c           add water and CDOM Line 651  c           add water and CDOM
651              bt_k(k,ilam) = bw(ilam) + bctot + bpart_k(k,ilam)              bt_k(k,ilam) = bw(ilam) + bctot + bpart_k(k,ilam)
652              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)
653              bb_k(k,ilam) = MAX(darwin_bbmin, bb_k(k,ilam))              bb_k(k,ilam) = MAX(darwin_bbmin, bb_k(k,ilam))
654    c           initialize output variables
655                Edz(ilam,k) = 0.0
656                Esz(ilam,k) = 0.0
657                Euz(ilam,k) = 0.0
658                Estop(ilam,k) = 0.0
659                Eutop(ilam,k) = 0.0
660                amp1(ilam,k) = 0.0
661                amp2(ilam,k) = 0.0
662            ENDDO            ENDDO
663           ENDDO           ENDDO
664    
665  #ifdef DAR_RADTRANS_ITERATIVE           IF (darwin_radtrans_niter.GE.0) THEN
666           call MONOD_RADTRANS_ITER(             call MONOD_RADTRANS_ITER(
667       I                dz_k,rmud,Edwsf,Eswsf,a_k,bt_k,bb_k,       I                dz_k,rmud,Edwsf,Eswsf,a_k,bt_k,bb_k,
668       I                darwin_radtrans_kmax,darwin_radtrans_niter,       I                darwin_radtrans_kmax,darwin_radtrans_niter,
669       O                Edz,Esz,Euz,Eutop,       O                Edz,Esz,Euz,Eutop,
670       O                tirrq,tirrwq,       O                tirrq,tirrwq,
671         O                amp1,amp2,
672       I                myThid)       I                myThid)
673  #else           ELSEIF (darwin_radtrans_niter.EQ.-1) THEN
674  c dzlocal ?????  c dzlocal ?????
675           call MONOD_RADTRANS(             call MONOD_RADTRANS(
676       I                drF,rmud,Edwsf,Eswsf,a_k,bt_k,bb_k,       I                drF,rmud,Edwsf,Eswsf,a_k,bt_k,bb_k,
677       O                Edz,Esz,Euz,Eutop,       O                Edz,Esz,Euz,Eutop,
678       O                tirrq,tirrwq,       O                tirrq,tirrwq,
679       I                myThid)       I                myThid)
680             ELSE
681               call MONOD_RADTRANS_DIRECT(
682         I                dz_k,rmud,Edwsf,Eswsf,a_k,bt_k,bb_k,
683         I                darwin_radtrans_kmax,
684         O                Edz,Esz,Euz,Estop,Eutop,
685         O                tirrq,tirrwq,
686         O                amp1,amp2,
687         I                myThid)
688    #ifdef DAR_CHECK_IRR_CONT
689               IF( dz_k(1) .GT. 0.0 )THEN
690               DO ilam = 1,tlam
691               IF(Eswsf(ilam).GE.darwin_radmodThresh .OR.
692         &        Edwsf(ilam).GE.darwin_radmodThresh ) THEN
693                IF(ABS(Estop(ilam,1)-Eswsf(ilam)) .GT. discEs )THEN
694                  discEs = ABS(Estop(ilam,1)-Eswsf(ilam))
695                  idiscEs = i
696                  jdiscEs = j
697                  kdiscEs = 1
698                  ldiscEs = ilam
699                ENDIF
700                DO k=1,darwin_radtrans_kmax-1
701                 IF(ABS(Estop(ilam,k+1)-Esz(ilam,k)) .GT. discEs)THEN
702                  discEs = ABS(Estop(ilam,k+1)-Esz(ilam,k))
703                  idiscEs = i
704                  jdiscEs = j
705                  kdiscEs = k+1
706                  ldiscEs = ilam
707                 ENDIF
708                 IF(ABS(Eutop(ilam,k+1)-Euz(ilam,k)) .GT. discEu)THEN
709                  discEu = ABS(Eutop(ilam,k+1)-Euz(ilam,k))
710                  idiscEu = i
711                  jdiscEu = j
712                  kdiscEu = k+1
713                  ldiscEu = ilam
714                 ENDIF
715                ENDDO
716               ENDIF
717               ENDDO
718               ENDIF
719  #endif  #endif
720             ENDIF
721  c  c
722  c uses chl from prev timestep (as wavebands does)  c uses chl from prev timestep (as wavebands does)
723  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 665  c brute force... Line 764  c brute force...
764               psil  = max(Ptr(i,j,k,bi,bj,iPOSi ),0. _d 0)               psil  = max(Ptr(i,j,k,bi,bj,iPOSi ),0. _d 0)
765               NH4l  = max(Ptr(i,j,k,bi,bj,iNH4  ),0. _d 0)               NH4l  = max(Ptr(i,j,k,bi,bj,iNH4  ),0. _d 0)
766               NO2l  = max(Ptr(i,j,k,bi,bj,iNO2  ),0. _d 0)               NO2l  = max(Ptr(i,j,k,bi,bj,iNO2  ),0. _d 0)
767    #ifdef ALLOW_CDOM
768                 cdoml = max(Ptr(i,j,k,bi,bj,iCDOM  ),0. _d 0)
769    #endif
770  #ifdef ALLOW_CARBON  #ifdef ALLOW_CARBON
771               dicl  = max(Ptr(i,j,k,bi,bj,iDIC  ),0. _d 0)               dicl  = max(Ptr(i,j,k,bi,bj,iDIC  ),0. _d 0)
772               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 837  c ratio of maximum species
837                      Diver4(i,j,k)=Diver4(i,j,k)+1. _d 0                      Diver4(i,j,k)=Diver4(i,j,k)+1. _d 0
838                   endif                   endif
839                 enddo                 enddo
840    c           totphy > thresh0
841                endif
842    c Shannon and Simpson indices
843                Shannon(i,j,k) = 0. _d 0
844    c note: minimal valid value is 1, but we set to zero below threshold
845                Simpson(i,j,k) = 0. _d 0
846                if (totphy.gt.shannon_thresh) then
847                  do np=1,npmax
848                   if (Phy(np) .gt. 0. _d 0) then
849                    tmpphy(np) = Phy(np)/totphy
850                    Shannon(i,j,k)=Shannon(i,j,k)+tmpphy(np)*LOG(tmpphy(np))
851                    Simpson(i,j,k)=Simpson(i,j,k)+tmpphy(np)*tmpphy(np)
852                   endif
853                  enddo
854                  Shannon(i,j,k) = -Shannon(i,j,k)
855                  Simpson(i,j,k) = 1./Simpson(i,j,k)
856              endif              endif
857  #endif  #endif
858    
# Line 771  c for explicit sinking of particulate ma Line 889  c for explicit sinking of particulate ma
889               Slocal = salt(i,j,k,bi,bj)               Slocal = salt(i,j,k,bi,bj)
890  #endif  #endif
891    
892    c choice where to get pCO2 from
893    c taking from igsm dic run - fed through Tflux array
894    c               pCO2local=surfaceForcingT(i,j,bi,bj)
895    c or from darwin carbon module
896    #ifdef ALLOW_CARBON
897    #ifdef pH_3D
898                   pCO2local=pCO2(i,j,k,bi,bj)
899    #else
900                   pCO2local=pCO2(i,j,bi,bj)
901    #endif
902    #else
903                   pCO2local=280. _d -6
904    #endif
905    
906               freefu = max(freefe(i,j,k),0. _d 0)               freefu = max(freefe(i,j,k),0. _d 0)
907               if (k.eq.1) then               if (k.eq.1) then
908                 inputFel = inputFe(i,j,bi,bj)                 inputFel = inputFe(i,j,bi,bj)
# Line 815  c set tendencies to 0 Line 947  c set tendencies to 0
947                 dphychl(np)=0. _d 0                 dphychl(np)=0. _d 0
948               enddo               enddo
949  #endif  #endif
950    #ifdef ALLOW_CDOM
951                dcdoml=0. _d 0
952    #endif
953  #ifdef ALLOW_CARBON  #ifdef ALLOW_CARBON
954               ddicl=0. _d 0               ddicl=0. _d 0
955               ddocl=0. _d 0               ddocl=0. _d 0
# Line 842  c set other arguments to zero Line 977  c set other arguments to zero
977                  NfixPl(np)=0. _d 0                  NfixPl(np)=0. _d 0
978  #endif  #endif
979  #endif  #endif
980    #ifdef DAR_DIAG_PARW
981                    chl2cl(np)=0. _d 0
982    #endif
983    #ifdef DAR_DIAG_EK
984                    Ekl(np)=0. _d 0
985                    EkoverEl(np)=0. _d 0
986                    do ilam=1,tlam
987                      Ek_nll(np,ilam)=0. _d 0
988                      EkoverE_nll(np,ilam)=0. _d 0
989                    enddo
990    #endif
991               enddo               enddo
992    
993    
# Line 884  c ANNA pass extra variables if WAVEBANDS Line 1030  c ANNA pass extra variables if WAVEBANDS
1030       I                       pofeupl, psiupl,       I                       pofeupl, psiupl,
1031       I                       PARl,       I                       PARl,
1032       I                       Tlocal, Slocal,       I                       Tlocal, Slocal,
1033         I                       pCO2local,
1034       I                       freefu, inputFel,       I                       freefu, inputFel,
1035       I                       bottom, dzlocal,       I                       bottom, dzlocal,
1036       O                       Rstarl, RNstarl,       O                       Rstarl, RNstarl,
# Line 910  c ANNA pass extra variables if WAVEBANDS Line 1057  c ANNA pass extra variables if WAVEBANDS
1057  #endif  #endif
1058  #ifdef GEIDER  #ifdef GEIDER
1059       O                       phychl,       O                       phychl,
1060    #ifdef DAR_DIAG_EK
1061         I                       Ekl, EkoverEl,
1062    #endif
1063    #ifdef DAR_DIAG_PARW
1064         I                       chl2cl,
1065    #endif
1066  #ifdef DYNAMIC_CHL  #ifdef DYNAMIC_CHL
1067       I                       dphychl,       I                       dphychl,
1068       I                       chlup,       I                       chlup,
1069    #ifdef DAR_DIAG_EK
1070         O                       accliml,
1071    #endif
1072    #endif
1073    #ifdef ALLOW_CDOM
1074         O                       dcdoml,
1075         I                       cdoml,
1076  #endif  #endif
1077  #ifdef WAVEBANDS  #ifdef WAVEBANDS
1078       I                       PARw_k(1,k),       I                       PARw_k(1,k),
1079    #ifdef DAR_DIAG_EK
1080         I                       Ek_nll, EkoverE_nll,
1081    #endif
1082  #endif  #endif
1083  #endif  #endif
1084  #ifdef ALLOW_PAR_DAY  #ifdef ALLOW_PAR_DAY
# Line 1059  c ---- end steph's alternative Line 1222  c ---- end steph's alternative
1222                  dfetl=dfetl+fet_flx(i,j,k,bi,bj)                  dfetl=dfetl+fet_flx(i,j,k,bi,bj)
1223                  dsil=dsil+si_flx(i,j,k,bi,bj)                  dsil=dsil+si_flx(i,j,k,bi,bj)
1224  #endif  #endif
1225  c  
1226    #ifdef ALLOW_OBCS
1227              IF (useOBCS) THEN
1228                dpo4l  = dpo4l *maskInC(i,j,bi,bj)
1229                dno3l  = dno3l *maskInC(i,j,bi,bj)
1230                dfetl  = dfetl *maskInC(i,j,bi,bj)
1231                dsil   = dsil  *maskInC(i,j,bi,bj)
1232                ddopl  = ddopl *maskInC(i,j,bi,bj)
1233                ddonl  = ddonl *maskInC(i,j,bi,bj)
1234                ddofel = ddofel*maskInC(i,j,bi,bj)
1235                dpopl  = dpopl *maskInC(i,j,bi,bj)
1236                dponl  = dponl *maskInC(i,j,bi,bj)
1237                dpofel = dpofel*maskInC(i,j,bi,bj)
1238                dpsil  = dpsil *maskInC(i,j,bi,bj)
1239                dnh4l  = dnh4l *maskInC(i,j,bi,bj)
1240                dno2l  = dno2l *maskInC(i,j,bi,bj)
1241                DO nz = 1,nzmax
1242                 dzoop (nz) = dzoop (nz)*maskInC(i,j,bi,bj)
1243                 dzoon (nz) = dzoon (nz)*maskInC(i,j,bi,bj)
1244                 dzoofe(nz) = dzoofe(nz)*maskInC(i,j,bi,bj)
1245                 dzoosi(nz) = dzoosi(nz)*maskInC(i,j,bi,bj)
1246                ENDDO
1247                DO np = 1,npmax
1248                 dPhy(np) = dPhy(np)*maskInC(i,j,bi,bj)
1249    #ifdef GEIDER
1250    #ifdef DYNAMIC_CHL
1251                 dphychl(np) = dphychl(np)*maskInC(i,j,bi,bj)
1252    #endif
1253    #endif
1254                ENDDO
1255    #ifdef ALLOW_CDOM
1256                dcdoml = dcdoml*maskInC(i,j,bi,bj)
1257    #endif
1258    #ifdef ALLOW_CARBON
1259                ddicl = ddicl*maskInC(i,j,bi,bj)
1260                ddocl = ddocl*maskInC(i,j,bi,bj)
1261                dpocl = dpocl*maskInC(i,j,bi,bj)
1262                dpicl = dpicl*maskInC(i,j,bi,bj)
1263                dalkl = dalkl*maskInC(i,j,bi,bj)
1264                do2l  = do2l *maskInC(i,j,bi,bj)
1265                DO nz = 1,nzmax
1266                 dzoocl(nz) = dzoocl(nz)*maskInC(i,j,bi,bj)
1267                ENDDO
1268    #endif
1269              ENDIF
1270    #endif
1271    
1272  c now update main tracer arrays  c now update main tracer arrays
1273            dtplankton = PTRACERS_dTLev(k)/float(nsubtime)            dtplankton = PTRACERS_dTLev(k)/float(nsubtime)
1274            Ptr(i,j,k,bi,bj,iPO4 ) = Ptr(i,j,k,bi,bj,iPO4) +            Ptr(i,j,k,bi,bj,iPO4 ) = Ptr(i,j,k,bi,bj,iPO4) +
# Line 1120  c in darwin_plankton this is stored for Line 1329  c in darwin_plankton this is stored for
1329  #endif  #endif
1330  #endif  #endif
1331           ENDDO           ENDDO
1332    #ifdef ALLOW_CDOM
1333              Ptr(i,j,k,bi,bj,iCDOM ) = Ptr(i,j,k,bi,bj,iCDOM ) +
1334         &                                  dtplankton*dcdoml
1335    #endif
1336  #ifdef ALLOW_CARBON  #ifdef ALLOW_CARBON
1337            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 ) +
1338       &                                  dtplankton*ddicl       &                                  dtplankton*ddicl
# Line 1250  c    &                       deltaTclock Line 1463  c    &                       deltaTclock
1463       &                           phychl(np)*dtplankton       &                           phychl(np)*dtplankton
1464               enddo               enddo
1465  #endif  #endif
1466    #ifdef DAR_DIAG_PARW
1467                do ilam=1,tlam
1468                   PARwave(i,j,k,bi,bj,ilam)=PARwave(i,j,k,bi,bj,ilam)+
1469         &                           PARw_k(ilam,k)*dtplankton
1470                enddo
1471                do np=1,npmax
1472                  chl2cave(i,j,k,bi,bj,np)=chl2cave(i,j,k,bi,bj,np)+
1473         &                          chl2cl(np)*dtplankton
1474                enddo
1475    #endif
1476  #ifdef DAR_DIAG_ACDOM  #ifdef DAR_DIAG_ACDOM
1477  c            print*,'acdom',k,acdom_k(k,darwin_diag_acdom_ilam)  c            print*,'acdom',k,acdom_k(k,darwin_diag_acdom_ilam)
1478               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 1271  Coj            no Eu at surface (yet) Line 1494  Coj            no Eu at surface (yet)
1494                 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)+
1495       &                                 Euz(ilam,k-1)*dtplankton       &                                 Euz(ilam,k-1)*dtplankton
1496                endif                endif
1497                  Estave(i,j,k,bi,bj,ilam)=Estave(i,j,k,bi,bj,ilam)+
1498         &                                 Estop(ilam,k)*dtplankton
1499                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)+
1500       &                                 Eutop(ilam,k)*dtplankton       &                                 Eutop(ilam,k)*dtplankton
1501               enddo               enddo
1502  #endif  #endif
1503    #ifdef DAR_DIAG_IRR_AMPS
1504                 do ilam = 1,tlam
1505                   amp1ave(i,j,k,bi,bj,ilam)=amp1ave(i,j,k,bi,bj,ilam)+
1506         &                                 amp1(ilam,k)*dtplankton
1507                   amp2ave(i,j,k,bi,bj,ilam)=amp2ave(i,j,k,bi,bj,ilam)+
1508         &                                 amp2(ilam,k)*dtplankton
1509                 enddo
1510    #endif
1511  #ifdef DAR_DIAG_ABSORP  #ifdef DAR_DIAG_ABSORP
1512               do ilam = 1,tlam               do ilam = 1,tlam
1513                 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 1532  Coj            no Eu at surface (yet)
1532       &                                 bbpart_k(k,ilam)*dtplankton       &                                 bbpart_k(k,ilam)*dtplankton
1533               enddo               enddo
1534  #endif  #endif
1535    #ifdef DAR_RADTRANS
1536                 if (k.eq.1) then
1537                   rmudave(i,j,bi,bj)=rmudave(i,j,bi,bj)+
1538         &                                 rmud*dtplankton
1539                 endif
1540    #endif
1541    #ifdef DAR_DIAG_EK
1542                do np=1,npmax
1543                 Ekave(i,j,k,bi,bj,np)=Ekave(i,j,k,bi,bj,np)+
1544         &                        Ekl(np)*dtplankton
1545                 EkoverEave(i,j,k,bi,bj,np)=EkoverEave(i,j,k,bi,bj,np)+
1546         &                        EkoverEl(np)*dtplankton
1547                 acclimave(i,j,k,bi,bj,np)=acclimave(i,j,k,bi,bj,np)+
1548         &                        accliml(np)*dtplankton
1549                 do ilam=1,tlam
1550                    Ek_nlave(i,j,k,bi,bj,np,ilam)=
1551         &                        Ek_nlave(i,j,k,bi,bj,np,ilam)+
1552         &                        Ek_nll(np,ilam)*dtplankton
1553                    EkoverE_nlave(i,j,k,bi,bj,np,ilam)=
1554         &                        EkoverE_nlave(i,j,k,bi,bj,np,ilam)+
1555         &                        EkoverE_nll(np,ilam)*dtplankton
1556                 enddo
1557                enddo
1558    #endif
1559  #ifdef DAR_DIAG_RSTAR  #ifdef DAR_DIAG_RSTAR
1560               do np=1,npmax               do np=1,npmax
1561                 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 1344  Coj            no Eu at surface (yet) Line 1601  Coj            no Eu at surface (yet)
1601       &                              FluxCO2(i,j,bi,bj)*dtplankton       &                              FluxCO2(i,j,bi,bj)*dtplankton
1602                SUROave(i,j,bi,bj)   =SUROave(i,j,bi,bj)+                SUROave(i,j,bi,bj)   =SUROave(i,j,bi,bj)+
1603       &                              flxO2(i,j)*dtplankton       &                              flxO2(i,j)*dtplankton
1604                 endif
1605    #ifdef pH_3D
1606                  pCO2ave(i,j,k,bi,bj)   =pCO2ave(i,j,k,bi,bj)+
1607         &                              pCO2(i,j,k,bi,bj)*dtplankton
1608                  pHave(i,j,k,bi,bj)     =pHave(i,j,k,bi,bj)+
1609         &                              pH(i,j,k,bi,bj)*dtplankton
1610    #else
1611                 if (k.eq.1) then
1612                pCO2ave(i,j,bi,bj)   =pCO2ave(i,j,bi,bj)+                pCO2ave(i,j,bi,bj)   =pCO2ave(i,j,bi,bj)+
1613       &                              pCO2(i,j,bi,bj)*dtplankton       &                              pCO2(i,j,bi,bj)*dtplankton
1614                pHave(i,j,bi,bj)     =pHave(i,j,bi,bj)+                pHave(i,j,bi,bj)     =pHave(i,j,bi,bj)+
1615       &                              pH(i,j,bi,bj)*dtplankton       &                              pH(i,j,bi,bj)*dtplankton
1616               endif               endif
1617  #endif  #endif
1618    #endif
1619            endif              endif  
1620  c end if hFac>0  c end if hFac>0
1621    
# Line 1404  C       reset the other slot for averagi Line 1670  C       reset the other slot for averagi
1670  C itistime  C itistime
1671  #endif  #endif
1672    
1673    #ifdef DAR_CHECK_IRR_CONT
1674           i = myXGlobalLo-1+(bi-1)*sNx+idiscEs
1675           j = myYGlobalLo-1+(bj-1)*sNy+jdiscEs
1676           write(6,'(I4.4,X,A,4(X,I4),1PE24.16)')myProcId,'max Es disc',
1677         &                                   i,j,kdiscEs,ldiscEs,discEs
1678           i = myXGlobalLo-1+(bi-1)*sNx+idiscEu
1679           j = myYGlobalLo-1+(bj-1)*sNy+jdiscEu
1680           write(6,'(I4.4,X,A,4(X,I4),1PE24.16)')myProcId,'max Eu disc',
1681         &                                   i,j,kdiscEu,ldiscEu,discEu
1682    #endif
1683    
1684  COJ fill diagnostics  COJ fill diagnostics
1685  #ifdef ALLOW_DIAGNOSTICS  #ifdef ALLOW_DIAGNOSTICS
1686         IF ( useDiagnostics ) THEN         IF ( useDiagnostics ) THEN
# Line 1446  c ANNA end TAVE Line 1723  c ANNA end TAVE
1723          WRITE(diagname,'(A8)') 'Diver4  '          WRITE(diagname,'(A8)') 'Diver4  '
1724          CALL DIAGNOSTICS_FILL( Diver4(1-Olx,1-Oly,1), diagname,          CALL DIAGNOSTICS_FILL( Diver4(1-Olx,1-Oly,1), diagname,
1725       &                         0,Nr,2,bi,bj,myThid )       &                         0,Nr,2,bi,bj,myThid )
1726            WRITE(diagname,'(A8)') 'Shannon '
1727            CALL DIAGNOSTICS_FILL( Shannon(1-Olx,1-Oly,1), diagname,
1728         &                         0,Nr,2,bi,bj,myThid )
1729            WRITE(diagname,'(A8)') 'Simpson '
1730            CALL DIAGNOSTICS_FILL( Simpson(1-Olx,1-Oly,1), diagname,
1731         &                         0,Nr,2,bi,bj,myThid )
1732  #endif  #endif
1733  #ifdef ALLOW_DIAZ  #ifdef ALLOW_DIAZ
1734  #ifdef DAR_DIAG_NFIXP  #ifdef DAR_DIAG_NFIXP
# Line 1477  c ANNA end TAVE Line 1760  c ANNA end TAVE
1760       &                         0,1,2,bi,bj,myThid )       &                         0,1,2,bi,bj,myThid )
1761          CALL DIAGNOSTICS_FILL( flxO2(1-Olx,1-Oly), 'DICOFLX ',          CALL DIAGNOSTICS_FILL( flxO2(1-Olx,1-Oly), 'DICOFLX ',
1762       &                         0,1,2,bi,bj,myThid )       &                         0,1,2,bi,bj,myThid )
1763    #ifdef pH_3D
1764            CALL DIAGNOSTICS_FILL( pCO2(1-Olx,1-Oly,1,bi,bj), 'DICPCO2 ',
1765         &                         0,Nr,2,bi,bj,myThid )
1766            CALL DIAGNOSTICS_FILL( pH(1-Olx,1-Oly,1,bi,bj), 'DICPHAV ',
1767         &                         0,Nr,2,bi,bj,myThid )
1768    #else
1769          CALL DIAGNOSTICS_FILL( pCO2(1-Olx,1-Oly,bi,bj), 'DICPCO2 ',          CALL DIAGNOSTICS_FILL( pCO2(1-Olx,1-Oly,bi,bj), 'DICPCO2 ',
1770       &                         0,1,2,bi,bj,myThid )       &                         0,1,2,bi,bj,myThid )
1771          CALL DIAGNOSTICS_FILL( pH(1-Olx,1-Oly,bi,bj), 'DICPHAV ',          CALL DIAGNOSTICS_FILL( pH(1-Olx,1-Oly,bi,bj), 'DICPHAV ',
1772       &                         0,1,2,bi,bj,myThid )       &                         0,1,2,bi,bj,myThid )
1773    #endif
1774  #endif /* ALLOW_CARBON */  #endif /* ALLOW_CARBON */
1775         ENDIF         ENDIF
1776  #endif /* ALLOW_DIAGNOSTICS */  #endif /* ALLOW_DIAGNOSTICS */
# Line 1493  c determine iron partitioning  - solve f Line 1783  c determine iron partitioning  - solve f
1783  c  c
1784  #ifdef ALLOW_TIMEAVE  #ifdef ALLOW_TIMEAVE
1785  c save averages  c save averages
1786         do k=1,nR           dar_timeave(bi,bj) = dar_timeave(bi,bj) + dtplankton
          dar_timeave(bi,bj,k)=dar_timeave(bi,bj,k)  
      &                         +dtplankton  
1787  #ifdef ALLOW_CARBON  #ifdef ALLOW_CARBON
1788           dic_timeave(bi,bj,k)=dic_timeave(bi,bj,k)           dic_timeave(bi,bj) = dic_timeave(bi,bj) + dtplankton
      &                         +dtplankton  
1789  #endif  #endif
        enddo  
1790  #endif  #endif
1791  c  c
1792  c -----------------------------------------------------  c -----------------------------------------------------
Legend:
Removed from v.1.4  
changed lines
  Added in v.1.18
  ViewVC Help
Powered by ViewVC 1.1.22