pkg/darwin/darwin_generate_random.F

8fbfd1f382 Oliv*0001 #include "DARWIN_OPTIONS.h"
                0002
                0003 CBOP
                0004 C !ROUTINE: DARWIN_GENERATE_RANDOM
                0005
                0006 C !INTERFACE: ==========================================================
                0007       SUBROUTINE DARWIN_GENERATE_RANDOM( myThid )
                0008
                0009 C !DESCRIPTION:
                0010 C     Generate parameters for plankton types using a "Monte Carlo"
                0011 C     approach.
                0012 C
                0013 C     Mick Follows, Scott Grant Fall/Winter 2005
                0014 C     Stephanie Dutkiewicz Spring/Summer 2005
                0015 C     Anna Hickman Summer 2008
                0016 C
                0017 C     DARWIN_TWO_SPECIES_SETUP
                0018 C      1=large, 2=small
                0019 C     DARWIN_NINE_SPECIES_SETUP
                0020 C      1=diatom, 2=other large, 3=syn, 4=hl pro, 5=ll pro, 6=trich,
                0021 C      7=uni diaz, 8=small euk, 9=cocco
                0022
                0023 C !USES: ===============================================================
                0024       IMPLICIT NONE
                0025 #include "EEPARAMS.h"
                0026 #ifdef ALLOW_RADTRANS
                0027 #include "RADTRANS_SIZE.h"
                0028 #endif
                0029 #ifdef ALLOW_DARWIN
                0030 #include "DARWIN_SIZE.h"
                0031 #include "DARWIN_INDICES.h"
                0032 #include "DARWIN_RADTRANS.h"
                0033 #include "DARWIN_PARAMS.h"
                0034 #include "DARWIN_TRAITPARAMS.h"
                0035 #include "DARWIN_TRAITS.h"
                0036 #endif
                0037
                0038 C !INPUT PARAMETERS: ===================================================
                0039 C  myThid               :: thread number
                0040       INTEGER myThid
                0041 CEOP
                0042
                0043 #ifdef ALLOW_DARWIN
                0044
                0045 C !FUNCTIONS: ==========================================================
                0046       _RL DARWIN_RANDOM
                0047       EXTERNAL DARWIN_RANDOM
                0048       _RL DARWIN_RANDOM_NORMAL
                0049       EXTERNAL DARWIN_RANDOM_NORMAL
                0050
                0051 C !LOCAL VARIABLES: ====================================================
                0052 C     msgBuf    - Informational/error meesage buffer
                0053       CHARACTER*(MAX_LEN_MBUF) msgBuf
                0054       integer iUnit1, iUnit2
                0055       INTEGER iMinPrey, iMaxPrey, iMinPred, iMaxPred
                0056       _RL RandNoSize
                0057       _RL RandNoDiatom
                0058       _RL RandNoDiazo
                0059       _RL RandNoGrow
                0060       _RL RandNoMort
                0061       _RL RandNoNsrc
                0062       _RL RandNoAPType3
                0063       _RL RandNoAPType2
                0064       _RL RandNoAPType1
                0065       _RL RandNoAPType4
                0066       _RL RandNoTemp
                0067       _RL RandNoKsat
                0068       _RL RandNoKsatPAR
                0069       _RL RandNoKinhPAR
                0070       _RL RandNoDummy
                0071       _RL RandNoGrowGeider
                0072       _RL RandNoYield
                0073       _RL RandNoChl2C
                0074
                0075       _RL growthdays
                0076       _RL mortdays
                0077       _RL pday
                0078       _RL year
                0079       _RL month
                0080       _RL fiveday
                0081       _RL rtime
                0082       _RL dm
                0083       _RL dmzoo(nplank)
                0084       _RL volp
                0085       _RL sm
                0086       _RL PI
                0087       INTEGER np,nz,l,i,iopt
                0088       INTEGER signvar
                0089       PARAMETER ( PI    = 3.14159265358979323844D0   )
                0090
                0091 C     used to be global in monod pkg
                0092       _RL phyto_esd(nplank)
                0093       _RL phyto_vol(nplank)
                0094
                0095       INTEGER physize(nplank)
                0096       INTEGER nsource(nplank)
                0097       INTEGER diacoc(nplank)
                0098
                0099 C ======================================================================
                0100
                0101 c length of day (seconds)
                0102       pday = 86400.0 _d 0
                0103
                0104       DO np = 1, nPlank
                0105        kexcc(np) = 0.0 _d 0
                0106        kexcn(np) = 0.0 _d 0
                0107        kexcp(np) = 0.0 _d 0
                0108        kexcsi(np) = 0.0 _d 0
                0109        kexcfe(np) = 0.0 _d 0
                0110
                0111        Qnmax(np) = 0.0 _d 0
                0112        Qnmin(np) = 0.0 _d 0
                0113        Qpmax(np) = 0.0 _d 0
                0114        Qpmin(np) = 0.0 _d 0
                0115        Qsimax(np) = 0.0 _d 0
                0116        Qsimin(np) = 0.0 _d 0
                0117        Qfemax(np) = 0.0 _d 0
                0118        Qfemin(np) = 0.0 _d 0
                0119
                0120        vmaxFeT(np) = 0.0 _d 0
                0121        vmaxNH4(np) = 0.0 _d 0
                0122        vmaxNO2(np) = 0.0 _d 0
                0123        vmaxNO3(np) = 0.0 _d 0
                0124        vmaxN(np) = 0.0 _d 0
                0125        vmaxPO4(np) = 0.0 _d 0
                0126        vmaxSiO2(np) = 0.0 _d 0
                0127
                0128 C      initialize discrete traits to "unset"
                0129        physize(np)    = -1
                0130        diacoc(np)     = -1
                0131        diazo(np) = -1
                0132        nsource(np)    = -1
                0133
                0134        isPhoto(np) = 0
                0135        isPrey(np) = 0
                0136        isPred(np) = 0
                0137       ENDDO
                0138
                0139 C ======================================================================
                0140 c     phytoplankton
                0141 C ======================================================================
                0142
                0143       DO np = 1, nPhoto
                0144
                0145       isPhoto(np) = 1
                0146       isPrey(np) = 1
                0147
                0148       Xmin(np) = phymin
                0149
                0150 C ======================================================================
                0151 C RANDOM NUMBERS
                0152
                0153 C pre-compute random numbers and discrete traits
                0154
                0155 C phyto either "small" (physize(np)=0.0) or "big" (physize(np)=1.0)
                0156 C at this point independent of whether diatom or coccolithophor or not
                0157       RandNoSize = darwin_random(myThid)
                0158       IF (physize(np).LT.0) THEN
                0159        IF(RandNoSize .GT. 0.500 _d 0)then
                0160         physize(np) = 1
                0161        ELSE
                0162         physize(np) = 0
                0163        ENDIF
                0164       ENDIF
                0165 #ifdef DARWIN_TWO_SPECIES_SETUP
                0166       IF (np.EQ.1) physize(np) = 1
                0167       IF (np.EQ.2) physize(np) = 0
                0168 #endif
                0169 #ifdef DARWIN_NINE_SPECIES_SETUP
                0170       IF (np.LT.3.or.np.EQ.6.or.np.EQ.9) then
                0171         physize(np) = 1
                0172       ELSE
                0173         physize(np) = 0
                0174       ENDIF
                0175 #endif
                0176 c
                0177 c phyto either diatoms (diacoc=1.0) and use silica or cocolithophor
                0178 c (diacoc=2.0) and produce PIC or neither (diacoc=0.0)
                0179 c if they are large
                0180       IF (physize(np).EQ.1) then
                0181         RandNoDiatom = darwin_random(myThid)
                0182         IF (diacoc(np) .LT. 0) THEN
                0183          IF(RandNoDiatom .GT. 0.500 _d 0)then
                0184           diacoc(np) = 1
                0185          ELSE
                0186           diacoc(np) = 0
                0187          ENDIF
                0188 c         IF(RandNo .GT. 0.670 _d 0)then
                0189 c           diacoc(np) = 1
                0190 c         ENDIF
                0191 c         IF(RandNo .GT. 0.330 _d 0 .AND. RandNo. le. 0.67 _d 0)then
                0192 c           diacoc(np) = 2
                0193 c         ENDIF
                0194 c         IF (RandNo .LE. 0.330 _d 0) then
                0195 c           diacoc(np) = 0
                0196 c         ENDIF
                0197         ENDIF
                0198       ELSE
                0199          diacoc(np) = 0
                0200       ENDIF
                0201 #ifdef DARWIN_TWO_SPECIES_SETUP
                0202       diacoc(np) = 0
                0203 #endif
                0204 #ifdef DARWIN_NINE_SPECIES_SETUP
                0205       IF (np.EQ.1) then
                0206         diacoc(np) = 1
                0207       ELSE
                0208         diacoc(np) = 0
                0209       ENDIF
                0210       IF (np.EQ.9) then
                0211         diacoc(np) = 2
                0212       ENDIF
                0213 #endif
                0214
                0215 c phyto either diazotrophs (diazo=1.0) or not (diazo=0.0)
                0216       RandNoDiazo = darwin_random(myThid)
                0217       IF (diazo(np) .LT. 0) THEN
                0218        IF(RandNoDiazo .GT. 0.6700 _d 0)then
                0219         diazo(np) = 1
                0220        ELSE
                0221         diazo(np) = 0
                0222        ENDIF
                0223       ENDIF
                0224 c TEST ...........................................
                0225 #ifndef DARWIN_ALLOW_DIAZ
                0226       diazo(np) = 0
                0227 #endif
                0228 c TEST ...........................................
                0229 #ifdef DARWIN_TWO_SPECIES_SETUP
                0230       diazo(np) = 0
                0231 #endif
                0232 #ifdef DARWIN_NINE_SPECIES_SETUP
                0233       IF (np.GT.5.AND.np.LT.8) then
                0234          diazo(np) = 1
                0235       ELSE
                0236          diazo(np) = 0
                0237       ENDIF
                0238 #endif
                0239
                0240       RandNoGrow = darwin_random(myThid)
                0241       RandNoMort = darwin_random(myThid)
                0242 c nutrient source
                0243       IF(diazo(np) .ne. 1)then
                0244          RandNoNsrc = darwin_random(myThid)
                0245          IF (nsource(np) .LT. 0) THEN
                0246           IF (physize(np).EQ.1) then
                0247            nsource(np) = 3
                0248           ELSE
                0249            IF(RandNoNsrc .GT. 0.670 _d 0)then
                0250              nsource(np) = 1
                0251            ELSEif(RandNoNsrc .LT. 0.33 _d 0)then
                0252              nsource(np) = 2
                0253            ELSE
                0254              nsource(np) = 3
                0255            ENDIF
                0256 c ANNA shift bias away from pros. Now equal chance of being HL, LL, Syn, Euk.
                0257 c ANNA i.e. now 50% chance of being Pro (nsource 1 or 2, with 50% change of each being HL)
                0258 c ANNA i.e. and 50% chance of being non-Pro (nsource 3, with 50% chance of non-pro being Syn)
                0259 c           IF(RandNo .GT. 0.50 _d 0)then
                0260 c             nsource(np) = 3
                0261 c           ELSEif(RandNo .LT. 0.25 _d 0)then
                0262 c             nsource(np) = 2
                0263 c           ELSE
                0264 c             nsource(np) = 1
                0265 c           ENDIF
                0266           ENDIF
                0267          ENDIF
                0268       ELSE
                0269          nsource(np) = 0
                0270       ENDIF
                0271 #ifdef DARWIN_TWO_SPECIES_SETUP
                0272       nsource(np) = 3
                0273 #endif
                0274 #ifdef DARWIN_NINE_SPECIES_SETUP
                0275       IF (np.LT.4) then
                0276         nsource(np) = 3
                0277       ENDIF
                0278       nsource(4)=2
                0279       nsource(5)=1
                0280       IF (np.GT.5.AND.np.LT.8) then
                0281         nsource(np) = 0
                0282       ENDIF
                0283       IF (np.GT.7) then
                0284         nsource(np) = 3
                0285       ENDIF
                0286 #endif
                0287
                0288 c.....................................................
                0289 c ANNA make selections for RADTRANS
                0290 c.....................................................
                0291 #ifdef ALLOW_RADTRANS
                0292 c for now, choice of four absorption spectra types
                0293 c pros get either 'HL' or 'LL'
                0294 c small others get 'syn' or 'euk'
                0295 c large get 'euk'
                0296 c each 'type', once assigned, gets given actual values in wavebands_init_vari.F
                0297
                0298 c ANNA_Q could use tricho abs and scattering spectra (Subramanian et al. 1999)
                0299 c ANNA_Q think diaz is turned off for now
                0300 c Diaz will be 0 if not defined, and will have nsource = 0.
                0301       print*,'nopt',nopt,np,nsource(np),aptype(np)
                0302       IF (nOpt.EQ.4) then
                0303        RandNoAPType3 = darwin_random(myThid)
                0304        RandNoAPType2 = darwin_random(myThid)
                0305        RandNoAPType1 = darwin_random(myThid)
                0306        IF (nsource(np).EQ.0) then   !if diazo
                0307         IF (physize(np).EQ.1.0d0) then !if BIG
                0308         aptype(np) = 1                !euk (assume diatom association)
                0309         ELSE                           !or
                0310         aptype(np) = 2                !syn (for now - tricho has billins)
                0311         ENDIF
                0312        ENDIF
                0313
                0314        IF (nsource(np).EQ.3) then !if all three sources (NO3)
                0315         IF (physize(np).EQ.1.0d0) then !if BIG
                0316         aptype(np) = 1                !euk
                0317         ELSE                           !if SMALL
                0318          IF (RandNoAPType3.GT.0.500d0) then
                0319          aptype(np) = 1               !euk
                0320          ELSE                          !or
                0321          aptype(np) = 2               !Syn
                0322          ENDIF
                0323         ENDIF
                0324        ENDIF
                0325
                0326        IF (nsource(np).EQ.2) then !if NH4 only
                0327         IF (RandNoAPType2.GT.0.500d0) then
                0328         aptype(np) = 3               !Pro HL
                0329         ELSE                          !or
                0330         aptype(np) = 4               !Pro LL
                0331         ENDIF
                0332        ENDIF
                0333
                0334        IF (nsource(np).EQ.1) then !if NH4 & NO2
                0335         IF (RandNoAPType1.GT.0.500d0) then
                0336         aptype(np) = 3               !Pro HL
                0337         ELSE                          !or
                0338         aptype(np) = 4               !Pro LL
                0339         ENDIF
                0340        ENDIF
                0341        print*,'ap',np,nsource(np),aptype(np)
                0342       ENDIF
                0343 c
                0344       IF (nOpt.EQ.12) then
                0345        RandNoAPType4 = darwin_random(myThid)
                0346        IF (nsource(np).EQ.0) then   !if diazo
                0347         IF (physize(np).EQ.1.0d0) then !if BIG
                0348          IF (diacoc(np).EQ.1.0d0) then
                0349            aptype(np) = 5                !diatom association
                0350          ENDIF
                0351          IF (diacoc(np).EQ.0.0d0) then
                0352            aptype(np) = 7                !tricho
                0353          ENDIF
                0354          IF (diacoc(np).EQ.2.0d0) then
                0355            aptype(np) = 6                !coccolithopher(?)
                0356          ENDIF
                0357         ELSE                           !or
                0358          aptype(np) = 1                !unicellular (whould be 8 -
                0359                                         !but currently zero)
                0360         ENDIF
                0361        ENDIF
                0362
                0363        IF (nsource(np).EQ.3) then !if all three sources (NO3)
                0364         IF (physize(np).EQ.1.0d0) then !if BIG
                0365          IF (diacoc(np).EQ.1.0d0) then
                0366            aptype(np) = 5                !diatom
                0367          ENDIF
                0368          IF (diacoc(np).EQ.0.0d0) then
                0369            aptype(np) = 9                !Lg Euk
                0370          ENDIF
                0371          IF (diacoc(np).EQ.2.0d0) then
                0372            aptype(np) = 6                !coccolithopher
                0373          ENDIF
                0374         ELSE                           !if SMALL
                0375          IF (RandNoAPType4.GT.0.500d0) then
                0376          aptype(np) = 1               !euk
                0377          ELSE                          !or
                0378          aptype(np) = 2               !Syn
                0379          ENDIF
                0380         ENDIF
                0381        ENDIF
                0382       ENDIF
                0383
                0384 #ifdef DARWIN_TWO_SPECIES_SETUP
                0385       IF (np.EQ.1) aptype(np) = 10
                0386       IF (np.EQ.2) aptype(np) = 10
                0387 #endif
                0388 #ifdef DARWIN_NINE_SPECIES_SETUP
                0389       IF (np.EQ.1) aptype(np) = 5
                0390       IF (np.EQ.2) aptype(np) = 9
                0391       IF (np.EQ.3) aptype(np) = 2
                0392       IF (np.EQ.4) aptype(np) = 3
                0393       IF (np.EQ.5) aptype(np) = 4
                0394       IF (np.EQ.6) aptype(np) = 7
                0395       IF (np.EQ.7) aptype(np) = 8
                0396       IF (np.EQ.8) aptype(np) = 1
                0397       IF (np.EQ.9) aptype(np) = 6
                0398       aptype(np) = 10
                0399 #endif
                0400
                0401       print*,'aptype',np,aptype(np)
                0402       iopt = aptype(np)
                0403       IF (1 .LE. iopt .AND. iopt .LE. nOpt) THEN
                0404         DO l = 1, nlam
                0405          aphy_chl(np,l) = aphy_chl_type(iopt,l)
                0406          aphy_chl_ps(np,l) = aphy_chl_ps_type(iopt,l)
f61b1017e2 Oliv*0407          aphy_mgC(np,l) = aphy_mgC_type(iopt,l)
8fbfd1f382 Oliv*0408          bphy_mgC(np,l) = bphy_mgC_type(iopt,l)
                0409          bbphy_mgC(np,l) = bbphy_mgC_type(iopt,l)
                0410         ENDDO
                0411       ELSE
                0412         WRITE(msgBuf,'(A,2I4)')'invalid optical phyto type:',np,iopt
                0413         CALL PRINT_ERROR( msgBuf, myThid )
                0414         STOP 'ABNORMAL END: S/R DARWIN_READTRAITS'
                0415       ENDIF
                0416
                0417 #else
                0418 c ANNA number of RandNo's carreid out MUST MATCH regardless of wavebands or not.
                0419 C ANNA the number of RandNo statements here MUST MATCH the number done above
                0420
                0421 c        RandNo = darwin_random(myThid)
                0422 c        RandNo = darwin_random(myThid)
                0423 c        RandNo = darwin_random(myThid)
                0424
                0425 #endif
                0426 c ANNA ENDIF
                0427
                0428
                0429       RandNoTemp = darwin_random(myThid)
                0430       RandNoKsat = darwin_random(myThid)
                0431 #ifndef DARWIN_ALLOW_GEIDER
                0432       IF(physize(np) .EQ. 1)then
                0433          RandNoKsatPAR = darwin_random_normal(myThid)
                0434          RandNoKinhPAR = darwin_random_normal(myThid)
                0435       ELSE
                0436 c QQ remove someday
                0437          RandNoDummy = darwin_random(myThid)
                0438          RandNoKsatPAR = darwin_random_normal(myThid)
                0439          RandNoKinhPAR = darwin_random_normal(myThid)
                0440       ENDIF
                0441 #endif
                0442 #ifdef DARWIN_ALLOW_GEIDER
                0443       RandNoGrowGeider = darwin_random(myThid)
                0444       RandNoYield = darwin_random(myThid)
                0445       RandNoChl2C = darwin_random(myThid)
                0446 #endif
                0447
                0448 C ======================================================================
                0449
                0450 c size of phytoplankton
                0451       IF(physize(np).EQ. 1)then
                0452         dm = 10. _d 0  ! diameter (micrometer)
                0453       ELSE
                0454         dm = 1. _d 0  ! diameter (micrometer)
                0455       ENDIF
                0456 c phytoplankton volume in micrometers cubed
                0457       volp=4. _d 0/3. _d 0 *PI*(dm/2. _d 0)**3 _d 0
                0458 c
                0459 c common block variables (in m and m3)
                0460       phyto_esd(np)=dm* 1. _d -6
                0461       phyto_vol(np)=volp* 1. _d -18
                0462
                0463 c growth rates
                0464 c big/small phyto growth rates..
                0465       IF(physize(np) .EQ. 1)then
                0466         growthdays = Biggrow +RandNoGrow*Biggrowrange
                0467       ELSE
                0468         growthdays = Smallgrow +RandNoGrow*Smallgrowrange
                0469       ENDIF
                0470 #ifdef DARWIN_TWO_SPECIES_SETUP
                0471       IF(physize(np) .EQ. 1)then
                0472         growthdays = Biggrow
                0473       ELSE
                0474         growthdays = Smallgrow
                0475       ENDIF
                0476 #endif
                0477 #ifdef DARWIN_NINE_SPECIES_SETUP
                0478       IF(physize(np) .EQ. 1)then
                0479         growthdays = Biggrow
                0480       ELSE
                0481         growthdays = Smallgrow
                0482       ENDIF
                0483 #endif
                0484 c but diazotrophs always slower due to energetics
                0485       IF(diazo(np) .EQ. 1) then
                0486           growthdays = growthdays * diaz_growfac
                0487       ENDIF
                0488 c cocco have slower growth than other large
                0489       IF (diacoc(np).EQ.2. _d 0) then
                0490          growthdays= growthdays * cocco_growfac
                0491       ENDIF
                0492 c diatom has faster thatn other large
                0493       IF (diacoc(np).EQ.1. _d 0) then
                0494          growthdays= growthdays * diatom_growfac
                0495       ENDIF
                0496 c now convert to a growth rate
                0497       IF (growthdays.GT.0. _d 0) then
                0498        PCmax(np) = 1.0 _d 0/(growthdays*pday)
                0499       ELSE
                0500        PCmax(np) = 0. _d 0
                0501       ENDIF
                0502
                0503 c mortality and export fraction rates
                0504 c big/small phyto mortality rates..
                0505       IF(physize(np) .EQ. 1)then
                0506         mortdays = Bigmort +RandnoMort*Bigmortrange
                0507         ExportFracMort(np)=Bigexport
                0508         ExportFracMort2(np)=Bigexport
                0509 #ifdef DARWIN_ALLOW_EXUDE
                0510         ExportFracExude(np)=Bigexport
                0511 #else
                0512         ExportFracExude(np)=DARWIN_UNINIT_RL
                0513 #endif
                0514       ELSE
                0515         mortdays = Smallmort +RandNoMort*Smallmortrange
                0516         ExportFracMort(np)=Smallexport
                0517         ExportFracMort2(np)=Smallexport
                0518 #ifdef DARWIN_ALLOW_EXUDE
                0519         ExportFracExude(np)=Smallexport
                0520 #else
                0521         ExportFracExude(np)=DARWIN_UNINIT_RL
                0522 #endif
                0523       ENDIF
59e7a82054 Oliv*0524       FracExudeC(np) = 0.3 _d 0
8fbfd1f382 Oliv*0525 #ifdef DARWIN_TWO_SPECIES_SETUP
                0526       IF(physize(np) .EQ. 1)then
                0527         mortdays = Bigmort
                0528       ELSE
                0529         mortdays = Smallmort
                0530       ENDIF
                0531 #endif
                0532 #ifdef DARWIN_NINE_SPECIES_SETUP
                0533       IF(physize(np) .EQ. 1)then
                0534         mortdays = Bigmort
                0535       ELSE
                0536         mortdays = Smallmort
                0537       ENDIF
                0538 #endif
                0539
                0540 c now convert to a mortality rate
                0541       IF (mortdays.GT.0. _d 0) then
                0542         mort(np) = 1.0 _d 0/(mortdays*pday)
                0543       ELSE
                0544         mort(np) = 0. _d 0
                0545       ENDIF
                0546       mort2(np) = 0.0 _d 0
                0547
                0548 C phytoplankton do not have temperature-dependent mortality
                0549       tempMort(np) = 0
                0550       tempMort2(np) = 0
                0551
                0552       respRate(np) = 0.0 _d 0
                0553
                0554
                0555 c..........................................................
                0556 c generate phyto Temperature Function parameters
                0557 c.......................................................
                0558       phytoTempCoeff(np) = tempcoeff1
                0559       phytoTempExp1(np) = tempcoeff3
59e7a82054 Oliv*0560       phytoTempAe(np) = 0.0438 _d 0
8fbfd1f382 Oliv*0561       IF(physize(np) .EQ. 1)then
                0562         phytoTempExp2(np) = tempcoeff2_big
                0563       ELSE
                0564         phytoTempExp2(np) = tempcoeff2_small
                0565       ENDIF
                0566
                0567 #ifdef DARWIN_TEMP_RANGE
                0568 cswd    phytoTempOptimum(np) = 30.0 _d 0 - RandNo*28.0 _d 0
                0569       phytoTempOptimum(np) = tempmax - RandNoTemp*temprange
                0570       phytoDecayPower(np) = tempdecay
                0571 #else
                0572       phytoTempOptimum(np) = 0. _d 0
                0573       phytoDecayPower(np) = 0. _d 0
                0574 #endif
                0575
                0576 c stoichiometric ratios for each functional group of phyto
                0577 c relative to phosphorus  - the base currency nutrient
                0578 c set Si:P
                0579       IF(diacoc(np) .EQ. 1)then
                0580         R_SiC(np) =  val_R_SiC_diatom
                0581       ELSE
                0582         R_SiC(np) = 0.0 _d 0
                0583       ENDIF
                0584       IF(diacoc(np) .EQ. 2)then
                0585         R_PICPOC(np) =  val_R_PICPOC
                0586       ELSE
                0587         R_PICPOC(np) = 0.0 _d 0
                0588       ENDIF
                0589 c set N:P and iron requirement according to diazotroph status
                0590       IF(diazo(np) .EQ. 1)then
                0591         R_NC(np) = val_R_NC_diaz
                0592         R_FeC(np) =  val_R_FeC_diaz
                0593       ELSE
                0594         R_NC(np) = val_R_NC
                0595         R_FeC(np) = val_R_FeC
                0596       ENDIF
                0597 c set C:P ratio
                0598         R_PC(np) = val_R_PC
                0599 c set sinking rates according to allometry
                0600       IF(physize(np) .EQ. 1)then
                0601          biosink(np) = BigSink
                0602       ELSE
                0603          biosink(np) = SmallSink
                0604       ENDIF
                0605       bioswim(np) = 0.0 _d 0
                0606 c half-saturation coeffs
                0607
                0608       IF(physize(np) .EQ. 1)then
                0609          ksatPO4(np) = BigPsat + RandNoKsat*BigPsatrange
                0610       ELSE
                0611 c          ksatPO4(np) = SmallPsat + RandNoKsat*SmallPsatrange
                0612 c          if (nsource(np).LT.3) then
                0613 c            ksatPO4(np) = ksatPO4(np)*prochlPsat
                0614 c           ENDIF
                0615          IF (nsource(np).EQ.3) then
                0616            ksatPO4(np) = SmallPsat + RandNoKsat*SmallPsatrange
                0617          ENDIF
                0618          IF (nsource(np).EQ..0) then
                0619 c            ksatPO4(np) = SmallPsat + RandNoKsat*SmallPsatrange
                0620            ksatPO4(np) = UniDzPsat + RandNoKsat*UniDzPsatrange
                0621          ENDIF
                0622          IF (nsource(np).EQ.2.or.nsource(np).EQ.1) then
                0623            ksatPO4(np) = ProcPsat + RandNoKsat*ProcPsatrange
                0624          ENDIF
                0625       ENDIF
                0626       IF (diacoc(np) .EQ. 2) THEN
                0627          ksatPO4(np) = CoccoPsat + RandNoKsat*CoccoPsatrange
                0628       ENDIF
                0629 #ifdef DARWIN_TWO_SPECIES_SETUP
                0630       IF(physize(np) .EQ. 1)then
                0631          ksatPO4(np) = BigPsat
                0632       ELSE
                0633          ksatPO4(np) = SmallPsat
                0634       ENDIF
                0635 #endif
                0636 #ifdef DARWIN_NINE_SPECIES_SETUP
                0637       IF(physize(np) .EQ. 1)then
                0638          ksatPO4(np) = BigPsat
                0639       ELSE
                0640          ksatPO4(np) = SmallPsat
                0641       ENDIF
                0642       IF (nsource(np).EQ.2.or.nsource(np).EQ.1) then
                0643          ksatPO4(np) = ProcPsat
                0644       ENDIF
                0645       IF (diacoc(np) .EQ. 2) then
                0646          ksatPO4(np) = ksatPO4(np)/0.8 _d 0
                0647       ENDIF
                0648 #endif
                0649
                0650       ksatNO3(np) = ksatPO4(np)*R_NC(np)/R_PC(np)
                0651       ksatNO2(np) = ksatNO3(np)*ksatNO2fac
                0652 c Made ksatNH4 smaller since it is the preferred source
                0653       ksatNH4(np) = ksatNO3(np)*ksatNH4fac
                0654       ksatFeT(np) = ksatPO4(np)*R_FeC(np)/R_PC(np)
                0655       ksatSiO2(np) = val_ksatsio2
                0656       amminhib(np) = val_amminhib
                0657
                0658       acclimtimescl(np) = val_acclimtimescl
                0659
                0660 #ifndef DARWIN_ALLOW_GEIDER
                0661       R_ChlC(np) = val_R_ChlC
                0662
                0663 cNEW Light parameters:
                0664 c     ksatPAR {0.1 - 1.3}
                0665 c     0.35=Av High Light Adapted, 0.8=Av Low Light Adapted
                0666 c     kinhPAR  {0.0 - 3.0}
                0667 c     0.5 =Av High Light Adapted, 2.0=Av Low Light Adapted
                0668 c High Light Groups for Large size:
                0669       IF(physize(np) .EQ. 1)then
                0670          ksatPAR(np) = abs(Bigksatpar+Bigksatparstd*RandNoKsatPAR)
                0671
                0672          kinhPAR(np) = abs(BigkinhPAR+BigkinhPARstd*RandNoKinhPAR)
                0673       ELSE
                0674 c QQ remove someday
                0675 c Low Light Groups for Small size:
                0676          ksatPAR(np) = abs(smallksatpar+smallksatparstd*RandNoKsatPAR)
                0677
                0678          kinhPAR(np) = abs(smallkinhPAR+smallkinhPARstd*RandNoKinhPAR)
                0679       ENDIF
                0680 #ifdef DARWIN_TWO_SPECIES_SETUP
                0681       IF(physize(np) .EQ. 1)then
                0682          ksatPAR(np) = abs(Bigksatpar)
                0683          kinhPAR(np) = abs(Bigkinhpar)
                0684       ELSE
                0685          ksatPAR(np) = abs(smallksatpar)
                0686          kinhPAR(np) = abs(smallkinhpar)
                0687       ENDIF
                0688 #endif
                0689 #ifdef DARWIN_NINE_SPECIES_SETUP
                0690       IF(physize(np) .EQ. 1)then
                0691          ksatPAR(np) = abs(Bigksatpar)
                0692          kinhPAR(np) = abs(BigkinhPAR)
                0693       ELSE
                0694          ksatPAR(np) = abs(smallksatpar)
                0695          kinhPAR(np) = abs(smallkinhPAR)
                0696       ENDIF
                0697       IF (np.EQ.5) then
                0698         kinhPAR(np) = abs(LLProkinhPAR)
                0699       ENDIF
                0700       IF (np.EQ.9) then
                0701         kinhPAR(np) = abs(CoccokinhPAR)
                0702       ENDIF
                0703 #endif
                0704 #endif
                0705
                0706 #ifdef DARWIN_ALLOW_GEIDER
                0707 c big/small phyto growth rates..
                0708       IF(physize(np) .EQ. 1)then
                0709         growthdays = Biggrow +RandNoGrowGeider*Biggrowrange
                0710       ELSE
                0711         growthdays = Smallgrow +RandNoGrowGeider*Smallgrowrange
                0712       ENDIF
                0713 #ifdef DARWIN_TWO_SPECIES_SETUP
                0714       IF(physize(np) .EQ. 1)then
                0715         growthdays = Biggrow
                0716       ELSE
                0717         growthdays = Smallgrow
                0718       ENDIF
                0719 #endif
                0720 #ifdef DARWIN_NINE_SPECIES_SETUP
                0721       IF(physize(np) .EQ. 1)then
                0722         growthdays = Biggrow
                0723       ELSE
                0724         growthdays = Smallgrow
                0725       ENDIF
                0726 #endif
                0727 c but diazotrophs always slower due to energetics
                0728       IF(diazo(np) .EQ. 1) then
                0729           growthdays = growthdays * diaz_growfac
                0730       ENDIF
                0731 c cocco have slower growth than other large
                0732       IF (diacoc(np).EQ.2. _d 0) then
                0733          growthdays= growthdays * cocco_growfac
                0734       ENDIF
                0735 c diatom has faster thatn other large
                0736       IF (diacoc(np).EQ.1. _d 0) then
                0737          growthdays= growthdays * diatom_growfac
                0738       ENDIF
                0739 c now convert to a growth rate
                0740       IF (growthdays.GT.0. _d 0) then
                0741        pcmax(np) = 1.0 _d 0/(growthdays*pday)
                0742       ELSE
                0743        pcmax(np) = 0. _d 0
                0744       ENDIF
                0745 c
                0746 c photo-inhibition
                0747 #ifdef ALLOW_RADTRANS
                0748 c only LL Pro are inhibited
                0749        IF (aptype(np).EQ.4) then
                0750           inhibGeider(np) = inhibcoef_geid_val
                0751        ELSE
                0752           inhibGeider(np) = 0. _d 0
                0753        ENDIF
                0754 #else
                0755 c no inhibition
                0756        IF(physize(np) .EQ. 1)then
                0757          inhibGeider(np) = 0. _d 0
                0758        ELSE
                0759          inhibGeider(np) = 0. _d 0  !inhibcoef_geid_val
                0760        ENDIF
                0761 #endif
                0762 c
                0763
                0764 c big/small phyto PI slope (chl specific)
                0765 c        IF(physize(np) .EQ. 1)then
                0766 c          alphachl(np) = Bigalphachl +Randno*Bigalphachlrange
                0767 c        ELSE
                0768 c          alphachl(np) = Smallalphachl +RandNo*Smallalphachlrange
                0769 c        ENDIF
                0770
                0771 c ANNA gieder via mQyield instead of alpha
                0772 c big/small phyto Maximum Quantum Yield
                0773       IF(physize(np) .EQ. 1)then
                0774         mQyield(np) = BigmQyield +RandNoYield*BigmQyieldrange
                0775       ELSE
                0776         mQyield(np) = SmallmQyield +RandNoYield*SmallmQyieldrange
                0777       ENDIF
                0778 #ifdef DARWIN_TWO_SPECIES_SETUP
                0779       IF(physize(np) .EQ. 1)then
                0780         mQyield(np) = BigmQyield
                0781       ELSE
                0782         mQyield(np) = SmallmQyield
                0783       ENDIF
                0784 #endif
                0785 #ifdef DARWIN_NINE_SPECIES_SETUP
                0786       IF(physize(np) .EQ. 1)then
                0787         mQyield(np) = BigmQyield
                0788       ELSE
                0789         mQyield(np) = SmallmQyield
                0790       ENDIF
                0791 #endif
                0792 #ifdef ALLOW_RADTRANS
                0793 c ANNA for wavebands only, re-set mQyield to be constant for all np's
                0794 c ANNA i.e. let alpha vary only with aphy_chl_ps
                0795 c ANNA value is mean of vals for big and small.
                0796         mQyield(np) = 4.0 _d -5
                0797 #endif
                0798
                0799
                0800 c big/small phyto C:Chl max
                0801       IF(physize(np) .EQ. 1)then
                0802         chl2cmax(np) = Bigchl2cmax +RandnoChl2C*Bigchl2cmaxrange
                0803       ELSE
                0804         chl2cmax(np) = Smallchl2cmax +RandNoChl2C*Smallchl2cmaxrange
                0805       ENDIF
                0806 #ifdef DARWIN_TWO_SPECIES_SETUP
                0807       IF(physize(np) .EQ. 1)then
                0808         chl2cmax(np) = Bigchl2cmax
                0809       ELSE
                0810         chl2cmax(np) = Smallchl2cmax
                0811       ENDIF
                0812 #endif
                0813 #ifdef DARWIN_NINE_SPECIES_SETUP
                0814       IF(physize(np) .EQ. 1)then
                0815         chl2cmax(np) = Bigchl2cmax
                0816       ELSE
                0817         chl2cmax(np) = Smallchl2cmax
                0818       ENDIF
                0819 #endif
                0820 c ANNA chl2cmin added
                0821 c       chl2cmin(np) = 0.003  _d 0 * 12. _d 0  ! mg Chl a/mmol C
                0822
                0823 #endif /* DARWIN_ALLOW_GEIDER */
                0824
                0825       print*,'nsource',np,nsource(np)
                0826       IF (nsource(np) .EQ. 3) THEN
                0827         useNH4(np) = 1
                0828         useNO2(np) = 1
                0829         useNO3(np) = 1
                0830         combNO(np) = 1
                0831       ELSEIF (nsource(np) .EQ. 2) THEN
                0832         useNH4(np) = 1
                0833         useNO2(np) = 0
                0834         useNO3(np) = 0
                0835         combNO(np) = 0
                0836       ELSEIF (nsource(np) .EQ. 1) THEN
                0837         useNH4(np) = 1
                0838         useNO2(np) = 1
                0839         useNO3(np) = 0
                0840         combNO(np) = 0
                0841       ELSE
                0842         useNH4(np) = 0
                0843         useNO2(np) = 0
                0844         useNO3(np) = 0
                0845         combNO(np) = 0
                0846       ENDIF
                0847
                0848       IF (diacoc(np) .NE. 1) THEN
                0849         ksatSiO2(np) = 0.0
                0850         vmaxSiO2(np) = 0.0
                0851         R_SiC(np) = 0.0
                0852       ENDIF
                0853
                0854       IF (diacoc(np) .EQ. 1) THEN
                0855         hasSi(np) = 1
                0856         hasPIC(np) = 0
                0857       ELSEIF (diacoc(np) .EQ. 2) THEN
                0858         hasSi(np) = 0
                0859         hasPIC(np) = 1
                0860       ELSE
                0861         hasSi(np) = 0
                0862         hasPIC(np) = 0
                0863       ENDIF
                0864
                0865       ENDDO  ! np
                0866
                0867       DO np = 1, nplank
                0868
                0869         bactType(np) = 0
                0870         isAerobic(np) = 0
                0871         isDenit(np) = 0
                0872
                0873         yield(np) = 1.0 _d 0
                0874         yieldO2(np) = 1.0 _d 0
                0875         yieldNO3(np) = 1.0 _d 0
                0876
                0877         ksatPON(np) = 1.0 _d 0
                0878         ksatPOC(np) = 1.0 _d 0
                0879         ksatPOP(np) = 1.0 _d 0
                0880         ksatPOFe(np) = 1.0 _d 0
                0881         ksatDON(np) = 1.0 _d 0
                0882         ksatDOC(np) = 1.0 _d 0
                0883         ksatDOP(np) = 1.0 _d 0
                0884         ksatDOFe(np) = 1.0 _d 0
                0885
59e7a82054 Oliv*0886         hetTempAe(np) = 0.0438 _d 0
                0887         hetTempExp2(np) = 0.001 _d 0
                0888         hetTempOptimum(np) = 2 _d 0
                0889         hetDecayPower(np) = 4 _d 0
                0890
8fbfd1f382 Oliv*0891       ENDDO  ! np
                0892
                0893 #ifdef ALLOW_RADTRANS
                0894       DO np = nPhoto + 1, nplank
                0895         DO l = 1, nlam
                0896          bphy_mgC(np,l) = 0 _d 0
                0897          bbphy_mgC(np,l) = 0 _d 0
                0898         ENDDO
                0899       ENDDO
                0900 #endif
                0901
                0902
                0903 C ======================================================================
                0904 c     zooplankton
                0905 C ======================================================================
                0906
                0907       DO nz = 1, nplank
                0908        DO np = 1, nplank
                0909         palat(np,nz) = 0 _d 0
                0910         asseff(np,nz) = 0 _d 0
                0911         ExportFracPreyPred(np,nz) = 0 _d 0
                0912        ENDDO
                0913       ENDDO
                0914
                0915       iMinPrey = 1
                0916       iMaxPrey = nPhoto
                0917       iMinPred = nPhoto + 1
                0918       iMaxPred = nplank
                0919       DO nz = iMinPred, iMaxPred
                0920        isPred(nz) = 1
                0921       ENDDO
                0922
                0923       IF ( oldTwoGrazers ) THEN
                0924 c assume zoo(1) = small, zoo(2) = big
                0925
                0926        IF ( iMaxPred-iMinPred .NE. 1 ) THEN
                0927         WRITE(msgBuf,'(2A)') 'DARWIN_GENERATE: ',
                0928      &    'must have exactly 2 predators when oldTwoGrazers=.TRUE.'
                0929         CALL PRINT_ERROR( msgBuf , 1)
                0930         STOP 'ABNORMAL END: S/R DARWIN_GENERATE'
                0931        ENDIF
                0932        physize(iMinPred) = 0
                0933        physize(iMaxPred) = 1
                0934        grazemax(iMinPred) = GrazeFast
                0935        grazemax(iMaxPred) = GrazeFast
                0936        ExportFracMort(iMinPred) = ZooexfacSmall
                0937        ExportFracMort(iMaxPred) = ZooexfacBig
                0938        ExportFracMort2(iMinPred) = ZooexfacSmall
                0939        ExportFracMort2(iMaxPred) = ZooexfacBig
                0940 #ifdef DARWIN_ALLOW_EXUDE
                0941        ExportFracExude(iMinPred) = ZooexfacSmall
                0942        ExportFracExude(iMaxPred) = ZooexfacBig
                0943 #else
                0944        ExportFracExude(iMinPred) = DARWIN_UNINIT_RL
                0945        ExportFracExude(iMaxPred) = DARWIN_UNINIT_RL
                0946 #endif
59e7a82054 Oliv*0947        FracExudeC(iMinPred) = 0.3 _d 0
                0948        FracExudeC(iMaxPred) = 0.3 _d 0
8fbfd1f382 Oliv*0949        mort(iMinPred) = ZoomortSmall
                0950        mort(iMaxPred) = ZoomortBig
                0951        mort2(iMinPred) = ZoomortSmall2
                0952        mort2(iMaxPred) = ZoomortBig2
                0953        DO np = iMinPrey, iMaxPrey
                0954         ExportFracPreyPred(np,iMinPred) = ExGrazFracSmall
                0955         ExportFracPreyPred(np,iMaxPred) = ExGrazFracBig
                0956        ENDDO
                0957        dmzoo(iMinPred) = 30. _d 0   ! diameter (micrometer)
                0958        dmzoo(iMaxPred) = 300. _d 0  ! diameter (micrometer)
                0959 c palatibity according to "allometry"
                0960 c big grazers preferentially eat big phyto etc...
                0961        DO nz = iMinPred, iMaxPred
                0962         DO np = iMinPrey, iMaxPrey
                0963           IF (physize(nz).EQ.physize(np)) then
                0964             palat(np,nz) = palathi
                0965             asseff(np,nz) = GrazeEffmod
                0966           ELSE
                0967             palat(np,nz) = palatlo
                0968             IF (physize(np).EQ.0. _d 0) then
                0969               asseff(np,nz) = GrazeEffhi
                0970             ELSE
                0971               asseff(np,nz) = GrazeEfflow
                0972             ENDIF
                0973           ENDIF
                0974 c diatoms even less palatible
                0975           IF (diacoc(np).EQ.1. _d 0) then
                0976             palat(np,nz)= palat(np,nz)*diatomgraz
                0977           ENDIF
                0978 c coccolithophes less palatible
                0979           IF (diacoc(np).EQ.2. _d 0) then
                0980             palat(np,nz)= palat(np,nz)*coccograz
                0981           ENDIF
                0982 c other large phyto less palatible
                0983           IF (diacoc(np).EQ.0. _d 0 .AND.physize(np).EQ.1. _d 0) then
                0984             palat(np,nz)= palat(np,nz)*olargegraz
                0985           ENDIF
                0986 c need something in here for tricho
                0987         ENDDO
                0988        ENDDO
                0989
                0990 c     not oldTwoGrazers
                0991       ELSE
                0992
                0993        DO nz = iMinPred, iMaxPred
                0994         grazemax(nz) = GrazeRate
                0995         ExportFracMort(nz) = Zooexfac
                0996         ExportFracMort2(nz) = Zooexfac
                0997 #ifdef DARWIN_ALLOW_EXUDE
                0998         ExportFracExude(nz) = Zooexfac
                0999 #else
                1000         ExportFracExude(nz) = DARWIN_UNINIT_RL
                1001 #endif
59e7a82054 Oliv*1002         FracExudeC(nz) = 0.3 _d 0
8fbfd1f382 Oliv*1003         mort(nz) = Zoomort
                1004         mort2(nz) = Zoomort2
                1005         dmzoo(nz) = ZooDM
                1006         DO np = iMinPrey, iMaxPrey
                1007          palat(np,nz) = val_palat
                1008          asseff(np,nz) = val_ass_eff
                1009          ExportFracPreyPred(np,nz) = ExGrazFrac
                1010         ENDDO
                1011        ENDDO
                1012
                1013 c     oldTwoGrazers
                1014       ENDIF
                1015
                1016 c
                1017       DO nz = iMinPred, iMaxPred
                1018         R_NC(nz) = val_R_NC_zoo
                1019         R_PC(nz) = val_R_PC_zoo
                1020         R_SiC(nz) = val_R_SiC_zoo
                1021         R_FeC(nz) = val_R_FeC_zoo
                1022         R_ChlC(nz) = val_R_ChlC_zoo
                1023         R_PICPOC(nz) = val_R_PICPOC_zoo
                1024
                1025         Xmin(nz) = 0.0 _d 0
                1026
                1027         kgrazesat(nz) = kgrazesat_val
                1028 C zooplankton do have temperature-dependent mortality
4cb675e367 Oliv*1029         tempMort(nz) = 0
                1030         tempMort2(nz) = 0
                1031         tempGraz(nz) = 0
59e7a82054 Oliv*1032         grazTempAe(nz) = 0.0438 _d 0
                1033         grazTempExp2(nz) = 0.001 _d 0
                1034         grazTempOptimum(nz) = 2 _d 0
                1035         grazDecayPower(nz) = 4 _d 0
8fbfd1f382 Oliv*1036
                1037         respRate(nz) = 0.0 _d 0
                1038         biosink(nz) = 0.0 _d 0
                1039         bioswim(nz) = 0.0 _d 0
                1040
                1041 c zooplankton volume in micrometers cubed
                1042         volp = 4. _d 0/3. _d 0 *PI*(dmzoo(nz)/2. _d 0)**3 _d 0
                1043 c
                1044 c common block variables (in m and m3)
                1045         phyto_esd(nz) = dmzoo(nz)* 1. _d -6
                1046         phyto_vol(nz) = volp* 1. _d -18
                1047       ENDDO
                1048
                1049 #endif /* ALLOW_DARWIN */
                1050
                1051       RETURN
                1052       END
                1053