- Timestamp:
- 2016-08-03T16:59:29+02:00 (8 years ago)
- Location:
- branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P5Z
- Files:
-
- 7 edited
-
p5zagg.F90 (modified) (3 diffs)
-
p5zpoc.F90 (modified) (16 diffs)
-
p5zprod.F90 (modified) (1 diff)
-
p5zrem.F90 (modified) (2 diffs)
-
p5zsed.F90 (modified) (1 diff)
-
p5zsink.F90 (modified) (3 diffs)
-
p5zsms.F90 (modified) (1 diff)
Legend:
- Unmodified
- Added
- Removed
-
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.
