- Timestamp:
- 2016-08-03T16:59:29+02:00 (8 years ago)
- Location:
- branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z
- Files:
-
- 15 edited
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 !!----------------------------------------------------------------------
Note: See TracChangeset
for help on using the changeset viewer.