Changeset 6841 for branches/CNRS
- Timestamp:
- 2016-08-03T16:59:29+02:00 (8 years ago)
- Location:
- branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES
- Files:
-
- 22 edited
-
P4Z/p4zagg.F90 (modified) (3 diffs)
-
P4Z/p4zche.F90 (modified) (6 diffs)
-
P4Z/p4zfechem.F90 (modified) (4 diffs)
-
P4Z/p4zligand.F90 (modified) (2 diffs)
-
P4Z/p4zlim.F90 (modified) (1 diff)
-
P4Z/p4zlys.F90 (modified) (1 diff)
-
P4Z/p4zmeso.F90 (modified) (1 diff)
-
P4Z/p4zmort.F90 (modified) (1 diff)
-
P4Z/p4zopt.F90 (modified) (9 diffs)
-
P4Z/p4zpoc.F90 (modified) (16 diffs)
-
P4Z/p4zprod.F90 (modified) (10 diffs)
-
P4Z/p4zrem.F90 (modified) (2 diffs)
-
P4Z/p4zsed.F90 (modified) (1 diff)
-
P4Z/p4zsink.F90 (modified) (5 diffs)
-
P4Z/p4zsms.F90 (modified) (1 diff)
-
P5Z/p5zagg.F90 (modified) (3 diffs)
-
P5Z/p5zpoc.F90 (modified) (16 diffs)
-
P5Z/p5zprod.F90 (modified) (1 diff)
-
P5Z/p5zrem.F90 (modified) (2 diffs)
-
P5Z/p5zsed.F90 (modified) (1 diff)
-
P5Z/p5zsink.F90 (modified) (3 diffs)
-
P5Z/p5zsms.F90 (modified) (1 diff)
Legend:
- Unmodified
- Added
- Removed
-
branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zagg.F90
r6453 r6841 58 58 IF( nn_timing == 1 ) CALL timing_start('p4z_agg') 59 59 60 ! Initialization of some global variables61 ! ---------------------------------------62 prodpoc(:,:,:) = 0.63 conspoc(:,:,:) = 0.64 prodgoc(:,:,:) = 0.65 consgoc(:,:,:) = 0.66 60 ! 67 61 ! Exchange between organic matter compartments due to coagulation/disaggregation … … 100 94 zaggdoc2 = ( 3.53E3 * zfact + 0.1 * zstep ) * trb(ji,jj,jk,jpgoc) * 0.3 * trb(ji,jj,jk,jpdoc) 101 95 ! tranfer of DOC to POC due to brownian motion 102 zaggdoc3 = ( 5095. * 0. * trb(ji,jj,jk,jppoc) +114. * 0.3 * trb(ji,jj,jk,jpdoc) ) *zstep * 0.3 * trb(ji,jj,jk,jpdoc)96 zaggdoc3 = ( 114. * 0.3 * trb(ji,jj,jk,jpdoc) ) *zstep * 0.3 * trb(ji,jj,jk,jpdoc) 103 97 104 98 ! Update the trends … … 109 103 tra(ji,jj,jk,jpdoc) = tra(ji,jj,jk,jpdoc) - zaggdoc - zaggdoc2 - zaggdoc3 110 104 ! 111 prodpoc(ji,jj,jk) = prodpoc(ji,jj,jk) + zaggdoc + zaggdoc3 112 conspoc(ji,jj,jk) = conspoc(ji,jj,jk) - zagg 105 conspoc(ji,jj,jk) = conspoc(ji,jj,jk) - zagg + zaggdoc + zaggdoc3 113 106 prodgoc(ji,jj,jk) = prodgoc(ji,jj,jk) + zagg + zaggdoc2 114 107 ! -
branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zche.F90
r6455 r6841 150 150 !!--------------------------------------------------------------------- 151 151 INTEGER :: ji, jj, jk 152 REAL(wp) :: ztkel, ztkel1, zt , z t2 , zsal , zsal2 , zbuf1 , zbuf2152 REAL(wp) :: ztkel, ztkel1, zt , zsal , zsal2 , zbuf1 , zbuf2 153 153 REAL(wp) :: ztgg , ztgg2, ztgg3 , ztgg4 , ztgg5 154 154 REAL(wp) :: zpres, ztc , zcl , zcpexp, zoxy , zcpexp2 … … 452 452 DO jk = 1, jpk 453 453 DO jj = 1, jpj 454 DO ji = 1, jp j454 DO ji = 1, jpi 455 455 p_alkcb = trb(ji,jj,jk,jptal) * 1000. / (rhop(ji,jj,jk) + rtrn) 456 456 p_dictot = trb(ji,jj,jk,jpdic) * 1000. / (rhop(ji,jj,jk) + rtrn) … … 533 533 REAL(wp) :: znumer_po4, zdnumer_po4, zdenom_po4, zalk_po4, zdalk_po4 534 534 REAL(wp) :: znumer_sil, zdnumer_sil, zdenom_sil, zalk_sil, zdalk_sil 535 REAL(wp) :: znumer_nh4, zdnumer_nh4, zdenom_nh4, zalk_nh4, zdalk_nh4536 REAL(wp) :: znumer_h2s, zdnumer_h2s, zdenom_h2s, zalk_h2s, zdalk_h2s537 535 REAL(wp) :: znumer_so4, zdnumer_so4, zdenom_so4, zalk_so4, zdalk_so4 538 536 REAL(wp) :: znumer_flu, zdnumer_flu, zdenom_flu, zalk_flu, zdalk_flu … … 556 554 DO jk = 1, jpk 557 555 DO jj = 1, jpj 558 DO ji = 1, jp j556 DO ji = 1, jpi 559 557 IF (rmask(ji,jj,jk) == 1.) THEN 560 558 p_alktot = trb(ji,jj,jk,jptal) * 1000. / (rhop(ji,jj,jk) + rtrn) … … 589 587 DO jk = 1, jpk 590 588 DO jj = 1, jpj 591 DO ji = 1, jp j589 DO ji = 1, jpi 592 590 IF (rmask(ji,jj,jk) == 1.) THEN 593 591 zfact = rhop(ji,jj,jk) / 1000. + rtrn … … 789 787 !! *** ROUTINE p4z_che_alloc *** 790 788 !!---------------------------------------------------------------------- 791 #if defined key_ligand792 789 INTEGER :: ierr(3) ! Local variables 793 #else794 INTEGER :: ierr(2) ! Local variables795 #endif796 790 !!---------------------------------------------------------------------- 797 791 -
branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zfechem.F90
r6453 r6841 73 73 REAL(wp) :: ztrc, zdust 74 74 #if ! defined key_kriest 75 REAL(wp) :: zdenom , zdenom275 REAL(wp) :: zdenom2 76 76 #endif 77 77 REAL(wp), POINTER, DIMENSION(:,:,:) :: zTL1, zFe3, ztotlig, precip … … 79 79 REAL(wp), POINTER, DIMENSION(:,: ) :: zstrn, zstrn2 80 80 REAL(wp) :: zzFeL1, zzFeL2, zzFe2, zzFeP, zzFe3, zzstrn2 81 REAL(wp) :: zrum, zcodel, zargu, z val, zlight81 REAL(wp) :: zrum, zcodel, zargu, zlight 82 82 REAL(wp) :: zkox, zkph1, zkph2, zph, zionic, ztligand 83 83 REAL(wp) :: za, zb, zc, zkappa1, zkappa2, za0, za1, za2 … … 172 172 zkox = zkox * MAX( 1.e-6, zoxy) / ( chemo2(ji,jj,jk) + rtrn ) 173 173 ! PHOTOREDUCTION of complexed iron : Tagliabue and Arrigo (2006) 174 zkph2 = MAX( 0., 15. * zlight / ( zlight + 2. ) ) 174 zkph2 = MAX( 0., 15. * zlight / ( zlight + 2. ) ) * (1. - fr_i(ji,jj)) 175 175 zkph1 = zkph2 / 5. 176 176 ! pass the dfe concentration from PISCES … … 325 325 ! ---------------------------------------------------------------- 326 326 zlam1a = ( 0.369 * 0.3 * trb(ji,jj,jk,jpdoc) + 102.4 * trb(ji,jj,jk,jppoc) ) * xdiss(ji,jj,jk) & 327 & + ( 114. * 0.3 * trb(ji,jj,jk,jpdoc) + 5.09E3* trb(ji,jj,jk,jppoc) )327 & + ( 114. * 0.3 * trb(ji,jj,jk,jpdoc) + 0. * trb(ji,jj,jk,jppoc) ) 328 328 zaggdfea = zlam1a * zstep * zfecoll 329 329 #if defined key_ligand -
branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zligand.F90
r6453 r6841 89 89 ! photochem loss of weak ligand 90 90 zlablgw = MAX( 0., trn(ji,jj,jk, jpfer) * plig(ji,jj,jk) ) 91 zlgwpr = prlgw * zstep * etot(ji,jj,jk) * trn(ji,jj,jk,jplgw) 91 zlgwpr = prlgw * zstep * etot(ji,jj,jk) * trn(ji,jj,jk,jplgw) * (1. - fr_i(ji,jj)) 92 92 tra(ji,jj,jk,jplgw) = tra(ji,jj,jk,jplgw) + zlgwp - zlgwr - zlgwpr 93 93 END DO … … 104 104 ! dissolution rate is maximal in the presence of light and 105 105 ! lower in the aphotici zone 106 zrfepa = rfep * (1- EXP(-1*etot(ji,jj,jk) / 25. ) ) ! 25 Wm-2 constant 106 ! ! 25 Wm-2 constant 107 zrfepa = rfep * (1- EXP(-1*etot(ji,jj,jk) / 25. ) ) * (1.- fr_i(ji,jj)) 107 108 zrfepa = MAX( (zrfepa / 10.0), zrfepa ) ! min of 10 days lifetime 108 109 zfepr = rfep * zstep * trn(ji,jj,jk,jpfep) -
branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zlim.F90
r6453 r6841 187 187 & / ( concnno3 * concnnh4 + concnnh4 * trb(ji,jj,jk,jpno3) + concnno3 * trb(ji,jj,jk,jpnh4) ) 188 188 zlim2 = trb(ji,jj,jk,jppo4) / ( trb(ji,jj,jk,jppo4) + concnnh4 ) 189 zlim3 = trb(ji,jj,jk,jpfer) / ( trb(ji,jj,jk,jpfer) + 5.E-11)189 zlim3 = trb(ji,jj,jk,jpfer) / ( trb(ji,jj,jk,jpfer) + 1.E-10 ) 190 190 ztem1 = MAX( 0., tsn(ji,jj,jk,jp_tem) ) 191 191 ztem2 = tsn(ji,jj,jk,jp_tem) - 10. -
branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zlys.F90
r6453 r6841 65 65 REAL(wp) :: zdispot, zfact, zcalcon 66 66 REAL(wp) :: zomegaca, zexcess, zexcess0 67 REAL(wp) :: zrfact268 67 CHARACTER (len=25) :: charout 69 68 REAL(wp), POINTER, DIMENSION(:,:,:) :: zco3, zcaldiss, zhinit, zhi -
branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zmeso.F90
r6453 r6841 244 244 & + zgraztotf * unass2 - zfracfe 245 245 zfracal = trb(ji,jj,jk,jpcal) / (trb(ji,jj,jk,jppoc) + trb(ji,jj,jk,jpgoc) + rtrn ) 246 zgrazcal = zgrazffeg* (1. - part2) * zfracal246 zgrazcal = (zgrazffeg + zgrazpoc) * (1. - part2) * zfracal 247 247 #endif 248 248 -
branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zmort.F90
r6453 r6841 85 85 DO jj = 1, jpj 86 86 DO ji = 1, jpi 87 zcompaph = MAX( ( trb(ji,jj,jk,jpphy) - 1e- 8), 0.e0 )87 zcompaph = MAX( ( trb(ji,jj,jk,jpphy) - 1e-9 ), 0.e0 ) 88 88 zstep = xstep 89 89 # if defined key_degrad -
branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zopt.F90
r6453 r6841 78 78 INTEGER :: irgb 79 79 REAL(wp) :: zchl 80 REAL(wp) :: z c0 , zc1 , zc2, zc3, z1_dep81 REAL(wp), POINTER, DIMENSION(:,: ) :: zdepmoy, zetmp1, zetmp2, zetmp3, zetmp4, zqsr100 80 REAL(wp) :: z1_dep 81 REAL(wp), POINTER, DIMENSION(:,: ) :: zdepmoy, zetmp1, zetmp2, zetmp3, zetmp4, zqsr100, zqsr_corr 82 82 #if defined key_pisces_quota 83 83 REAL(wp), POINTER, DIMENSION(:,: ) :: zetmp5 … … 89 89 ! 90 90 ! Allocate temporary workspace 91 CALL wrk_alloc( jpi, jpj, zqsr100, zdepmoy, zetmp1, zetmp2, zetmp3, zetmp4 )91 CALL wrk_alloc( jpi, jpj, zqsr100, zdepmoy, zetmp1, zetmp2, zetmp3, zetmp4, zqsr_corr ) 92 92 CALL wrk_alloc( jpi, jpj, jpk, zpar, ze0, ze1, ze2, ze3 ) 93 93 #if defined key_pisces_quota … … 128 128 zqsr100(:,:) = 0.01 * qsr_mean(:,:) ! daily mean qsr 129 129 ! 130 CALL p4z_opt_par( kt, qsr_mean, ze1, ze2, ze3 ) 130 zqsr_corr(:,:) = qsr_mean(:,:) / ( 1. - fr_i(:,:) + rtrn ) 131 CALL p4z_opt_par( kt, zqsr_corr, ze1, ze2, ze3 ) 131 132 ! 132 133 DO jk = 1, nksrp … … 139 140 END DO 140 141 ! 141 CALL p4z_opt_par( kt, qsr, ze1, ze2, ze3 ) 142 zqsr_corr(:,:) = qsr(:,:) / ( 1. - fr_i(:,:) + rtrn ) 143 CALL p4z_opt_par( kt, zqsr_corr, ze1, ze2, ze3 ) 142 144 ! 143 145 DO jk = 1, nksrp … … 149 151 zqsr100(:,:) = 0.01 * qsr(:,:) 150 152 ! 151 CALL p4z_opt_par( kt, qsr, ze1, ze2, ze3 ) 153 zqsr_corr(:,:) = qsr(:,:) / ( 1. - fr_i(:,:) + rtrn ) 154 CALL p4z_opt_par( kt, zqsr_corr, ze1, ze2, ze3 ) 152 155 ! 153 156 DO jk = 1, nksrp … … 174 177 ENDIF 175 178 ! !* Euphotic depth and level 176 neln(:,:) = 1 ! ------------------------ 177 heup(:,:) = 300. 179 neln(:,:) = 1 ! ------------------------ 180 heup(:,:) = fsdepw(:,:,2) 181 heup_01(:,:) = fsdepw(:,:,2) 178 182 179 183 DO jk = 2, nksrp … … 185 189 heup(ji,jj) = fsdepw(ji,jj,jk+1) ! Euphotic layer depth 186 190 ENDIF 191 IF( etot_ndcy(ji,jj,jk) * tmask(ji,jj,jk) >= 0.50 ) THEN 192 ! ! nb: ensure the compatibility with nmld_trc definition in trd_mld_trc_zint 193 heup_01(ji,jj) = fsdepw(ji,jj,jk+1) ! Euphotic layer depth 194 ENDIF 187 195 END DO 188 196 END DO 189 197 END DO 190 198 ! 191 heup(:,:) = MIN( 300., heup(:,:) ) 199 heup(:,:) = MIN( 300., heup (:,:) ) 200 heup_01(:,:) = MIN( 300., heup_01(:,:) ) 192 201 ! !* mean light over the mixed layer 193 202 zdepmoy(:,:) = 0.e0 ! ------------------------------- … … 254 263 ENDIF 255 264 ! 256 CALL wrk_dealloc( jpi, jpj, zqsr100, zdepmoy, zetmp1, zetmp2, zetmp3, zetmp4 )265 CALL wrk_dealloc( jpi, jpj, zqsr100, zdepmoy, zetmp1, zetmp2, zetmp3, zetmp4, zqsr_corr ) 257 266 CALL wrk_dealloc( jpi, jpj, jpk, zpar, ze0, ze1, ze2, ze3 ) 258 267 #if defined key_pisces_quota … … 349 358 !! * local declarations 350 359 INTEGER :: ji,jj 351 REAL(wp) :: zcoef352 360 !!--------------------------------------------------------------------- 353 361 ! -
branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zpoc.F90
r6453 r6841 30 30 PUBLIC p4z_poc ! called in p4zbio.F90 31 31 PUBLIC p4z_poc_init ! called in trcsms_pisces.F90 32 PUBLIC alngam 33 PUBLIC gamain 32 34 33 35 !! * Shared module variables 34 36 REAL(wp), PUBLIC :: xremip !: remineralisation rate of POC 35 37 INTEGER , PUBLIC :: jcpoc !: number of lability classes 36 38 REAL(wp), PUBLIC :: rshape !: shape factor of the gamma distribution 37 39 38 40 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: alphan, reminp … … 62 64 INTEGER :: ji, jj, jk, jn 63 65 REAL(wp) :: zremip, zremig, zdep, zorem2, zofer 64 REAL(wp) :: z tem, zsizek, zsizek1, alphat, remint65 REAL(wp) :: solgoc, zpoc , zpoctot, zremif66 REAL(wp) :: zsizek, zsizek1, alphat, remint 67 REAL(wp) :: solgoc, zpoc 66 68 #if ! defined key_kriest 67 69 REAL(wp) :: zofer2, zofer3 … … 100 102 DO jn = 1, jcpoc 101 103 alphag(:,:,:,jn) = alphan(jn) 104 alphap(:,:,:,jn) = alphan(jn) 102 105 END DO 103 106 … … 118 121 ! ------------------------------------------------------------ 119 122 ! 120 IF( fsdept(ji,jj,jk) > =zdep ) THEN123 IF( fsdept(ji,jj,jk) > zdep ) THEN 121 124 alphat = 0. 122 125 remint = 0. 123 !124 ! The first level just below the mixed layer needs a125 ! specific treatment because lability is supposed constant126 ! everywhere within the mixed layer. This means that127 ! change in lability in the bottom part of the previous cell128 ! should not be computed129 ! ----------------------------------------------------------130 126 ! 131 IF ( fsdept(ji,jj,jk-1) < zdep ) THEN 127 IF ( fsdept(ji,jj,jk-1) <= zdep ) THEN 128 ! 129 ! The first level just below the mixed layer needs a 130 ! specific treatment because lability is supposed constant 131 ! everywhere within the mixed layer. This means that 132 ! change in lability in the bottom part of the previous cell 133 ! should not be computed 134 ! ---------------------------------------------------------- 135 ! 136 ! POC concentration is computed using the lagrangian 137 ! framework. It is only used for the lability param 138 zpoc = trb(ji,jj,jk-1,jpgoc) + consgoc(ji,jj,jk) * rday / rfact2 & 139 & * fse3t(ji,jj,jk) / 2. / (wsbio4(ji,jj,jk) + rtrn) 140 zpoc = max(0., zpoc) 141 ! 132 142 DO jn = 1, jcpoc 133 143 ! … … 139 149 zsizek = zdep / (wsbio2 + rtrn) * tgfunc(ji,jj,jk-1) 140 150 zsizek1 = fse3t(ji,jj,jk) / 2. / (wsbio4(ji,jj,jk) + rtrn) * tgfunc(ji,jj,jk) 141 ! POC concentration is computed using the lagrangian142 ! framework. It is only used for the lability param143 zpoc = trb(ji,jj,jk-1,jpgoc) + consgoc(ji,jj,jk) * rday / rfact2 &144 & * fse3t(ji,jj,jk) / 2. / (wsbio4(ji,jj,jk) + rtrn)145 zpoc = max(0., zpoc)146 151 ! the concentration of each lability class is calculated 147 152 ! as the sum of the different sources and sinks … … 161 166 ! --------------------------------------------------- 162 167 ! 168 zpoc = trb(ji,jj,jk-1,jpgoc) + consgoc(ji,jj,jk-1) * rday / rfact2 & 169 & * fse3t(ji,jj,jk-1) / 2. / (wsbio4(ji,jj,jk-1) + rtrn) + consgoc(ji,jj,jk) & 170 & * rday / rfact2 * fse3t(ji,jj,jk) / 2. / (wsbio4(ji,jj,jk) + rtrn) 171 zpoc = max(0., zpoc) 172 ! 163 173 DO jn = 1, jcpoc 164 174 zsizek = fse3t(ji,jj,jk-1) / 2. / (wsbio4(ji,jj,jk-1) + rtrn) * tgfunc(ji,jj,jk-1) 165 175 zsizek1 = fse3t(ji,jj,jk) / 2. / (wsbio4(ji,jj,jk) + rtrn) * tgfunc(ji,jj,jk) 166 zpoc = trb(ji,jj,jk-1,jpgoc) + consgoc(ji,jj,jk-1) * rday / rfact2 &167 & * fse3t(ji,jj,jk-1) / 2. / (wsbio4(ji,jj,jk-1) + rtrn) + consgoc(ji,jj,jk) &168 & * rday / rfact2 * fse3t(ji,jj,jk) / 2. / (wsbio4(ji,jj,jk) + rtrn)169 zpoc = max(0., zpoc)170 176 alphag(ji,jj,jk,jn) = alphag(ji,jj,jk-1,jn) * exp( -reminp(jn) & 171 177 & * ( zsizek + zsizek1 ) ) * zpoc + ( prodgoc(ji,jj,jk-1) & … … 273 279 & * totthick(ji,jj) + totcons(ji,jj) + wsbio + rtrn ) 274 280 alphat = alphat + alphap(ji,jj,jk,jn) 275 remint = remint + alphap(ji,jj,jk,jn) * reminp(jn)276 281 END DO 277 282 DO jn = 1, jcpoc 278 283 alphap(ji,jj,jk,jn) = alphap(ji,jj,jk,jn) / ( alphat + rtrn) 284 remint = remint + alphap(ji,jj,jk,jn) * reminp(jn) 279 285 END DO 280 286 ! Mean remineralization rate in the mixed layer 281 zremipoc(ji,jj,jk) = MIN(xremip, MAX(0., remint / ( alphat + rtrn)))287 zremipoc(ji,jj,jk) = MIN(xremip, MAX(0., remint )) 282 288 ENDIF 283 289 ENDIF … … 308 314 ! 309 315 IF ( fsdept(ji,jj,jk-1) <= zdep ) THEN 316 ! 317 ! Computation of the POC concentration using the 318 ! lagrangian algorithm 319 zpoc = trb(ji,jj,jk-1,jppoc) + conspoc(ji,jj,jk) * rday / rfact2 & 320 & * fse3t(ji,jj,jk) / 2. / (wsbio3(ji,jj,jk) + rtrn) 321 zpoc = max(0., zpoc) 322 ! 310 323 DO jn = 1, jcpoc 311 324 ! the scale factor is corrected with temperature 312 325 zsizek1 = fse3t(ji,jj,jk) / 2. / (wsbio3(ji,jj,jk) + rtrn) * tgfunc(ji,jj,jk) 313 ! Computation of the POC concentration using the314 ! lagrangian algorithm315 zpoc = trb(ji,jj,jk-1,jppoc) + conspoc(ji,jj,jk) * rday / rfact2 &316 & * fse3t(ji,jj,jk) / 2. / (wsbio3(ji,jj,jk) + rtrn)317 zpoc = max(0., zpoc)318 326 ! computation of the lability spectrum applying the 319 327 ! different sources and sinks … … 321 329 & * ( 1. - exp( -reminp(jn) * zsizek1 ) ) * rday / rfact2 & 322 330 & * alphag(ji,jj,jk,jn) / reminp(jn) / tgfunc(ji,jj,jk) 331 alphap(ji,jj,jk,jn) = MAX( 0., alphap(ji,jj,jk,jn) ) 323 332 alphat = alphat + alphap(ji,jj,jk,jn) 324 remint = remint + alphap(ji,jj,jk,jn) * reminp(jn)325 333 END DO 326 334 ELSE … … 331 339 ! -------------------------------------------------------- 332 340 ! 341 zpoc = trb(ji,jj,jk-1,jppoc) + conspoc(ji,jj,jk-1) * rday / rfact2 & 342 & * fse3t(ji,jj,jk-1) / 2. / (wsbio3(ji,jj,jk-1) + rtrn) + conspoc(ji,jj,jk) & 343 & * rday / rfact2 * fse3t(ji,jj,jk) / 2. / (wsbio3(ji,jj,jk) + rtrn) 344 zpoc = max(0., zpoc) 345 ! 333 346 DO jn = 1, jcpoc 334 347 zsizek = fse3t(ji,jj,jk-1) / 2. / (wsbio3(ji,jj,jk-1) + rtrn) * tgfunc(ji,jj,jk-1) 335 348 zsizek1 = fse3t(ji,jj,jk) / 2. / (wsbio3(ji,jj,jk) + rtrn) * tgfunc(ji,jj,jk) 336 zpoc = trb(ji,jj,jk-1,jppoc) + conspoc(ji,jj,jk-1) * rday / rfact2 &337 & * fse3t(ji,jj,jk-1) / 2. / (wsbio3(ji,jj,jk-1) + rtrn) + conspoc(ji,jj,jk) &338 & * rday / rfact2 * fse3t(ji,jj,jk) / 2. / (wsbio3(ji,jj,jk) + rtrn)339 zpoc = max(0., zpoc)340 349 alphap(ji,jj,jk,jn) = alphap(ji,jj,jk-1,jn) * exp( -reminp(jn) & 341 350 & * ( zsizek + zsizek1 ) ) * zpoc + ( prodpoc(ji,jj,jk-1) & … … 349 358 & + zorem3(ji,jj,jk) * ( 1. - exp( -reminp(jn) * zsizek1 ) ) * rday & 350 359 & / rfact2 * alphag(ji,jj,jk,jn) / reminp(jn) / tgfunc(ji,jj,jk) 360 alphap(ji,jj,jk,jn) = max(0., alphap(ji,jj,jk,jn) ) 351 361 alphat = alphat + alphap(ji,jj,jk,jn) 352 remint = remint + alphap(ji,jj,jk,jn) * reminp(jn)353 362 END DO 354 363 ENDIF … … 357 366 DO jn = 1, jcpoc 358 367 alphap(ji,jj,jk,jn) = alphap(ji,jj,jk,jn) / ( alphat + rtrn) 368 remint = remint + alphap(ji,jj,jk,jn) * reminp(jn) 359 369 END DO 360 370 ! Mean remineralization rate in the water column 361 zremipoc(ji,jj,jk) = MIN(xremip, MAX(0., remint / ( alphat + rtrn)))371 zremipoc(ji,jj,jk) = MIN(xremip, MAX(0., remint )) 362 372 ENDIF 363 373 ENDIF … … 436 446 INTEGER :: jn 437 447 REAL(wp) :: remindelta, reminup, remindown 438 NAMELIST/nampispoc/ xremip, jcpoc 448 INTEGER :: ifault 449 450 NAMELIST/nampispoc/ xremip, jcpoc, rshape 439 451 INTEGER :: ios ! Local integer output status for namelist read 440 452 … … 454 466 WRITE(numout,*) ' remineralisation rate of POC xremip =', xremip 455 467 WRITE(numout,*) ' Number of lability classes for POC jcpoc =', jcpoc 468 WRITE(numout,*) ' Shape factor of the gamma distribution rshape =', rshape 456 469 ENDIF 457 470 ! … … 463 476 ! 464 477 IF (jcpoc > 1) THEN 478 ! 465 479 remindelta = log(4. * 1000. ) / float(jcpoc-1) 466 reminup = xremip/400. * exp(remindelta) 467 alphan(1) = 1. - exp(-reminup/xremip) 468 reminp(1) = xremip - reminup * exp(-reminup/xremip) / alphan(1) 480 reminup = 1./ 400. * exp(remindelta) 481 ! 482 ! Discretization based on incomplete gamma functions 483 ! As incomplete gamma functions are not available in standard 484 ! fortran 95, they have been coded as functions in this module (gamain) 485 ! --------------------------------------------------------------------- 486 ! 487 alphan(1) = gamain(reminup, rshape, ifault) 488 reminp(1) = gamain(reminup, rshape+1.0, ifault) * xremip / alphan(1) 469 489 DO jn = 2, jcpoc-1 470 reminup = xremip/400. * exp(float(jn) * remindelta)471 remindown = xremip/ 400. * exp(float(jn-1) * remindelta)472 alphan(jn) = exp(-remindown /xremip) - exp(-reminup/xremip)473 reminp(jn) = xremip + (remindown * exp(-remindown /xremip) &474 & - reminup * exp(-reminup/xremip) )/ alphan(jn)490 reminup = 1./ 400. * exp(float(jn) * remindelta) 491 remindown = 1. / 400. * exp(float(jn-1) * remindelta) 492 alphan(jn) = gamain(reminup, rshape, ifault) - gamain(remindown, rshape, ifault) 493 reminp(jn) = gamain(reminup, rshape+1.0, ifault) - gamain(remindown, rshape+1.0, ifault) 494 reminp(jn) = reminp(jn) * xremip / alphan(jn) 475 495 END DO 476 remindown = xremip/400. * exp(float(jcpoc-1) * remindelta) 477 alphan(jcpoc) = exp(-remindown /xremip) 478 reminp(jcpoc) = xremip + remindown 496 remindown = 1. / 400. * exp(float(jcpoc-1) * remindelta) 497 alphan(jcpoc) = 1.0 - gamain(remindown, rshape, ifault) 498 reminp(jcpoc) = 1.0 - gamain(remindown, rshape+1.0, ifault) 499 reminp(jcpoc) = reminp(jcpoc) * xremip / alphan(jcpoc) 500 479 501 ELSE 480 502 alphan(jcpoc) = 1. … … 487 509 488 510 END SUBROUTINE p4z_poc_init 511 512 REAL FUNCTION alngam( xvalue, ifault ) 513 514 !*****************************************************************************80 515 ! 516 !! ALNGAM computes the logarithm of the gamma function. 517 ! 518 ! Modified: 519 ! 520 ! 13 January 2008 521 ! 522 ! Author: 523 ! 524 ! Allan Macleod 525 ! FORTRAN90 version by John Burkardt 526 ! 527 ! Reference: 528 ! 529 ! Allan Macleod, 530 ! Algorithm AS 245, 531 ! A Robust and Reliable Algorithm for the Logarithm of the Gamma Function, 532 ! Applied Statistics, 533 ! Volume 38, Number 2, 1989, pages 397-402. 534 ! 535 ! Parameters: 536 ! 537 ! Input, real ( kind = 8 ) XVALUE, the argument of the Gamma function. 538 ! 539 ! Output, integer ( kind = 4 ) IFAULT, error flag. 540 ! 0, no error occurred. 541 ! 1, XVALUE is less than or equal to 0. 542 ! 2, XVALUE is too big. 543 ! 544 ! Output, real ( kind = 8 ) ALNGAM, the logarithm of the gamma function of X. 545 ! 546 implicit none 547 548 real(wp), parameter :: alr2pi = 0.918938533204673E+00 549 integer:: ifault 550 real(wp), dimension ( 9 ) :: r1 = (/ & 551 -2.66685511495E+00, & 552 -24.4387534237E+00, & 553 -21.9698958928E+00, & 554 11.1667541262E+00, & 555 3.13060547623E+00, & 556 0.607771387771E+00, & 557 11.9400905721E+00, & 558 31.4690115749E+00, & 559 15.2346874070E+00 /) 560 real(wp), dimension ( 9 ) :: r2 = (/ & 561 -78.3359299449E+00, & 562 -142.046296688E+00, & 563 137.519416416E+00, & 564 78.6994924154E+00, & 565 4.16438922228E+00, & 566 47.0668766060E+00, & 567 313.399215894E+00, & 568 263.505074721E+00, & 569 43.3400022514E+00 /) 570 real(wp), dimension ( 9 ) :: r3 = (/ & 571 -2.12159572323E+05, & 572 2.30661510616E+05, & 573 2.74647644705E+04, & 574 -4.02621119975E+04, & 575 -2.29660729780E+03, & 576 -1.16328495004E+05, & 577 -1.46025937511E+05, & 578 -2.42357409629E+04, & 579 -5.70691009324E+02 /) 580 real(wp), dimension ( 5 ) :: r4 = (/ & 581 0.279195317918525E+00, & 582 0.4917317610505968E+00, & 583 0.0692910599291889E+00, & 584 3.350343815022304E+00, & 585 6.012459259764103E+00 /) 586 real (wp) :: x 587 real (wp) :: x1 588 real (wp) :: x2 589 real (wp), parameter :: xlge = 5.10E+05 590 real (wp), parameter :: xlgst = 1.0E+30 591 real (wp) :: xvalue 592 real (wp) :: y 593 594 x = xvalue 595 alngam = 0.0E+00 596 ! 597 ! Check the input. 598 ! 599 if ( xlgst <= x ) then 600 ifault = 2 601 return 602 end if 603 if ( x <= 0.0E+00 ) then 604 ifault = 1 605 return 606 end if 607 608 ifault = 0 609 ! 610 ! Calculation for 0 < X < 0.5 and 0.5 <= X < 1.5 combined. 611 ! 612 if ( x < 1.5E+00 ) then 613 614 if ( x < 0.5E+00 ) then 615 alngam = - log ( x ) 616 y = x + 1.0E+00 617 ! 618 ! Test whether X < machine epsilon. 619 ! 620 if ( y == 1.0E+00 ) then 621 return 622 end if 623 624 else 625 626 alngam = 0.0E+00 627 y = x 628 x = ( x - 0.5E+00 ) - 0.5E+00 629 630 end if 631 632 alngam = alngam + x * (((( & 633 r1(5) * y & 634 + r1(4) ) * y & 635 + r1(3) ) * y & 636 + r1(2) ) * y & 637 + r1(1) ) / (((( & 638 y & 639 + r1(9) ) * y & 640 + r1(8) ) * y & 641 + r1(7) ) * y & 642 + r1(6) ) 643 644 return 645 646 end if 647 ! 648 ! Calculation for 1.5 <= X < 4.0. 649 ! 650 if ( x < 4.0E+00 ) then 651 652 y = ( x - 1.0E+00 ) - 1.0E+00 653 654 alngam = y * (((( & 655 r2(5) * x & 656 + r2(4) ) * x & 657 + r2(3) ) * x & 658 + r2(2) ) * x & 659 + r2(1) ) / (((( & 660 x & 661 + r2(9) ) * x & 662 + r2(8) ) * x & 663 + r2(7) ) * x & 664 + r2(6) ) 665 ! 666 ! Calculation for 4.0 <= X < 12.0. 667 ! 668 else if ( x < 12.0E+00 ) then 669 670 alngam = (((( & 671 r3(5) * x & 672 + r3(4) ) * x & 673 + r3(3) ) * x & 674 + r3(2) ) * x & 675 + r3(1) ) / (((( & 676 x & 677 + r3(9) ) * x & 678 + r3(8) ) * x & 679 + r3(7) ) * x & 680 + r3(6) ) 681 ! 682 ! Calculation for 12.0 <= X. 683 ! 684 else 685 686 y = log ( x ) 687 alngam = x * ( y - 1.0E+00 ) - 0.5E+00 * y + alr2pi 688 689 if ( x <= xlge ) then 690 691 x1 = 1.0E+00 / x 692 x2 = x1 * x1 693 694 alngam = alngam + x1 * ( ( & 695 r4(3) * & 696 x2 + r4(2) ) * & 697 x2 + r4(1) ) / ( ( & 698 x2 + r4(5) ) * & 699 x2 + r4(4) ) 700 701 end if 702 703 end if 704 705 END FUNCTION alngam 706 707 REAL FUNCTION gamain( x, p, ifault ) 708 709 !*****************************************************************************80 710 ! 711 !! GAMAIN computes the incomplete gamma ratio. 712 ! 713 ! Discussion: 714 ! 715 ! A series expansion is used if P > X or X <= 1. Otherwise, a 716 ! continued fraction approximation is used. 717 ! 718 ! Modified: 719 ! 720 ! 17 January 2008 721 ! 722 ! Author: 723 ! 724 ! G Bhattacharjee 725 ! FORTRAN90 version by John Burkardt 726 ! 727 ! Reference: 728 ! 729 ! G Bhattacharjee, 730 ! Algorithm AS 32: 731 ! The Incomplete Gamma Integral, 732 ! Applied Statistics, 733 ! Volume 19, Number 3, 1970, pages 285-287. 734 ! 735 ! Parameters: 736 ! 737 ! Input, real ( kind = 8 ) X, P, the parameters of the incomplete 738 ! gamma ratio. 0 <= X, and 0 < P. 739 ! 740 ! Output, integer ( kind = 4 ) IFAULT, error flag. 741 ! 0, no errors. 742 ! 1, P <= 0. 743 ! 2, X < 0. 744 ! 3, underflow. 745 ! 4, error return from the Log Gamma routine. 746 ! 747 ! Output, real ( kind = 8 ) GAMAIN, the value of the incomplete 748 ! gamma ratio. 749 ! 750 implicit none 751 752 real (wp) a 753 real (wp), parameter :: acu = 1.0E-08 754 real (wp) an 755 real (wp) arg 756 real (wp) b 757 real (wp) dif 758 real (wp) factor 759 real (wp) g 760 real (wp) gin 761 integer i 762 integer ifault 763 real (wp), parameter :: oflo = 1.0E+37 764 real (wp) p 765 real (wp) pn(6) 766 real (wp) rn 767 real (wp) term 768 real (wp), parameter :: uflo = 1.0E-37 769 real (wp) x 770 ! 771 ! Check the input. 772 ! 773 if ( p <= 0.0E+00 ) then 774 ifault = 1 775 gamain = 0.0E+00 776 return 777 end if 778 779 if ( x < 0.0E+00 ) then 780 ifault = 2 781 gamain = 0.0E+00 782 return 783 end if 784 785 if ( x == 0.0E+00 ) then 786 ifault = 0 787 gamain = 0.0E+00 788 return 789 end if 790 791 g = alngam ( p, ifault ) 792 793 if ( ifault /= 0 ) then 794 ifault = 4 795 gamain = 0.0E+00 796 return 797 end if 798 799 arg = p * log ( x ) - x - g 800 801 if ( arg < log ( uflo ) ) then 802 ifault = 3 803 gamain = 0.0E+00 804 return 805 end if 806 807 ifault = 0 808 factor = exp ( arg ) 809 ! 810 ! Calculation by series expansion. 811 ! 812 if ( x <= 1.0E+00 .or. x < p ) then 813 814 gin = 1.0E+00 815 term = 1.0E+00 816 rn = p 817 818 do 819 820 rn = rn + 1.0E+00 821 term = term * x / rn 822 gin = gin + term 823 824 if ( term <= acu ) then 825 exit 826 end if 827 828 end do 829 830 gamain = gin * factor / p 831 return 832 833 end if 834 ! 835 ! Calculation by continued fraction. 836 ! 837 a = 1.0E+00 - p 838 b = a + x + 1.0E+00 839 term = 0.0E+00 840 841 pn(1) = 1.0E+00 842 pn(2) = x 843 pn(3) = x + 1.0E+00 844 pn(4) = x * b 845 846 gin = pn(3) / pn(4) 847 848 do 849 850 a = a + 1.0E+00 851 b = b + 2.0E+00 852 term = term + 1.0E+00 853 an = a * term 854 do i = 1, 2 855 pn(i+4) = b * pn(i+2) - an * pn(i) 856 end do 857 858 if ( pn(6) /= 0.0E+00 ) then 859 860 rn = pn(5) / pn(6) 861 dif = abs ( gin - rn ) 862 ! 863 ! Absolute error tolerance satisfied? 864 ! 865 if ( dif <= acu ) then 866 ! 867 ! Relative error tolerance satisfied? 868 ! 869 if ( dif <= acu * rn ) then 870 gamain = 1.0E+00 - factor * gin 871 exit 872 end if 873 874 end if 875 876 gin = rn 877 878 end if 879 880 do i = 1, 4 881 pn(i) = pn(i+2) 882 end do 883 if ( oflo <= abs ( pn(5) ) ) then 884 885 do i = 1, 4 886 pn(i) = pn(i) / oflo 887 end do 888 889 end if 890 891 end do 892 893 END FUNCTION gamain 489 894 490 895 #else -
branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zprod.F90
r6453 r6841 85 85 REAL(wp) :: zfact 86 86 CHARACTER (len=25) :: charout 87 REAL(wp), POINTER, DIMENSION(:,: ) :: z mixnano, zmixdiat, zstrn, zw2d87 REAL(wp), POINTER, DIMENSION(:,: ) :: zstrn, zw2d 88 88 REAL(wp), POINTER, DIMENSION(:,:,:) :: zpislopead, zpislopead2, zprdia, zprbio, zprdch, zprnch, zysopt, zw3d 89 89 REAL(wp), POINTER, DIMENSION(:,:,:) :: zprorca, zprorcad, zprofed, zprofen, zprochln, zprochld, zpronew, zpronewd 90 REAL(wp), POINTER, DIMENSION(:,:,:) :: zmxl_fac, zmxl_chl 90 91 !!--------------------------------------------------------------------- 91 92 ! … … 93 94 ! 94 95 ! Allocate temporary workspace 95 CALL wrk_alloc( jpi, jpj, zmixnano, zmixdiat, zstrn ) 96 CALL wrk_alloc( jpi, jpj, jpk, zpislopead, zpislopead2, zprdia, zprbio, zprdch, zprnch, zysopt ) 96 CALL wrk_alloc( jpi, jpj, zstrn ) 97 CALL wrk_alloc( jpi, jpj, jpk, zpislopead, zpislopead2, zprdia, zprbio, zprdch, zprnch, zysopt ) 98 CALL wrk_alloc( jpi, jpj, jpk, zmxl_fac, zmxl_chl ) 97 99 CALL wrk_alloc( jpi, jpj, jpk, zprorca, zprorcad, zprofed, zprofen, zprochln, zprochld, zpronew, zpronewd ) 98 100 ! … … 129 131 END DO 130 132 131 ! Impact of the day duration on phytoplankton growth133 ! Impact of the day duration and light intermittency on phytoplankton growth 132 134 DO jk = 1, jpkm1 133 135 DO jj = 1 ,jpj … … 135 137 IF( etot_ndcy(ji,jj,jk) > 1.E-3 ) THEN 136 138 zval = MAX( 1., zstrn(ji,jj) ) 137 zval = 1.5 * zval / ( 12. + zval ) 138 zprbio(ji,jj,jk) = prmax(ji,jj,jk) * zval 139 zprdia(ji,jj,jk) = zprbio(ji,jj,jk) 139 ! zval = 1.5 * zval / ( 12. + zval ) 140 ! zmxl_fac(ji,jj,jk) = zval * ( 1. - fr_i(ji,jj)) 141 zmxl_fac(ji,jj,jk) = zval 142 zmxl_chl(ji,jj,jk) = zval / 24. * ( 1. - fr_i(ji,jj)) 143 IF( fsdepw(ji,jj,jk+1) <= hmld(ji,jj) ) THEN 144 zval = MIN(1., heup_01(ji,jj) / ( hmld(ji,jj) + rtrn )) 145 zmxl_fac(ji,jj,jk) = zmxl_fac(ji,jj,jk) * zval 146 zmxl_chl(ji,jj,jk) = zmxl_chl(ji,jj,jk) * zval 147 ENDIF 148 IF (mig(ji) == 60 .and. mjg(jj) == 25 .and. jk == 1) THEn 149 write(0,*) 'plante ',zmxl_fac(ji,jj,jk),MAX( 1., zstrn(ji,jj) ), & 150 & -0.2 * zmxl_fac(ji,jj,jk) 151 ENDIF 152 zmxl_fac(ji,jj,jk) = ( 1. - exp( -0.2 * zmxl_fac(ji,jj,jk) ) ) * ( 1. - fr_i(ji,jj) ) 153 zmxl_chl(ji,jj,jk) = zmxl_chl(ji,jj,jk) * ( 1. - fr_i(ji,jj) ) 140 154 ENDIF 141 155 END DO … … 143 157 END DO 144 158 145 ! Maximum light intensity 146 WHERE( zstrn(:,:) < 1.e0 ) zstrn(:,:) = 24. 147 zstrn(:,:) = 24. / zstrn(:,:) 159 zprbio(:,:,:) = prmax(:,:,:) * zmxl_fac(:,:,:) 160 zprdia(:,:,:) = zprbio(:,:,:) 161 162 ! Computation of the P-I slope for nanos and diatoms 163 DO jk = 1, jpkm1 164 !CDIR NOVERRCHK 165 DO jj = 1, jpj 166 !CDIR NOVERRCHK 167 DO ji = 1, jpi 168 IF( etot_ndcy(ji,jj,jk) > 1.E-3 ) THEN 169 ztn = MAX( 0., tsn(ji,jj,jk,jp_tem) - 15. ) 170 zadap = xadap * ztn / ( 2.+ ztn ) 171 zconctemp = MAX( 0.e0 , trb(ji,jj,jk,jpdia) - xsizedia ) 172 zconctemp2 = trb(ji,jj,jk,jpdia) - zconctemp 173 ! 174 zpislopead (ji,jj,jk) = pislope * ( 1.+ zadap * EXP( -0.25 * enano(ji,jj,jk) ) ) & 175 & * trb(ji,jj,jk,jpnch) /( trb(ji,jj,jk,jpphy) * 12. + rtrn) 176 ! 177 zpislopead2(ji,jj,jk) = (pislope * zconctemp2 + pislope2 * zconctemp) / ( trb(ji,jj,jk,jpdia) + rtrn ) & 178 & * trb(ji,jj,jk,jpdch) /( trb(ji,jj,jk,jpdia) * 12. + rtrn) 179 ENDIF 180 END DO 181 END DO 182 END DO 148 183 149 184 IF( ln_newprod ) THEN 150 !CDIR NOVERRCHK151 185 DO jk = 1, jpkm1 152 !CDIR NOVERRCHK153 186 DO jj = 1, jpj 154 !CDIR NOVERRCHK155 187 DO ji = 1, jpi 156 ! Computation of the P-I slope for nanos and diatoms157 188 IF( etot_ndcy(ji,jj,jk) > 1.E-3 ) THEN 158 ztn = MAX( 0., tsn(ji,jj,jk,jp_tem) - 15. )159 zadap = xadap * ztn / ( 2.+ ztn )160 zconctemp = MAX( 0.e0 , trb(ji,jj,jk,jpdia) - xsizedia )161 zconctemp2 = trb(ji,jj,jk,jpdia) - zconctemp162 znanotot = enano(ji,jj,jk) * zstrn(ji,jj)163 zdiattot = ediat(ji,jj,jk) * zstrn(ji,jj)164 !165 zpislopead (ji,jj,jk) = pislope * ( 1.+ zadap * EXP( -znanotot ) ) &166 & * trb(ji,jj,jk,jpnch) /( trb(ji,jj,jk,jpphy) * 12. + rtrn)167 !168 zpislopead2(ji,jj,jk) = (pislope * zconctemp2 + pislope2 * zconctemp) / ( trb(ji,jj,jk,jpdia) + rtrn ) &169 & * trb(ji,jj,jk,jpdch) /( trb(ji,jj,jk,jpdia) * 12. + rtrn)170 171 189 ! Computation of production function for Carbon 172 190 ! --------------------------------------------- 173 zpislopen = zpislopead (ji,jj,jk) / ( ( r1_rday + bresp * r1_rday ) * rday + rtrn) 174 zpislope2n = zpislopead2(ji,jj,jk) / ( ( r1_rday + bresp * r1_rday ) * rday + rtrn) 175 zprbio(ji,jj,jk) = zprbio(ji,jj,jk) * ( 1.- EXP( -zpislopen * znanotot ) ) 176 zprdia(ji,jj,jk) = zprdia(ji,jj,jk) * ( 1.- EXP( -zpislope2n * zdiattot ) ) 177 191 zpislopen = zpislopead (ji,jj,jk) / ( ( r1_rday + bresp * r1_rday ) & 192 & * zmxl_fac(ji,jj,jk) * rday + rtrn) 193 zpislope2n = zpislopead2(ji,jj,jk) / ( ( r1_rday + bresp * r1_rday ) & 194 & * zmxl_fac(ji,jj,jk) * rday + rtrn) 195 zprbio(ji,jj,jk) = zprbio(ji,jj,jk) * ( 1.- EXP( -zpislopen * enano(ji,jj,jk) ) ) 196 zprdia(ji,jj,jk) = zprdia(ji,jj,jk) * ( 1.- EXP( -zpislope2n * ediat(ji,jj,jk) ) ) 178 197 ! Computation of production function for Chlorophyll 179 198 !-------------------------------------------------- 180 zmaxday = 1._wp / ( prmax(ji,jj,jk) * rday + rtrn ) 181 zprnch(ji,jj,jk) = prmax(ji,jj,jk) * ( 1.- EXP( -zpislopead (ji,jj,jk) * zmaxday * znanotot ) ) 182 zprdch(ji,jj,jk) = prmax(ji,jj,jk) * ( 1.- EXP( -zpislopead2(ji,jj,jk) * zmaxday * zdiattot ) ) 183 ENDIF 184 END DO 185 END DO 186 END DO 187 ELSE 188 !CDIR NOVERRCHK 189 DO jk = 1, jpkm1 190 !CDIR NOVERRCHK 191 DO jj = 1, jpj 192 !CDIR NOVERRCHK 193 DO ji = 1, jpi 194 195 ! Computation of the P-I slope for nanos and diatoms 196 IF( etot_ndcy(ji,jj,jk) > 1.E-3 ) THEN 197 ztn = MAX( 0., tsn(ji,jj,jk,jp_tem) - 15. ) 198 zadap = ztn / ( 2.+ ztn ) 199 zconctemp = MAX( 0.e0 , trb(ji,jj,jk,jpdia) - xsizedia ) 200 zconctemp2 = trb(ji,jj,jk,jpdia) - zconctemp 201 znanotot = enano(ji,jj,jk) * zstrn(ji,jj) 202 zdiattot = ediat(ji,jj,jk) * zstrn(ji,jj) 203 ! 204 zpislopead (ji,jj,jk) = pislope * ( 1.+ zadap * EXP( -znanotot ) ) 205 zpislopead2(ji,jj,jk) = (pislope * zconctemp2 + pislope2 * zconctemp) / ( trb(ji,jj,jk,jpdia) + rtrn ) 206 207 zpislopen = zpislopead(ji,jj,jk) * trb(ji,jj,jk,jpnch) & 208 & / ( trb(ji,jj,jk,jpphy) * 12. + rtrn ) & 209 & / ( prmax(ji,jj,jk) * rday * xlimphy(ji,jj,jk) + rtrn ) 210 211 zpislope2n = zpislopead2(ji,jj,jk) * trb(ji,jj,jk,jpdch) & 212 & / ( trb(ji,jj,jk,jpdia) * 12. + rtrn ) & 213 & / ( prmax(ji,jj,jk) * rday * xlimdia(ji,jj,jk) + rtrn ) 214 215 ! Computation of production function for Carbon 216 ! --------------------------------------------- 217 zprbio(ji,jj,jk) = zprbio(ji,jj,jk) * ( 1.- EXP( -zpislopen * znanotot ) ) 218 zprdia(ji,jj,jk) = zprdia(ji,jj,jk) * ( 1.- EXP( -zpislope2n * zdiattot ) ) 219 220 ! Computation of production function for Chlorophyll 221 !-------------------------------------------------- 199 zpislopen = zpislopead (ji,jj,jk) / ( prmax(ji,jj,jk) * zmxl_chl(ji,jj,jk) * rday + rtrn ) 200 zpislope2n = zpislopead2(ji,jj,jk) / ( prmax(ji,jj,jk) * zmxl_chl(ji,jj,jk) * rday + rtrn ) 222 201 zprnch(ji,jj,jk) = prmax(ji,jj,jk) * ( 1.- EXP( -zpislopen * enano(ji,jj,jk) ) ) 223 202 zprdch(ji,jj,jk) = prmax(ji,jj,jk) * ( 1.- EXP( -zpislope2n * ediat(ji,jj,jk) ) ) … … 226 205 END DO 227 206 END DO 207 ELSE 208 DO jk = 1, jpkm1 209 DO jj = 1, jpj 210 DO ji = 1, jpi 211 IF( etot_ndcy(ji,jj,jk) > 1.E-3 ) THEN 212 ! Computation of production function for Carbon 213 ! --------------------------------------------- 214 zpislopen = zpislopead(ji,jj,jk) / ( zprbio(ji,jj,jk) * rday * xlimphy(ji,jj,jk) + rtrn ) 215 zpislope2n = zpislopead2(ji,jj,jk) / ( zprdia(ji,jj,jk) * rday * xlimdia(ji,jj,jk) + rtrn ) 216 zprbio(ji,jj,jk) = zprbio(ji,jj,jk) * ( 1.- EXP( -zpislopen * enano(ji,jj,jk) ) ) 217 zprdia(ji,jj,jk) = zprdia(ji,jj,jk) * ( 1.- EXP( -zpislope2n * ediat(ji,jj,jk) ) ) 218 ! Computation of production function for Chlorophyll 219 !-------------------------------------------------- 220 zpislopen = zpislopen * zmxl_fac(ji,jj,jk) / ( zmxl_chl(ji,jj,jk) + rtrn ) 221 zpislope2n = zpislope2n * zmxl_fac(ji,jj,jk) / ( zmxl_chl(ji,jj,jk) + rtrn ) 222 zprnch(ji,jj,jk) = prmax(ji,jj,jk) * ( 1.- EXP( -zpislopen * enano(ji,jj,jk) ) ) 223 zprdch(ji,jj,jk) = prmax(ji,jj,jk) * ( 1.- EXP( -zpislope2n * ediat(ji,jj,jk) ) ) 224 ENDIF 225 END DO 226 END DO 227 END DO 228 228 ENDIF 229 229 … … 231 231 ! Computation of a proxy of the N/C ratio 232 232 ! --------------------------------------- 233 !CDIR NOVERRCHK 234 DO jk = 1, jpkm1 235 !CDIR NOVERRCHK 233 DO jk = 1, jpkm1 236 234 DO jj = 1, jpj 237 !CDIR NOVERRCHK238 235 DO ji = 1, jpi 239 236 zval = MIN( xnanopo4(ji,jj,jk), ( xnanonh4(ji,jj,jk) + xnanono3(ji,jj,jk) ) ) & … … 269 266 zysopt(ji,jj,jk) = grosip * zlim * zsilfac * zsilfac2 270 267 ENDIF 271 END DO272 END DO273 END DO274 275 ! Computation of the limitation term due to a mixed layer deeper than the euphotic depth276 DO jj = 1, jpj277 DO ji = 1, jpi278 zmxltst = MAX( 0.e0, hmld(ji,jj) - heup(ji,jj) )279 zmxlday = zmxltst * zmxltst * r1_rday280 zmixnano(ji,jj) = 1. - zmxlday / ( 2. + zmxlday )281 zmixdiat(ji,jj) = 1. - zmxlday / ( 4. + zmxlday )282 END DO283 END DO284 285 ! Mixed-layer effect on production286 DO jk = 1, jpkm1287 DO jj = 1, jpj288 DO ji = 1, jpi289 IF( fsdepw(ji,jj,jk+1) <= hmld(ji,jj) ) THEN290 zprbio(ji,jj,jk) = zprbio(ji,jj,jk) * zmixnano(ji,jj)291 zprdia(ji,jj,jk) = zprdia(ji,jj,jk) * zmixdiat(ji,jj)292 ENDIF293 268 END DO 294 269 END DO … … 330 305 END DO 331 306 332 IF( ln_newprod ) THEN 333 !CDIR NOVERRCHK 334 DO jk = 1, jpkm1 335 !CDIR NOVERRCHK 336 DO jj = 1, jpj 337 !CDIR NOVERRCHK 338 DO ji = 1, jpi 339 IF( fsdepw(ji,jj,jk+1) <= hmld(ji,jj) ) THEN 340 zprnch(ji,jj,jk) = zprnch(ji,jj,jk) * zmixnano(ji,jj) 341 zprdch(ji,jj,jk) = zprdch(ji,jj,jk) * zmixdiat(ji,jj) 342 ENDIF 343 IF( etot_ndcy(ji,jj,jk) > 1.E-3 ) THEN 344 ! production terms for nanophyto. ( chlorophyll ) 345 znanotot = enano(ji,jj,jk) * zstrn(ji,jj) 346 zprod = rday * zprorca(ji,jj,jk) * zprnch(ji,jj,jk) * xlimphy(ji,jj,jk) 347 zprochln(ji,jj,jk) = chlcmin * 12. * zprorca (ji,jj,jk) 348 zprochln(ji,jj,jk) = zprochln(ji,jj,jk) + (chlcnm-chlcmin) * 12. * zprod / & 307 DO jk = 1, jpkm1 308 DO jj = 1, jpj 309 DO ji = 1, jpi 310 IF( etot_ndcy(ji,jj,jk) > 1.E-3 ) THEN 311 ! production terms for nanophyto. ( chlorophyll ) 312 znanotot = enano(ji,jj,jk) / ( zmxl_chl(ji,jj,jk) + rtrn ) 313 zprod = rday * zprorca(ji,jj,jk) * zprnch(ji,jj,jk) * xlimphy(ji,jj,jk) 314 zprochln(ji,jj,jk) = chlcmin * 12. * zprorca (ji,jj,jk) 315 zprochln(ji,jj,jk) = zprochln(ji,jj,jk) + (chlcnm-chlcmin) * 12. * zprod / & 349 316 & ( zpislopead(ji,jj,jk) * znanotot +rtrn) 350 ! production terms for diatomees ( chlorophyll ) 351 zdiattot = ediat(ji,jj,jk) * zstrn(ji,jj) 352 zprod = rday * zprorcad(ji,jj,jk) * zprdch(ji,jj,jk) * xlimdia(ji,jj,jk) 353 zprochld(ji,jj,jk) = chlcmin * 12. * zprorcad(ji,jj,jk) 354 zprochld(ji,jj,jk) = zprochld(ji,jj,jk) + (chlcdm-chlcmin) * 12. * zprod / & 355 & ( zpislopead2(ji,jj,jk) * zdiattot +rtrn ) 356 ENDIF 357 END DO 358 END DO 359 END DO 360 ELSE 361 !CDIR NOVERRCHK 362 DO jk = 1, jpkm1 363 !CDIR NOVERRCHK 364 DO jj = 1, jpj 365 !CDIR NOVERRCHK 366 DO ji = 1, jpi 367 IF( etot_ndcy(ji,jj,jk) > 1.E-3 ) THEN 368 ! production terms for nanophyto. ( chlorophyll ) 369 znanotot = enano(ji,jj,jk) 370 zprod = rday * zprorca(ji,jj,jk) * zprnch(ji,jj,jk) * trb(ji,jj,jk,jpphy) * xlimphy(ji,jj,jk) 371 zprochln(ji,jj,jk) = chlcmin * 12. * zprorca (ji,jj,jk) 372 zprochln(ji,jj,jk) = zprochln(ji,jj,jk) + (chlcnm-chlcmin) * 144. * zprod & 373 & / ( zpislopead(ji,jj,jk) * trb(ji,jj,jk,jpnch) * znanotot +rtrn ) 374 ! production terms for diatomees ( chlorophyll ) 375 zdiattot = ediat(ji,jj,jk) 376 zprod = rday * zprorcad(ji,jj,jk) * zprdch(ji,jj,jk) * trb(ji,jj,jk,jpdia) * xlimdia(ji,jj,jk) 377 zprochld(ji,jj,jk) = chlcmin * 12. * zprorcad(ji,jj,jk) 378 zprochld(ji,jj,jk) = zprochld(ji,jj,jk) + (chlcdm-chlcmin) * 144. * zprod & 379 & / ( zpislopead2(ji,jj,jk) * trb(ji,jj,jk,jpdch) * zdiattot +rtrn ) 380 ENDIF 381 END DO 382 END DO 383 END DO 384 ENDIF 317 ! production terms for diatomees ( chlorophyll ) 318 zdiattot = ediat(ji,jj,jk) / ( zmxl_chl(ji,jj,jk) + rtrn ) 319 zprod = rday * zprorcad(ji,jj,jk) * zprdch(ji,jj,jk) * xlimdia(ji,jj,jk) 320 zprochld(ji,jj,jk) = chlcmin * 12. * zprorcad(ji,jj,jk) 321 zprochld(ji,jj,jk) = zprochld(ji,jj,jk) + (chlcdm-chlcmin) * 12. * zprod / & 322 & ( zpislopead2(ji,jj,jk) * zdiattot +rtrn ) 323 ENDIF 324 END DO 325 END DO 326 END DO 385 327 386 328 ! Update the arrays TRA which contain the biological sources and sinks … … 550 492 ENDIF 551 493 ! 552 CALL wrk_dealloc( jpi, jpj, zmixnano, zmixdiat, zstrn ) 553 CALL wrk_dealloc( jpi, jpj, jpk, zpislopead, zpislopead2, zprdia, zprbio, zprdch, zprnch, zysopt ) 494 CALL wrk_dealloc( jpi, jpj, zstrn ) 495 CALL wrk_dealloc( jpi, jpj, jpk, zpislopead, zpislopead2, zprdia, zprbio, zprdch, zprnch, zysopt ) 496 CALL wrk_dealloc( jpi, jpj, jpk, zmxl_fac, zmxl_chl ) 554 497 CALL wrk_dealloc( jpi, jpj, jpk, zprorca, zprorcad, zprofed, zprofen, zprochln, zprochld, zpronew, zpronewd ) 555 498 ! -
branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zrem.F90
r6453 r6841 163 163 ! below 2 umol/L. Inhibited at strong light 164 164 ! ---------------------------------------------------------- 165 zonitr =nitrif * zstep * trb(ji,jj,jk,jpnh4) / ( 1.+ emoy(ji,jj,jk) ) * ( 1.- nitrfac(ji,jj,jk) ) 165 zonitr = nitrif * zstep * trb(ji,jj,jk,jpnh4) * ( 1.- nitrfac(ji,jj,jk) ) & 166 & / ( 1.+ emoy(ji,jj,jk) ) * ( 1. + fr_i(ji,jj) * emoy(ji,jj,jk) ) 166 167 denitnh4(ji,jj,jk) = nitrif * zstep * trb(ji,jj,jk,jpnh4) * nitrfac(ji,jj,jk) 167 168 ! Update of the tracers trends … … 229 230 ! constant and specified in the namelist. 230 231 ! ---------------------------------------------------------- 231 zdep = MAX( hmld(ji,jj), heup (ji,jj) )232 zdep = MAX( hmld(ji,jj), heup_01(ji,jj) ) 232 233 zdep = MAX( 0., fsdept(ji,jj,jk) - zdep ) 233 234 ztem = MAX( tsn(ji,jj,1,jp_tem), 0. ) -
branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zsed.F90
r6455 r6841 393 393 zlight = ( 1.- EXP( -etot_ndcy(ji,jj,jk) / diazolight ) ) 394 394 nitrpot(ji,jj,jk) = MAX( 0.e0, ( 0.6 * tgfunc(ji,jj,jk) - 2.15 ) * r1_rday ) & 395 & * zfact * MIN( ztrfer, ztrpo4 ) * zlight 395 & * zfact * MIN( ztrfer, ztrpo4 ) * zlight * (1. - fr_i(ji,jj)) 396 396 zsoufer(ji,jj,jk) = zlight * 2E-11 / (2E-11 + biron(ji,jj,jk)) 397 397 END DO -
branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zsink.F90
r6453 r6841 95 95 INTEGER :: ji, jj, jk, jit 96 96 INTEGER :: iiter1, iiter2 97 REAL(wp) :: zfact, zwsmax, zmax , zstep97 REAL(wp) :: zfact, zwsmax, zmax 98 98 CHARACTER (len=25) :: charout 99 99 REAL(wp), POINTER, DIMENSION(:,:,:) :: zw3d … … 103 103 IF( nn_timing == 1 ) CALL timing_start('p4z_sink') 104 104 105 ! Initialization of some global variables 106 ! --------------------------------------- 107 prodpoc(:,:,:) = 0. 108 conspoc(:,:,:) = 0. 109 prodgoc(:,:,:) = 0. 110 consgoc(:,:,:) = 0. 105 111 ! 106 112 ! Sinking speeds of detritus is increased with depth as shown … … 110 116 DO jj = 1, jpj 111 117 DO ji = 1,jpi 112 zmax = MAX( heup (ji,jj), hmld(ji,jj) )118 zmax = MAX( heup_01(ji,jj), hmld(ji,jj) ) 113 119 zfact = MAX( 0., fsdepw(ji,jj,jk+1) - zmax ) / wsbio2scale 114 120 wsbio4(ji,jj,jk) = wsbio2 + MAX(0., ( wsbio2max - wsbio2 )) * zfact … … 119 125 ! limit the values of the sinking speeds to avoid numerical instabilities 120 126 wsbio3(:,:,:) = wsbio 121 wscal (:,:,:)= wsbio4(:,:,:)127 wscal(:,:,:) = wsbio4(:,:,:) 122 128 #if defined key_ligand 123 129 wsfep (:,:,:) = wfep … … 198 204 CALL p4z_sink2( wsbio4, sinking2, jpgoc, iiter2 ) 199 205 CALL p4z_sink2( wsbio4, sinkfer2, jpbfe, iiter2 ) 200 CALL p4z_sink2( ws bio4, sinksil , jpgsi, iiter2 )201 CALL p4z_sink2( wscal 206 CALL p4z_sink2( wscal, sinksil , jpgsi, iiter2 ) 207 CALL p4z_sink2( wscal, sinkcal , jpcal, iiter2 ) 202 208 END DO 203 209 -
branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zsms.F90
r6453 r6841 455 455 REAL(wp) :: zrdenittot, zsdenittot, znitrpottot 456 456 CHARACTER(LEN=100) :: cltxt 457 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zvol458 457 INTEGER :: jk 459 458 !!---------------------------------------------------------------------- -
branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P5Z/p5zagg.F90
r6453 r6841 60 60 ! 61 61 IF( nn_timing == 1 ) CALL timing_start('p5z_agg') 62 63 ! Initialization of some global variables 64 ! --------------------------------------- 65 prodpoc(:,:,:) = 0. 66 conspoc(:,:,:) = 0. 67 prodgoc(:,:,:) = 0. 68 consgoc(:,:,:) = 0. 69 62 ! 70 63 ! Exchange between organic matter compartments due to coagulation/disaggregation 71 64 ! --------------------------------------------------- … … 117 110 118 111 ! tranfer of DOC to POC due to brownian motion 119 zaggtmp = ( 5095. * trb(ji,jj,jk,jppoc) +114. * 0.3 * trb(ji,jj,jk,jpdoc) ) *zstep112 zaggtmp = ( 114. * 0.3 * trb(ji,jj,jk,jpdoc) ) *zstep 120 113 zaggdoc3 = zaggtmp * 0.3 * trb(ji,jj,jk,jpdoc) 121 114 zaggdon3 = zaggtmp * 0.3 * trb(ji,jj,jk,jpdon) … … 135 128 tra(ji,jj,jk,jpdop) = tra(ji,jj,jk,jpdop) - zaggdop - zaggdop2 - zaggdop3 136 129 ! 137 prodpoc(ji,jj,jk) = prodpoc(ji,jj,jk) + zaggdoc + zaggdoc3 138 conspoc(ji,jj,jk) = conspoc(ji,jj,jk) - zaggpoc 130 conspoc(ji,jj,jk) = conspoc(ji,jj,jk) - zagg + zaggdoc + zaggdoc3 139 131 prodgoc(ji,jj,jk) = prodgoc(ji,jj,jk) + zaggpoc + zaggdoc2 140 132 ! -
branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P5Z/p5zpoc.F90
r6453 r6841 32 32 PUBLIC p5z_poc ! called in p4zbio.F90 33 33 PUBLIC p5z_poc_init ! called in trcsms_pisces.F90 34 PUBLIC alngam 35 PUBLIC gamain 34 36 35 37 !! * Shared module variables … … 38 40 REAL(wp), PUBLIC :: xremipp !: remineralisation rate of DOP 39 41 INTEGER , PUBLIC :: jcpoc !: number of lability classes 42 REAL(wp), PUBLIC :: rshape !: shape factor of the gamma distribution 40 43 41 44 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: alphan, reminp … … 63 66 ! 64 67 INTEGER :: ji, jj, jk, jn 65 REAL(wp) :: zremip, zremig, zdep, z tem, zstep68 REAL(wp) :: zremip, zremig, zdep, zstep 66 69 REAL(wp) :: zopon, zopop, zopoc2, zopon2, zopop2, zofer 67 70 REAL(wp) :: zsizek, zsizek1, alphat, remint 68 REAL(wp) :: solgoc, zpoc , zpoctot, zremif71 REAL(wp) :: solgoc, zpoc 69 72 #if ! defined key_kriest 70 73 REAL(wp) :: zofer2, zofer3 … … 102 105 DO jn = 1, jcpoc 103 106 alphag(:,:,:,jn) = alphan(jn) 107 alphap(:,:,:,jn) = alphan(jn) 104 108 END DO 105 109 106 110 #if ! defined key_kriest 111 ! ----------------------------------------------------------------------- 112 ! Lability parameterization. This is the big particles part (GOC) 113 ! This lability parameterization can be activated only with the standard 114 ! particle scheme. Does not work with Kriest parameterization. 115 ! ----------------------------------------------------------------------- 107 116 DO jk = 2, jpkm1 108 117 DO jj = 1, jpj … … 110 119 IF (tmask(ji,jj,jk) == 1.) THEN 111 120 zdep = hmld(ji,jj) 112 IF( fsdept(ji,jj,jk) >= zdep ) THEN 121 ! 122 ! In the case of GOC, lability is constant in the mixed layer 123 ! It is computed only below the mixed layer depth 124 ! ------------------------------------------------------------ 125 ! 126 IF( fsdept(ji,jj,jk) > zdep ) THEN 113 127 alphat = 0. 114 128 remint = 0. 115 IF ( fsdept(ji,jj,jk-1) < zdep ) THEN 129 ! 130 IF ( fsdept(ji,jj,jk-1) <= zdep ) THEN 131 ! 132 ! The first level just below the mixed layer needs a 133 ! specific treatment because lability is supposed constant 134 ! everywhere within the mixed layer. This means that 135 ! change in lability in the bottom part of the previous cell 136 ! should not be computed 137 ! ---------------------------------------------------------- 138 ! 139 ! POC concentration is computed using the lagrangian 140 ! framework. It is only used for the lability param 141 zpoc = trb(ji,jj,jk-1,jpgoc) + consgoc(ji,jj,jk) * rday / rfact2 & 142 & * fse3t(ji,jj,jk) / 2. / (wsbio4(ji,jj,jk) + rtrn) 143 zpoc = max(0., zpoc) 144 ! 116 145 DO jn = 1, jcpoc 146 ! 147 ! Lagrangian based algorithm. The fraction of each 148 ! lability class is computed starting from the previous 149 ! level 150 ! ----------------------------------------------------- 151 ! 117 152 zsizek = zdep / (wsbio2 + rtrn) * tgfunc(ji,jj,jk-1) 118 153 zsizek1 = fse3t(ji,jj,jk) / 2. / (wsbio4(ji,jj,jk) + rtrn) * tgfunc(ji,jj,jk) 119 zpoc = trb(ji,jj,jk-1,jpgoc) + consgoc(ji,jj,jk) * rday / rfact2 & 120 & * fse3t(ji,jj,jk) / 2. / (wsbio4(ji,jj,jk) + rtrn) 121 zpoc = max(0., zpoc) 154 ! the concentration of each lability class is calculated 155 ! as the sum of the different sources and sinks 156 ! Please note that production of new GOC experiences 157 ! degradation 122 158 alphag(ji,jj,jk,jn) = alphan(jn) / (reminp(jn) * tgfunc(ji,jj,jk-1) ) & 123 159 & * (1. - exp( -reminp(jn) * zsizek ) ) * exp( -reminp(jn) * zsizek1 ) & … … 128 164 END DO 129 165 ELSE 166 ! 167 ! standard algorithm in the rest of the water column 168 ! See the comments in the previous block. 169 ! --------------------------------------------------- 170 ! 171 zpoc = trb(ji,jj,jk-1,jpgoc) + consgoc(ji,jj,jk-1) * rday / rfact2 & 172 & * fse3t(ji,jj,jk-1) / 2. / (wsbio4(ji,jj,jk-1) + rtrn) + consgoc(ji,jj,jk) & 173 & * rday / rfact2 * fse3t(ji,jj,jk) / 2. / (wsbio4(ji,jj,jk) + rtrn) 174 zpoc = max(0., zpoc) 175 ! 130 176 DO jn = 1, jcpoc 131 177 zsizek = fse3t(ji,jj,jk-1) / 2. / (wsbio4(ji,jj,jk-1) + rtrn) * tgfunc(ji,jj,jk-1) 132 178 zsizek1 = fse3t(ji,jj,jk) / 2. / (wsbio4(ji,jj,jk) + rtrn) * tgfunc(ji,jj,jk) 133 zpoc = trb(ji,jj,jk-1,jpgoc) + consgoc(ji,jj,jk-1) * rday / rfact2 &134 & * fse3t(ji,jj,jk-1) / 2. / (wsbio4(ji,jj,jk-1) + rtrn) + consgoc(ji,jj,jk) &135 & * rday / rfact2 * fse3t(ji,jj,jk) / 2. / (wsbio4(ji,jj,jk) + rtrn)136 zpoc = max(0., zpoc)137 179 alphag(ji,jj,jk,jn) = alphag(ji,jj,jk-1,jn) * exp( -reminp(jn) & 138 180 & * ( zsizek + zsizek1 ) ) * zpoc + ( prodgoc(ji,jj,jk-1) & … … 146 188 ENDIF 147 189 DO jn = 1, jcpoc 190 ! The contribution of each lability class at the current 191 ! level is computed 148 192 alphag(ji,jj,jk,jn) = alphag(ji,jj,jk,jn) / ( alphat + rtrn) 149 193 END DO 194 ! Computation of the mean remineralisation rate 150 195 zremigoc(ji,jj,jk) = MIN(xremipc, MAX(0., remint / ( alphat + rtrn) )) 151 196 ENDIF … … 195 240 ENDIF 196 241 197 totprod(:,:) = 0. 198 totthick(:,:) = 0. 199 totcons(:,:) = 0. 200 DO jk = 1, jpkm1 201 DO jj = 1, jpj 202 DO ji = 1, jpi 203 IF (tmask(ji,jj,jk) == 1.) THEN 204 zdep = hmld(ji,jj) 205 IF( fsdept(ji,jj,jk) <= zdep ) THEN 206 totprod(ji,jj) = totprod(ji,jj) + prodpoc(ji,jj,jk) * fse3t(ji,jj,jk) * rday/ rfact2 207 totthick(ji,jj) = totthick(ji,jj) + fse3t(ji,jj,jk)* tgfunc(ji,jj,jk) 208 totcons(ji,jj) = totcons(ji,jj) - conspoc(ji,jj,jk) * fse3t(ji,jj,jk) * rday/ rfact2 & 242 ! ------------------------------------------------------------------ 243 ! Lability parameterization for the small OM particles. This param 244 ! is based on the same theoretical background as the big particles. 245 ! However, because of its low sinking speed, lability is not supposed 246 ! to be equal to its initial value (the value of the freshly produced 247 ! organic matter). It is however uniform in the mixed layer. 248 ! ------------------------------------------------------------------- 249 ! 250 totprod(:,:) = 0. 251 totthick(:,:) = 0. 252 totcons(:,:) = 0. 253 ! intregrated production and consumption of POC in the mixed layer 254 ! ---------------------------------------------------------------- 255 ! 256 DO jk = 1, jpkm1 257 DO jj = 1, jpj 258 DO ji = 1, jpi 259 IF (tmask(ji,jj,jk) == 1.) THEN 260 zdep = hmld(ji,jj) 261 IF( fsdept(ji,jj,jk) <= zdep ) THEN 262 totprod(ji,jj) = totprod(ji,jj) + prodpoc(ji,jj,jk) * fse3t(ji,jj,jk) * rday/ rfact2 263 ! The temperature effect is included here 264 totthick(ji,jj) = totthick(ji,jj) + fse3t(ji,jj,jk)* tgfunc(ji,jj,jk) 265 totcons(ji,jj) = totcons(ji,jj) - conspoc(ji,jj,jk) * fse3t(ji,jj,jk) * rday/ rfact2 & 209 266 & / ( trb(ji,jj,jk,jppoc) + rtrn ) 210 ENDIF 211 ENDIF 212 END DO 213 END DO 214 END DO 215 267 ENDIF 268 ENDIF 269 END DO 270 END DO 271 END DO 272 273 ! Computation of the lability spectrum in the mixed layer. In the mixed 274 ! layer, this spectrum is supposed to be uniform. 275 ! --------------------------------------------------------------------- 216 276 DO jk = 1, jpkm1 217 277 DO jj = 1, jpj … … 223 283 IF( fsdept(ji,jj,jk) <= zdep ) THEN 224 284 DO jn = 1, jcpoc 285 ! For each lability class, the system is supposed to be 286 ! at equilibrium: Prod - Sink - w alphap = 0. 225 287 alphap(ji,jj,jk,jn) = totprod(ji,jj) * alphan(jn) / ( reminp(jn) & 226 288 & * totthick(ji,jj) + totcons(ji,jj) + wsbio + rtrn ) 227 289 alphat = alphat + alphap(ji,jj,jk,jn) 228 remint = remint + alphap(ji,jj,jk,jn) * reminp(jn)229 290 END DO 230 291 DO jn = 1, jcpoc 231 292 alphap(ji,jj,jk,jn) = alphap(ji,jj,jk,jn) / ( alphat + rtrn) 293 remint = remint + alphap(ji,jj,jk,jn) * reminp(jn) 232 294 END DO 233 zremipoc(ji,jj,jk) = MIN(xremipc, MAX(0., remint / ( alphat + rtrn) )) 295 ! Mean remineralization rate in the mixed layer 296 zremipoc(ji,jj,jk) = MIN(xremip, MAX(0., remint )) 234 297 ENDIF 235 298 ENDIF … … 237 300 END DO 238 301 END DO 239 302 ! 303 ! ----------------------------------------------------------------------- 304 ! The lability parameterization is used here. The code is here 305 ! almost identical to what is done for big particles. The only difference 306 ! is that an additional source from GOC to POC is included. This means 307 ! that since we need the lability spectrum of GOC, GOC spectrum 308 ! should be determined before. 309 ! ----------------------------------------------------------------------- 310 ! 240 311 DO jk = 2, jpkm1 241 312 DO jj = 1, jpj … … 246 317 alphat = 0. 247 318 remint = 0. 319 ! 320 ! Special treatment of the level just below the MXL 321 ! See the comments in the GOC section 322 ! --------------------------------------------------- 323 ! 248 324 IF ( fsdept(ji,jj,jk-1) <= zdep ) THEN 325 ! 326 ! Computation of the POC concentration using the 327 ! lagrangian algorithm 328 zpoc = trb(ji,jj,jk-1,jppoc) + conspoc(ji,jj,jk) * rday / rfact2 & 329 & * fse3t(ji,jj,jk) / 2. / (wsbio3(ji,jj,jk) + rtrn) 330 zpoc = max(0., zpoc) 331 ! 249 332 DO jn = 1, jcpoc 333 ! the scale factor is corrected with temperature 250 334 zsizek1 = fse3t(ji,jj,jk) / 2. / (wsbio3(ji,jj,jk) + rtrn) * tgfunc(ji,jj,jk) 251 zpoc = trb(ji,jj,jk-1,jppoc) + conspoc(ji,jj,jk) * rday / rfact2 & 252 & * fse3t(ji,jj,jk) / 2. / (wsbio3(ji,jj,jk) + rtrn) 253 zpoc = max(0., zpoc) 335 ! computation of the lability spectrum applying the 336 ! different sources and sinks 254 337 alphap(ji,jj,jk,jn) = alphap(ji,jj,jk-1,jn) + zorem3(ji,jj,jk) & 255 338 & * ( 1. - exp( -reminp(jn) * zsizek1 ) ) * rday / rfact2 & 256 339 & * alphag(ji,jj,jk,jn) / reminp(jn) / tgfunc(ji,jj,jk) 340 alphap(ji,jj,jk,jn) = MAX( 0., alphap(ji,jj,jk,jn) ) 257 341 alphat = alphat + alphap(ji,jj,jk,jn) 258 remint = remint + alphap(ji,jj,jk,jn) * reminp(jn)259 342 END DO 260 343 ELSE 344 ! 345 ! Lability parameterization for the interior of the ocean 346 ! This is very similar to what is done in the previous 347 ! block 348 ! -------------------------------------------------------- 349 ! 350 zpoc = trb(ji,jj,jk-1,jppoc) + conspoc(ji,jj,jk-1) * rday / rfact2 & 351 & * fse3t(ji,jj,jk-1) / 2. / (wsbio3(ji,jj,jk-1) + rtrn) + conspoc(ji,jj,jk) & 352 & * rday / rfact2 * fse3t(ji,jj,jk) / 2. / (wsbio3(ji,jj,jk) + rtrn) 353 zpoc = max(0., zpoc) 354 ! 261 355 DO jn = 1, jcpoc 262 356 zsizek = fse3t(ji,jj,jk-1) / 2. / (wsbio3(ji,jj,jk-1) + rtrn) * tgfunc(ji,jj,jk-1) 263 357 zsizek1 = fse3t(ji,jj,jk) / 2. / (wsbio3(ji,jj,jk) + rtrn) * tgfunc(ji,jj,jk) 264 zpoc = trb(ji,jj,jk-1,jppoc) + conspoc(ji,jj,jk-1) * rday / rfact2 &265 & * fse3t(ji,jj,jk-1) / 2. / (wsbio3(ji,jj,jk-1) + rtrn) + conspoc(ji,jj,jk) &266 & * rday / rfact2 * fse3t(ji,jj,jk) / 2. / (wsbio3(ji,jj,jk) + rtrn)267 zpoc = max(0., zpoc)268 358 alphap(ji,jj,jk,jn) = alphap(ji,jj,jk-1,jn) * exp( -reminp(jn) & 269 359 & * ( zsizek + zsizek1 ) ) * zpoc + ( prodpoc(ji,jj,jk-1) & … … 277 367 & + zorem3(ji,jj,jk) * ( 1. - exp( -reminp(jn) * zsizek1 ) ) * rday & 278 368 & / rfact2 * alphag(ji,jj,jk,jn) / reminp(jn) / tgfunc(ji,jj,jk) 369 alphap(ji,jj,jk,jn) = max(0., alphap(ji,jj,jk,jn) ) 279 370 alphat = alphat + alphap(ji,jj,jk,jn) 280 remint = remint + alphap(ji,jj,jk,jn) * reminp(jn)281 371 END DO 282 372 ENDIF 373 ! Normalization of the lability spectrum so that the 374 ! integral is equal to 1 283 375 DO jn = 1, jcpoc 284 376 alphap(ji,jj,jk,jn) = alphap(ji,jj,jk,jn) / ( alphat + rtrn) 377 remint = remint + alphap(ji,jj,jk,jn) * reminp(jn) 285 378 END DO 286 zremipoc(ji,jj,jk) = MIN(xremipc, MAX(0., remint / ( alphat + rtrn) )) 379 ! Mean remineralization rate in the water column 380 zremipoc(ji,jj,jk) = MIN(xremip, MAX(0., remint )) 287 381 ENDIF 288 382 ENDIF … … 364 458 INTEGER :: jn 365 459 REAL(wp) :: remindelta, reminup, remindown 366 NAMELIST/nampispoc/ xremipc, xremipn, xremipp, jcpoc 460 INTEGER :: ifault 461 462 NAMELIST/nampispoc/ xremipc, xremipn, xremipp, jcpoc, rshape 367 463 INTEGER :: ios ! Local integer output status for namelist read 368 464 … … 384 480 WRITE(numout,*) ' remineralisation rate of POP xremipp =', xremipp 385 481 WRITE(numout,*) ' Number of lability classes for POC jcpoc =', jcpoc 482 WRITE(numout,*) ' Shape factor of the gamma distribution rshape =', rshape 386 483 ENDIF 387 484 ! … … 394 491 ! 395 492 IF (jcpoc > 1) THEN 493 ! 396 494 remindelta = log(4. * 1000. ) / float(jcpoc-1) 397 reminup = xremipc/400. * exp(remindelta) 398 alphan(1) = 1. - exp(-reminup/xremipc) 399 reminp(1) = xremipc - reminup * exp(-reminup/xremipc) / alphan(1) 495 reminup = 1./ 400. * exp(remindelta) 496 ! 497 ! Discretization based on incomplete gamma functions 498 ! As incomplete gamma functions are not available in standard 499 ! fortran 95, they have been coded as functions in this module (gamain) 500 ! --------------------------------------------------------------------- 501 ! 502 alphan(1) = gamain(reminup, rshape, ifault) 503 reminp(1) = gamain(reminup, rshape+1.0, ifault) * xremip / alphan(1) 400 504 DO jn = 2, jcpoc-1 401 reminup = xremipc/400. * exp(float(jn) * remindelta)402 remindown = xremipc/ 400. * exp(float(jn-1) * remindelta)403 alphan(jn) = exp(-remindown /xremipc) - exp(-reminup/xremipc)404 reminp(jn) = xremipc + (remindown * exp(-remindown /xremipc) &405 & - reminup * exp(-reminup/xremipc) )/ alphan(jn)505 reminup = 1./ 400. * exp(float(jn) * remindelta) 506 remindown = 1. / 400. * exp(float(jn-1) * remindelta) 507 alphan(jn) = gamain(reminup, rshape, ifault) - gamain(remindown, rshape, ifault) 508 reminp(jn) = gamain(reminup, rshape+1.0, ifault) - gamain(remindown, rshape+1.0, ifault) 509 reminp(jn) = reminp(jn) * xremip / alphan(jn) 406 510 END DO 407 remindown = xremipc/400. * exp(float(jcpoc-1) * remindelta)408 alphan(jcpoc) = exp(-remindown /xremipc)409 reminp(jcpoc) = xremipc + remindown511 remindown = 1. / 400. * exp(float(jcpoc-1) * remindelta) 512 alphan(jcpoc) = 1.0 - gamain(remindown, rshape, ifault) 513 reminp(jcpoc) = reminp(jcpoc) * xremip / alphan(jcpoc) 410 514 ELSE 411 515 alphan(jcpoc) = 1. 412 reminp(jcpoc) = xremip c516 reminp(jcpoc) = xremip 413 517 ENDIF 414 518 … … 418 522 419 523 END SUBROUTINE p5z_poc_init 524 525 REAL FUNCTION alngam( xvalue, ifault ) 526 527 !*****************************************************************************80 528 ! 529 !! ALNGAM computes the logarithm of the gamma function. 530 ! 531 ! Modified: 532 ! 533 ! 13 January 2008 534 ! 535 ! Author: 536 ! 537 ! Allan Macleod 538 ! FORTRAN90 version by John Burkardt 539 ! 540 ! Reference: 541 ! 542 ! Allan Macleod, 543 ! Algorithm AS 245, 544 ! A Robust and Reliable Algorithm for the Logarithm of the Gamma Function, 545 ! Applied Statistics, 546 ! Volume 38, Number 2, 1989, pages 397-402. 547 ! 548 ! Parameters: 549 ! 550 ! Input, real ( kind = 8 ) XVALUE, the argument of the Gamma function. 551 ! 552 ! Output, integer ( kind = 4 ) IFAULT, error flag. 553 ! 0, no error occurred. 554 ! 1, XVALUE is less than or equal to 0. 555 ! 2, XVALUE is too big. 556 ! 557 ! Output, real ( kind = 8 ) ALNGAM, the logarithm of the gamma function of X. 558 ! 559 implicit none 560 561 real(wp), parameter :: alr2pi = 0.918938533204673E+00 562 integer:: ifault 563 real(wp), dimension ( 9 ) :: r1 = (/ & 564 -2.66685511495E+00, & 565 -24.4387534237E+00, & 566 -21.9698958928E+00, & 567 11.1667541262E+00, & 568 3.13060547623E+00, & 569 0.607771387771E+00, & 570 11.9400905721E+00, & 571 31.4690115749E+00, & 572 15.2346874070E+00 /) 573 real(wp), dimension ( 9 ) :: r2 = (/ & 574 -78.3359299449E+00, & 575 -142.046296688E+00, & 576 137.519416416E+00, & 577 78.6994924154E+00, & 578 4.16438922228E+00, & 579 47.0668766060E+00, & 580 313.399215894E+00, & 581 263.505074721E+00, & 582 43.3400022514E+00 /) 583 real(wp), dimension ( 9 ) :: r3 = (/ & 584 -2.12159572323E+05, & 585 2.30661510616E+05, & 586 2.74647644705E+04, & 587 -4.02621119975E+04, & 588 -2.29660729780E+03, & 589 -1.16328495004E+05, & 590 -1.46025937511E+05, & 591 -2.42357409629E+04, & 592 -5.70691009324E+02 /) 593 real(wp), dimension ( 5 ) :: r4 = (/ & 594 0.279195317918525E+00, & 595 0.4917317610505968E+00, & 596 0.0692910599291889E+00, & 597 3.350343815022304E+00, & 598 6.012459259764103E+00 /) 599 real (wp) :: x 600 real (wp) :: x1 601 real (wp) :: x2 602 real (wp), parameter :: xlge = 5.10E+05 603 real (wp), parameter :: xlgst = 1.0E+30 604 real (wp) :: xvalue 605 real (wp) :: y 606 607 x = xvalue 608 alngam = 0.0E+00 609 ! 610 ! Check the input. 611 ! 612 if ( xlgst <= x ) then 613 ifault = 2 614 return 615 end if 616 if ( x <= 0.0E+00 ) then 617 ifault = 1 618 return 619 end if 620 ifault = 0 621 ! 622 ! Calculation for 0 < X < 0.5 and 0.5 <= X < 1.5 combined. 623 ! 624 if ( x < 1.5E+00 ) then 625 626 if ( x < 0.5E+00 ) then 627 alngam = - log ( x ) 628 y = x + 1.0E+00 629 ! 630 ! Test whether X < machine epsilon. 631 ! 632 if ( y == 1.0E+00 ) then 633 return 634 end if 635 636 else 637 638 alngam = 0.0E+00 639 y = x 640 x = ( x - 0.5E+00 ) - 0.5E+00 641 642 end if 643 644 alngam = alngam + x * (((( & 645 r1(5) * y & 646 + r1(4) ) * y & 647 + r1(3) ) * y & 648 + r1(2) ) * y & 649 + r1(1) ) / (((( & 650 y & 651 + r1(9) ) * y & 652 + r1(8) ) * y & 653 + r1(7) ) * y & 654 + r1(6) ) 655 656 return 657 658 end if 659 ! 660 ! Calculation for 1.5 <= X < 4.0. 661 ! 662 if ( x < 4.0E+00 ) then 663 664 y = ( x - 1.0E+00 ) - 1.0E+00 665 666 alngam = y * (((( & 667 r2(5) * x & 668 + r2(4) ) * x & 669 + r2(3) ) * x & 670 + r2(2) ) * x & 671 + r2(1) ) / (((( & 672 x & 673 + r2(9) ) * x & 674 + r2(8) ) * x & 675 + r2(7) ) * x & 676 + r2(6) ) 677 ! 678 ! Calculation for 4.0 <= X < 12.0. 679 ! 680 else if ( x < 12.0E+00 ) then 681 682 alngam = (((( & 683 r3(5) * x & 684 + r3(4) ) * x & 685 + r3(3) ) * x & 686 + r3(2) ) * x & 687 + r3(1) ) / (((( & 688 x & 689 + r3(9) ) * x & 690 + r3(8) ) * x & 691 + r3(7) ) * x & 692 + r3(6) ) 693 ! 694 ! Calculation for 12.0 <= X. 695 ! 696 else 697 698 y = log ( x ) 699 alngam = x * ( y - 1.0E+00 ) - 0.5E+00 * y + alr2pi 700 701 if ( x <= xlge ) then 702 703 x1 = 1.0E+00 / x 704 x2 = x1 * x1 705 706 alngam = alngam + x1 * ( ( & 707 r4(3) * & 708 x2 + r4(2) ) * & 709 x2 + r4(1) ) / ( ( & 710 x2 + r4(5) ) * & 711 x2 + r4(4) ) 712 713 end if 714 end if 715 716 END FUNCTION alngam 717 718 REAL FUNCTION gamain( x, p, ifault ) 719 720 !*****************************************************************************80 721 ! 722 !! GAMAIN computes the incomplete gamma ratio. 723 ! 724 ! Discussion: 725 ! 726 ! A series expansion is used if P > X or X <= 1. Otherwise, a 727 ! continued fraction approximation is used. 728 ! 729 ! Modified: 730 ! 731 ! 17 January 2008 732 ! 733 ! Author: 734 ! 735 ! G Bhattacharjee 736 ! FORTRAN90 version by John Burkardt 737 ! 738 ! Reference: 739 ! 740 ! G Bhattacharjee, 741 ! Algorithm AS 32: 742 ! The Incomplete Gamma Integral, 743 ! Applied Statistics, 744 ! Volume 19, Number 3, 1970, pages 285-287. 745 ! 746 ! Parameters: 747 ! 748 ! Input, real ( kind = 8 ) X, P, the parameters of the incomplete 749 ! gamma ratio. 0 <= X, and 0 < P. 750 ! 751 ! Output, integer ( kind = 4 ) IFAULT, error flag. 752 ! 0, no errors. 753 ! 1, P <= 0. 754 ! 2, X < 0. 755 ! 3, underflow. 756 ! 4, error return from the Log Gamma routine. 757 ! 758 ! Output, real ( kind = 8 ) GAMAIN, the value of the incomplete 759 ! gamma ratio. 760 ! 761 implicit none 762 real (wp) a 763 real (wp), parameter :: acu = 1.0E-08 764 real (wp) an 765 real (wp) arg 766 real (wp) b 767 real (wp) dif 768 real (wp) factor 769 real (wp) g 770 real (wp) gin 771 integer i 772 integer ifault 773 real (wp), parameter :: oflo = 1.0E+37 774 real (wp) p 775 real (wp) pn(6) 776 real (wp) rn 777 real (wp) term 778 real (wp), parameter :: uflo = 1.0E-37 779 real (wp) x 780 ! 781 ! Check the input. 782 ! 783 if ( p <= 0.0E+00 ) then 784 ifault = 1 785 gamain = 0.0E+00 786 return 787 end if 788 789 if ( x < 0.0E+00 ) then 790 ifault = 2 791 gamain = 0.0E+00 792 return 793 end if 794 795 if ( x == 0.0E+00 ) then 796 ifault = 0 797 gamain = 0.0E+00 798 return 799 end if 800 801 g = alngam ( p, ifault ) 802 803 if ( ifault /= 0 ) then 804 ifault = 4 805 gamain = 0.0E+00 806 return 807 end if 808 809 arg = p * log ( x ) - x - g 810 if ( arg < log ( uflo ) ) then 811 ifault = 3 812 gamain = 0.0E+00 813 return 814 end if 815 816 ifault = 0 817 factor = exp ( arg ) 818 ! 819 ! Calculation by series expansion. 820 ! 821 if ( x <= 1.0E+00 .or. x < p ) then 822 823 gin = 1.0E+00 824 term = 1.0E+00 825 rn = p 826 827 do 828 829 rn = rn + 1.0E+00 830 term = term * x / rn 831 gin = gin + term 832 833 if ( term <= acu ) then 834 exit 835 end if 836 837 end do 838 839 gamain = gin * factor / p 840 return 841 842 end if 843 ! 844 ! Calculation by continued fraction. 845 ! 846 a = 1.0E+00 - p 847 b = a + x + 1.0E+00 848 term = 0.0E+00 849 850 pn(1) = 1.0E+00 851 pn(2) = x 852 pn(3) = x + 1.0E+00 853 pn(4) = x * b 854 855 gin = pn(3) / pn(4) 856 do 857 858 a = a + 1.0E+00 859 b = b + 2.0E+00 860 term = term + 1.0E+00 861 an = a * term 862 do i = 1, 2 863 pn(i+4) = b * pn(i+2) - an * pn(i) 864 end do 865 866 if ( pn(6) /= 0.0E+00 ) then 867 868 rn = pn(5) / pn(6) 869 dif = abs ( gin - rn ) 870 ! 871 ! Absolute error tolerance satisfied? 872 ! 873 if ( dif <= acu ) then 874 ! 875 ! Relative error tolerance satisfied? 876 ! 877 if ( dif <= acu * rn ) then 878 gamain = 1.0E+00 - factor * gin 879 exit 880 end if 881 882 end if 883 884 gin = rn 885 886 end if 887 888 do i = 1, 4 889 pn(i) = pn(i+2) 890 end do 891 if ( oflo <= abs ( pn(5) ) ) then 892 893 do i = 1, 4 894 pn(i) = pn(i) / oflo 895 end do 896 897 end if 898 899 end do 900 901 END FUNCTION gamain 420 902 421 903 #else -
branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P5Z/p5zprod.F90
r6453 r6841 161 161 IF( etot(ji,jj,jk) > 1.E-3 ) THEN 162 162 zval = MAX( 1., zstrn(ji,jj) ) 163 zval = 1.5 * zval / (12. + zval) 163 zval = 1.5 * zval / (12. + zval) * (1. - fr_i(ji,jj)) 164 164 zprbio(ji,jj,jk) = prmaxn(ji,jj,jk) * zval 165 165 zprpic(ji,jj,jk) = prmaxp(ji,jj,jk) * zval -
branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P5Z/p5zrem.F90
r6453 r6841 187 187 ! below 2 umol/L. Inhibited at strong light 188 188 ! ---------------------------------------------------------- 189 zonitr = nitrif * zstep * trb(ji,jj,jk,jpnh4) * (0.2 + 0.8 / ( 1.+ emoy(ji,jj,jk) ) )&190 & * ( 1.- nitrfac(ji,jj,jk) )189 zonitr = nitrif * zstep * trb(ji,jj,jk,jpnh4) * ( 1.- nitrfac(ji,jj,jk) ) & 190 & * ( 0.2 + 0.8 / ( 1.+ emoy(ji,jj,jk) ) * ( 1. + fr_i(ji,jj) * emoy(ji,jj,jk) ) ) 191 191 zdenitnh4 = nitrif * zstep * trb(ji,jj,jk,jpnh4) * nitrfac(ji,jj,jk) 192 192 ! 193 193 ! Update of the tracers trends 194 194 ! ---------------------------- … … 255 255 ! constant and specified in the namelist. 256 256 ! ---------------------------------------------------------- 257 zdep = MAX( hmld(ji,jj), heup (ji,jj) )257 zdep = MAX( hmld(ji,jj), heup_01(ji,jj) ) 258 258 zdep = MAX( 0., fsdept(ji,jj,jk) - zdep ) 259 259 ztem = MAX( tsn(ji,jj,1,jp_tem), 0. ) -
branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P5Z/p5zsed.F90
r6453 r6841 414 414 ztrdop(ji,jj,jk) = trb(ji,jj,jk,jpdop) / ( 1E-6 + trb(ji,jj,jk,jpdop) ) * (1. - ztrpo4(ji,jj,jk)) 415 415 ztrdp = ztrpo4(ji,jj,jk) + ztrdop(ji,jj,jk) 416 zlight = ( 1.- EXP( -etot(ji,jj,jk) / diazolight ) ) 416 zlight = ( 1.- EXP( -etot(ji,jj,jk) / diazolight ) ) * (1. - fr_i(ji,jj)) 417 417 nitrpot(ji,jj,jk) = zmudia * r1_rday * zfact * MIN( ztrfer, ztrdp ) * zlight 418 418 zsoufer(ji,jj,jk) = zlight * 2E-11 / (2E-11 + biron(ji,jj,jk)) -
branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P5Z/p5zsink.F90
r6453 r6841 98 98 INTEGER :: ji, jj, jk, jit 99 99 INTEGER :: iiter1, iiter2 100 REAL(wp) :: zfact, zwsmax, zmax , zstep100 REAL(wp) :: zfact, zwsmax, zmax 101 101 CHARACTER (len=25) :: charout 102 102 REAL(wp), POINTER, DIMENSION(:,:,:) :: zw3d … … 106 106 IF( nn_timing == 1 ) CALL timing_start('p5z_sink') 107 107 108 ! Initialization of some global variables 109 ! --------------------------------------- 110 prodpoc(:,:,:) = 0. 111 conspoc(:,:,:) = 0. 112 prodgoc(:,:,:) = 0. 113 consgoc(:,:,:) = 0. 108 114 ! 109 115 ! Sinking speeds of detritus is increased with depth as shown … … 113 119 DO jj = 1, jpj 114 120 DO ji = 1,jpi 115 zmax = MAX( heup (ji,jj), hmld(ji,jj) )121 zmax = MAX( heup_01(ji,jj), hmld(ji,jj) ) 116 122 zfact = MAX( 0., fsdepw(ji,jj,jk+1) - zmax ) / wsbio2scale 117 123 wsbio4(ji,jj,jk) = wsbio2 + MAX(0., ( wsbio2max - wsbio2 )) * zfact -
branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P5Z/p5zsms.F90
r6453 r6841 466 466 REAL(wp) :: zrdenittot, zsdenittot, znitrpottot 467 467 CHARACTER(LEN=100) :: cltxt 468 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zvol469 468 INTEGER :: jk 470 469 !!
Note: See TracChangeset
for help on using the changeset viewer.
