- Timestamp:
- 2016-11-30T17:56:53+01:00 (8 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2016/dev_merge_2016/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zsink.F90
r6140 r7403 9 9 !! 3.5 ! 2012-07 (O. Aumont) Introduce potential time-splitting 10 10 !!---------------------------------------------------------------------- 11 #if defined key_pisces12 !!----------------------------------------------------------------------13 11 !! p4z_sink : Compute vertical flux of particulate matter due to gravitational sinking 14 12 !! p4z_sink_init : Unitialisation of sinking speed parameters … … 29 27 PUBLIC p4z_sink_alloc 30 28 31 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: wsbio3 !: POC sinking speed32 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: wsbio4 !: GOC sinking speed33 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: wscal !: Calcite and BSi sinking speeds34 35 29 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: sinking, sinking2 !: POC sinking fluxes 36 30 ! ! (different meanings depending on the parameterization) 31 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: sinkingn, sinking2n !: POC sinking fluxes 32 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: sinkingp, sinking2p !: POC sinking fluxes 37 33 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: sinkcal, sinksil !: CaCO3 and BSi sinking fluxes 38 34 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: sinkfer !: Small BFe sinking fluxes 39 #if ! defined key_kriest40 35 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: sinkfer2 !: Big iron sinking fluxes 41 #endif 36 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: sinkfep !: Fep sinking fluxes 42 37 43 38 INTEGER :: ik100 44 45 #if defined key_kriest46 REAL(wp) :: xkr_sfact !: Sinking factor47 REAL(wp) :: xkr_stick !: Stickiness48 REAL(wp) :: xkr_nnano !: Nbr of cell in nano size class49 REAL(wp) :: xkr_ndiat !: Nbr of cell in diatoms size class50 REAL(wp) :: xkr_nmicro !: Nbr of cell in microzoo size class51 REAL(wp) :: xkr_nmeso !: Nbr of cell in mesozoo size class52 REAL(wp) :: xkr_naggr !: Nbr of cell in aggregates size class53 54 REAL(wp) :: xkr_frac55 56 REAL(wp), PUBLIC :: xkr_dnano !: Size of particles in nano pool57 REAL(wp), PUBLIC :: xkr_ddiat !: Size of particles in diatoms pool58 REAL(wp), PUBLIC :: xkr_dmicro !: Size of particles in microzoo pool59 REAL(wp), PUBLIC :: xkr_dmeso !: Size of particles in mesozoo pool60 REAL(wp), PUBLIC :: xkr_daggr !: Size of particles in aggregates pool61 REAL(wp), PUBLIC :: xkr_wsbio_min !: min vertical particle speed62 REAL(wp), PUBLIC :: xkr_wsbio_max !: max vertical particle speed63 64 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: xnumm !: maximum number of particles in aggregates65 #endif66 39 67 40 !!---------------------------------------------------------------------- … … 72 45 CONTAINS 73 46 74 #if ! defined key_kriest75 47 !!---------------------------------------------------------------------- 76 48 !! 'standard sinking parameterisation' ??? … … 91 63 REAL(wp) :: zagg1, zagg2, zagg3, zagg4 92 64 REAL(wp) :: zagg , zaggfe, zaggdoc, zaggdoc2, zaggdoc3 93 REAL(wp) :: zfact, zwsmax, zmax , zstep65 REAL(wp) :: zfact, zwsmax, zmax 94 66 CHARACTER (len=25) :: charout 95 67 REAL(wp), POINTER, DIMENSION(:,:,:) :: zw3d … … 98 70 ! 99 71 IF( nn_timing == 1 ) CALL timing_start('p4z_sink') 72 73 74 ! Initialization of some global variables 75 ! --------------------------------------- 76 prodpoc(:,:,:) = 0. 77 conspoc(:,:,:) = 0. 78 prodgoc(:,:,:) = 0. 79 consgoc(:,:,:) = 0. 80 100 81 ! 101 82 ! Sinking speeds of detritus is increased with depth as shown … … 105 86 DO jj = 1, jpj 106 87 DO ji = 1,jpi 107 zmax = MAX( heup (ji,jj), hmld(ji,jj) )108 zfact = MAX( 0., gdepw_n(ji,jj,jk+1) - zmax ) / 5000._wp109 wsbio4(ji,jj,jk) = wsbio2 + ( 200.- wsbio2) * zfact88 zmax = MAX( heup_01(ji,jj), hmld(ji,jj) ) 89 zfact = MAX( 0., gdepw_n(ji,jj,jk+1) - zmax ) / wsbio2scale 90 wsbio4(ji,jj,jk) = wsbio2 + MAX(0., ( wsbio2max - wsbio2 )) * zfact 110 91 END DO 111 92 END DO … … 114 95 ! limit the values of the sinking speeds to avoid numerical instabilities 115 96 wsbio3(:,:,:) = wsbio 116 wscal (:,:,:) = wsbio4(:,:,:) 97 117 98 ! 118 99 ! OA This is (I hope) a temporary solution for the problem that may … … 155 136 IF( tmask(ji,jj,jk) == 1 ) THEN 156 137 zwsmax = 0.5 * e3t_n(ji,jj,jk) / xstep 157 wsbio3(ji,jj,jk) = MIN( wsbio3(ji,jj,jk), zwsmax * FLOAT( iiter1) )158 wsbio4(ji,jj,jk) = MIN( wsbio4(ji,jj,jk), zwsmax * FLOAT( iiter2) )138 wsbio3(ji,jj,jk) = MIN( wsbio3(ji,jj,jk), zwsmax * REAL( iiter1, wp ) ) 139 wsbio4(ji,jj,jk) = MIN( wsbio4(ji,jj,jk), zwsmax * REAL( iiter2, wp ) ) 159 140 ENDIF 160 141 END DO 161 142 END DO 162 143 END DO 144 145 wscal (:,:,:) = wsbio4(:,:,:) 163 146 164 147 ! Initializa to zero all the sinking arrays … … 185 168 END DO 186 169 187 ! Exchange between organic matter compartments due to coagulation/disaggregation 188 ! --------------------------------------------------- 189 DO jk = 1, jpkm1 190 DO jj = 1, jpj 191 DO ji = 1, jpi 192 ! 193 zstep = xstep 194 # if defined key_degrad 195 zstep = zstep * facvol(ji,jj,jk) 196 # endif 197 zfact = zstep * xdiss(ji,jj,jk) 198 ! Part I : Coagulation dependent on turbulence 199 zagg1 = 25.9 * zfact * trb(ji,jj,jk,jppoc) * trb(ji,jj,jk,jppoc) 200 zagg2 = 4452. * zfact * trb(ji,jj,jk,jppoc) * trb(ji,jj,jk,jpgoc) 201 202 ! Part II : Differential settling 203 204 ! Aggregation of small into large particles 205 zagg3 = 47.1 * zstep * trb(ji,jj,jk,jppoc) * trb(ji,jj,jk,jpgoc) 206 zagg4 = 3.3 * zstep * trb(ji,jj,jk,jppoc) * trb(ji,jj,jk,jppoc) 207 208 zagg = zagg1 + zagg2 + zagg3 + zagg4 209 zaggfe = zagg * trb(ji,jj,jk,jpsfe) / ( trb(ji,jj,jk,jppoc) + rtrn ) 210 211 ! Aggregation of DOC to POC : 212 ! 1st term is shear aggregation of DOC-DOC 213 ! 2nd term is shear aggregation of DOC-POC 214 ! 3rd term is differential settling of DOC-POC 215 zaggdoc = ( ( 0.369 * 0.3 * trb(ji,jj,jk,jpdoc) + 102.4 * trb(ji,jj,jk,jppoc) ) * zfact & 216 & + 2.4 * zstep * trb(ji,jj,jk,jppoc) ) * 0.3 * trb(ji,jj,jk,jpdoc) 217 ! transfer of DOC to GOC : 218 ! 1st term is shear aggregation 219 ! 2nd term is differential settling 220 zaggdoc2 = ( 3.53E3 * zfact + 0.1 * zstep ) * trb(ji,jj,jk,jpgoc) * 0.3 * trb(ji,jj,jk,jpdoc) 221 ! tranfer of DOC to POC due to brownian motion 222 zaggdoc3 = ( 5095. * trb(ji,jj,jk,jppoc) + 114. * 0.3 * trb(ji,jj,jk,jpdoc) ) *zstep * 0.3 * trb(ji,jj,jk,jpdoc) 223 224 ! Update the trends 225 tra(ji,jj,jk,jppoc) = tra(ji,jj,jk,jppoc) - zagg + zaggdoc + zaggdoc3 226 tra(ji,jj,jk,jpgoc) = tra(ji,jj,jk,jpgoc) + zagg + zaggdoc2 227 tra(ji,jj,jk,jpsfe) = tra(ji,jj,jk,jpsfe) - zaggfe 228 tra(ji,jj,jk,jpbfe) = tra(ji,jj,jk,jpbfe) + zaggfe 229 tra(ji,jj,jk,jpdoc) = tra(ji,jj,jk,jpdoc) - zaggdoc - zaggdoc2 - zaggdoc3 230 ! 231 END DO 232 END DO 233 END DO 234 170 IF( ln_p5z ) THEN 171 sinkingn (:,:,:) = 0.e0 172 sinking2n(:,:,:) = 0.e0 173 sinkingp (:,:,:) = 0.e0 174 sinking2p(:,:,:) = 0.e0 175 176 ! Compute the sedimentation term using p4zsink2 for all the sinking particles 177 ! ----------------------------------------------------- 178 DO jit = 1, iiter1 179 CALL p4z_sink2( wsbio3, sinkingn , jppon, iiter1 ) 180 CALL p4z_sink2( wsbio3, sinkingp , jppop, iiter1 ) 181 END DO 182 183 DO jit = 1, iiter2 184 CALL p4z_sink2( wsbio4, sinking2n, jpgon, iiter2 ) 185 CALL p4z_sink2( wsbio4, sinking2p, jpgop, iiter2 ) 186 END DO 187 ENDIF 188 189 IF( ln_ligand ) THEN 190 wsfep (:,:,:) = wfep 191 DO jk = 1,jpkm1 192 DO jj = 1, jpj 193 DO ji = 1, jpi 194 IF( tmask(ji,jj,jk) == 1 ) THEN 195 zwsmax = 0.5 * e3t_n(ji,jj,jk) / xstep 196 wsfep(ji,jj,jk) = MIN( wsfep(ji,jj,jk), zwsmax * REAL( iiter1, wp ) ) 197 ENDIF 198 END DO 199 END DO 200 END DO 201 ! 202 sinkfep(:,:,:) = 0.e0 203 DO jit = 1, iiter1 204 CALL p4z_sink2( wsfep, sinkfep , jpfep, iiter1 ) 205 END DO 206 ENDIF 235 207 236 208 ! Total carbon export per year … … 281 253 CALL wrk_dealloc( jpi, jpj, jpk, zw3d ) 282 254 ENDIF 283 ELSE284 IF( ln_diatrc ) THEN285 zfact = 1.e3 * rfact2r286 trc2d(:,:,jp_pcs0_2d + 4) = sinking (:,:,ik100) * zfact * tmask(:,:,1)287 trc2d(:,:,jp_pcs0_2d + 5) = sinking2(:,:,ik100) * zfact * tmask(:,:,1)288 trc2d(:,:,jp_pcs0_2d + 6) = sinkfer (:,:,ik100) * zfact * tmask(:,:,1)289 trc2d(:,:,jp_pcs0_2d + 7) = sinkfer2(:,:,ik100) * zfact * tmask(:,:,1)290 trc2d(:,:,jp_pcs0_2d + 8) = sinksil (:,:,ik100) * zfact * tmask(:,:,1)291 trc2d(:,:,jp_pcs0_2d + 9) = sinkcal (:,:,ik100) * zfact * tmask(:,:,1)292 ENDIF293 255 ENDIF 294 256 ! … … 320 282 ! 321 283 END SUBROUTINE p4z_sink_init 322 323 #else324 !!----------------------------------------------------------------------325 !! 'Kriest sinking parameterisation' key_kriest ???326 !!----------------------------------------------------------------------327 328 SUBROUTINE p4z_sink ( kt, knt )329 !!---------------------------------------------------------------------330 !! *** ROUTINE p4z_sink ***331 !!332 !! ** Purpose : Compute vertical flux of particulate matter due to333 !! gravitational sinking - Kriest parameterization334 !!335 !! ** Method : - ???336 !!---------------------------------------------------------------------337 !338 INTEGER, INTENT(in) :: kt, knt339 !340 INTEGER :: ji, jj, jk, jit, niter1, niter2341 REAL(wp) :: zagg1, zagg2, zagg3, zagg4, zagg5, zfract, zaggsi, zaggsh342 REAL(wp) :: zagg , zaggdoc, zaggdoc1, znumdoc343 REAL(wp) :: znum , zeps, zfm, zgm, zsm344 REAL(wp) :: zdiv , zdiv1, zdiv2, zdiv3, zdiv4, zdiv5345 REAL(wp) :: zval1, zval2, zval3, zval4346 REAL(wp) :: zfact347 INTEGER :: ik1348 CHARACTER (len=25) :: charout349 REAL(wp), POINTER, DIMENSION(:,:,:) :: znum3d350 REAL(wp), POINTER, DIMENSION(:,:,:) :: zw3d351 REAL(wp), POINTER, DIMENSION(:,: ) :: zw2d352 !!---------------------------------------------------------------------353 !354 IF( nn_timing == 1 ) CALL timing_start('p4z_sink')355 !356 CALL wrk_alloc( jpi, jpj, jpk, znum3d )357 !358 ! Initialisation of variables used to compute Sinking Speed359 ! ---------------------------------------------------------360 361 znum3d(:,:,:) = 0.e0362 zval1 = 1. + xkr_zeta363 zval2 = 1. + xkr_zeta + xkr_eta364 zval3 = 1. + xkr_eta365 366 ! Computation of the vertical sinking speed : Kriest et Evans, 2000367 ! -----------------------------------------------------------------368 369 DO jk = 1, jpkm1370 DO jj = 1, jpj371 DO ji = 1, jpi372 IF( tmask(ji,jj,jk) /= 0.e0 ) THEN373 znum = trb(ji,jj,jk,jppoc) / ( trb(ji,jj,jk,jpnum) + rtrn ) / xkr_massp374 ! -------------- To avoid sinking speed over 50 m/day -------375 znum = MIN( xnumm(jk), znum )376 znum = MAX( 1.1 , znum )377 znum3d(ji,jj,jk) = znum378 !------------------------------------------------------------379 zeps = ( zval1 * znum - 1. )/ ( znum - 1. )380 zfm = xkr_frac**( 1. - zeps )381 zgm = xkr_frac**( zval1 - zeps )382 zdiv = MAX( 1.e-4, ABS( zeps - zval2 ) ) * SIGN( 1., ( zeps - zval2 ) )383 zdiv1 = zeps - zval3384 wsbio3(ji,jj,jk) = xkr_wsbio_min * ( zeps - zval1 ) / zdiv &385 & - xkr_wsbio_max * zgm * xkr_eta / zdiv386 wsbio4(ji,jj,jk) = xkr_wsbio_min * ( zeps-1. ) / zdiv1 &387 & - xkr_wsbio_max * zfm * xkr_eta / zdiv1388 IF( znum == 1.1) wsbio3(ji,jj,jk) = wsbio4(ji,jj,jk)389 ENDIF390 END DO391 END DO392 END DO393 394 wscal(:,:,:) = MAX( wsbio3(:,:,:), 30._wp )395 396 ! INITIALIZE TO ZERO ALL THE SINKING ARRAYS397 ! -----------------------------------------398 399 sinking (:,:,:) = 0.e0400 sinking2(:,:,:) = 0.e0401 sinkcal (:,:,:) = 0.e0402 sinkfer (:,:,:) = 0.e0403 sinksil (:,:,:) = 0.e0404 405 ! Compute the sedimentation term using p4zsink2 for all the sinking particles406 ! -----------------------------------------------------407 408 niter1 = niter1max409 niter2 = niter2max410 411 DO jit = 1, niter1412 CALL p4z_sink2( wsbio3, sinking , jppoc, niter1 )413 CALL p4z_sink2( wsbio3, sinkfer , jpsfe, niter1 )414 CALL p4z_sink2( wscal , sinksil , jpgsi, niter1 )415 CALL p4z_sink2( wscal , sinkcal , jpcal, niter1 )416 END DO417 418 DO jit = 1, niter2419 CALL p4z_sink2( wsbio4, sinking2, jpnum, niter2 )420 END DO421 422 ! Exchange between organic matter compartments due to coagulation/disaggregation423 ! ---------------------------------------------------424 425 zval1 = 1. + xkr_zeta426 zval2 = 1. + xkr_eta427 zval3 = 3. + xkr_eta428 zval4 = 4. + xkr_eta429 430 DO jk = 1,jpkm1431 DO jj = 1,jpj432 DO ji = 1,jpi433 IF( tmask(ji,jj,jk) /= 0.e0 ) THEN434 435 znum = trb(ji,jj,jk,jppoc)/(trb(ji,jj,jk,jpnum)+rtrn) / xkr_massp436 !-------------- To avoid sinking speed over 50 m/day -------437 znum = min(xnumm(jk),znum)438 znum = MAX( 1.1,znum)439 !------------------------------------------------------------440 zeps = ( zval1 * znum - 1.) / ( znum - 1.)441 zdiv = MAX( 1.e-4, ABS( zeps - zval3) ) * SIGN( 1., zeps - zval3 )442 zdiv1 = MAX( 1.e-4, ABS( zeps - 4. ) ) * SIGN( 1., zeps - 4. )443 zdiv2 = zeps - 2.444 zdiv3 = zeps - 3.445 zdiv4 = zeps - zval2446 zdiv5 = 2.* zeps - zval4447 zfm = xkr_frac**( 1.- zeps )448 zsm = xkr_frac**xkr_eta449 450 ! Part I : Coagulation dependant on turbulence451 ! ----------------------------------------------452 453 zagg1 = 0.163 * trb(ji,jj,jk,jpnum)**2 &454 & * 2.*( (zfm-1.)*(zfm*xkr_mass_max**3-xkr_mass_min**3) &455 & * (zeps-1)/zdiv1 + 3.*(zfm*xkr_mass_max-xkr_mass_min) &456 & * (zfm*xkr_mass_max**2-xkr_mass_min**2) &457 & * (zeps-1.)**2/(zdiv2*zdiv3))458 zagg2 = 2*0.163*trb(ji,jj,jk,jpnum)**2*zfm* &459 & ((xkr_mass_max**3+3.*(xkr_mass_max**2 &460 & *xkr_mass_min*(zeps-1.)/zdiv2 &461 & +xkr_mass_max*xkr_mass_min**2*(zeps-1.)/zdiv3) &462 & +xkr_mass_min**3*(zeps-1)/zdiv1) &463 & -zfm*xkr_mass_max**3*(1.+3.*((zeps-1.)/ &464 & (zeps-2.)+(zeps-1.)/zdiv3)+(zeps-1.)/zdiv1))465 466 zagg3 = 0.163*trb(ji,jj,jk,jpnum)**2*zfm**2*8. * xkr_mass_max**3467 468 ! Aggregation of small into large particles469 ! Part II : Differential settling470 ! ----------------------------------------------471 472 zagg4 = 2.*3.141*0.125*trb(ji,jj,jk,jpnum)**2* &473 & xkr_wsbio_min*(zeps-1.)**2 &474 & *(xkr_mass_min**2*((1.-zsm*zfm)/(zdiv3*zdiv4) &475 & -(1.-zfm)/(zdiv*(zeps-1.)))- &476 & ((zfm*zfm*xkr_mass_max**2*zsm-xkr_mass_min**2) &477 & *xkr_eta)/(zdiv*zdiv3*zdiv5) )478 479 zagg5 = 2.*3.141*0.125*trb(ji,jj,jk,jpnum)**2 &480 & *(zeps-1.)*zfm*xkr_wsbio_min &481 & *(zsm*(xkr_mass_min**2-zfm*xkr_mass_max**2) &482 & /zdiv3-(xkr_mass_min**2-zfm*zsm*xkr_mass_max**2) &483 & /zdiv)484 485 !486 ! Fractionnation by swimming organisms487 ! ------------------------------------488 489 zfract = 2.*3.141*0.125*trb(ji,jj,jk,jpmes)*12./0.12/0.06**3*trb(ji,jj,jk,jpnum) &490 & * (0.01/xkr_mass_min)**(1.-zeps)*0.1**2 &491 & * 10000.*xstep492 493 ! Aggregation of DOC to small particles494 ! --------------------------------------495 496 zaggdoc = 0.83 * trb(ji,jj,jk,jpdoc) * xstep * xdiss(ji,jj,jk) * trb(ji,jj,jk,jpdoc) &497 & + 0.005 * 231. * trb(ji,jj,jk,jpdoc) * xstep * trb(ji,jj,jk,jpdoc)498 zaggdoc1 = 271. * trb(ji,jj,jk,jppoc) * xstep * xdiss(ji,jj,jk) * trb(ji,jj,jk,jpdoc) &499 & + 0.02 * 16706. * trb(ji,jj,jk,jppoc) * xstep * trb(ji,jj,jk,jpdoc)500 501 # if defined key_degrad502 zagg1 = zagg1 * facvol(ji,jj,jk)503 zagg2 = zagg2 * facvol(ji,jj,jk)504 zagg3 = zagg3 * facvol(ji,jj,jk)505 zagg4 = zagg4 * facvol(ji,jj,jk)506 zagg5 = zagg5 * facvol(ji,jj,jk)507 zaggdoc = zaggdoc * facvol(ji,jj,jk)508 zaggdoc1 = zaggdoc1 * facvol(ji,jj,jk)509 # endif510 zaggsh = ( zagg1 + zagg2 + zagg3 ) * rfact2 * xdiss(ji,jj,jk) / 1000.511 zaggsi = ( zagg4 + zagg5 ) * xstep / 10.512 zagg = 0.5 * xkr_stick * ( zaggsh + zaggsi )513 !514 znumdoc = trb(ji,jj,jk,jpnum) / ( trb(ji,jj,jk,jppoc) + rtrn )515 tra(ji,jj,jk,jppoc) = tra(ji,jj,jk,jppoc) + zaggdoc + zaggdoc1516 tra(ji,jj,jk,jpnum) = tra(ji,jj,jk,jpnum) + zfract + zaggdoc / xkr_massp - zagg517 tra(ji,jj,jk,jpdoc) = tra(ji,jj,jk,jpdoc) - zaggdoc - zaggdoc1518 519 ENDIF520 END DO521 END DO522 END DO523 524 ! Total primary production per year525 t_oce_co2_exp = t_oce_co2_exp + glob_sum( ( sinking(:,:,ik100) * e1e2t(:,:) * tmask(:,:,1) )526 !527 IF( lk_iomput ) THEN528 IF( knt == nrdttrc ) THEN529 CALL wrk_alloc( jpi, jpj, zw2d )530 CALL wrk_alloc( jpi, jpj, jpk, zw3d )531 zfact = 1.e+3 * rfact2r ! conversion from mol/l/kt to mol/m3/s532 !533 IF( iom_use( "EPC100" ) ) THEN534 zw2d(:,:) = sinking(:,:,ik100) * zfact * tmask(:,:,1) ! Export of carbon at 100m535 CALL iom_put( "EPC100" , zw2d )536 ENDIF537 IF( iom_use( "EPN100" ) ) THEN538 zw2d(:,:) = sinking2(:,:,ik100) * zfact * tmask(:,:,1) ! Export of number of aggregates ?539 CALL iom_put( "EPN100" , zw2d )540 ENDIF541 IF( iom_use( "EPCAL100" ) ) THEN542 zw2d(:,:) = sinkcal(:,:,ik100) * zfact * tmask(:,:,1) ! Export of calcite at 100m543 CALL iom_put( "EPCAL100" , zw2d )544 ENDIF545 IF( iom_use( "EPSI100" ) ) THEN546 zw2d(:,:) = sinksil(:,:,ik100) * zfact * tmask(:,:,1) ! Export of bigenic silica at 100m547 CALL iom_put( "EPSI100" , zw2d )548 ENDIF549 IF( iom_use( "EXPC" ) ) THEN550 zw3d(:,:,:) = sinking(:,:,:) * zfact * tmask(:,:,:) ! Export of carbon in the water column551 CALL iom_put( "EXPC" , zw3d )552 ENDIF553 IF( iom_use( "EXPN" ) ) THEN554 zw3d(:,:,:) = sinking(:,:,:) * zfact * tmask(:,:,:) ! Export of carbon in the water column555 CALL iom_put( "EXPN" , zw3d )556 ENDIF557 IF( iom_use( "EXPCAL" ) ) THEN558 zw3d(:,:,:) = sinkcal(:,:,:) * zfact * tmask(:,:,:) ! Export of calcite559 CALL iom_put( "EXPCAL" , zw3d )560 ENDIF561 IF( iom_use( "EXPSI" ) ) THEN562 zw3d(:,:,:) = sinksil(:,:,:) * zfact * tmask(:,:,:) ! Export of bigenic silica563 CALL iom_put( "EXPSI" , zw3d )564 ENDIF565 IF( iom_use( "XNUM" ) ) THEN566 zw3d(:,:,:) = znum3d(:,:,:) * tmask(:,:,:) ! Number of particles on aggregats567 CALL iom_put( "XNUM" , zw3d )568 ENDIF569 IF( iom_use( "WSC" ) ) THEN570 zw3d(:,:,:) = wsbio3(:,:,:) * tmask(:,:,:) ! Sinking speed of carbon particles571 CALL iom_put( "WSC" , zw3d )572 ENDIF573 IF( iom_use( "WSN" ) ) THEN574 zw3d(:,:,:) = wsbio4(:,:,:) * tmask(:,:,:) ! Sinking speed of particles number575 CALL iom_put( "WSN" , zw3d )576 ENDIF577 !578 CALL wrk_dealloc( jpi, jpj, zw2d )579 CALL wrk_dealloc( jpi, jpj, jpk, zw3d )580 ELSE581 IF( ln_diatrc ) THEN582 zfact = 1.e3 * rfact2r583 trc2d(:,: ,jp_pcs0_2d + 4) = sinking (:,:,ik100) * zfact * tmask(:,:,1)584 trc2d(:,: ,jp_pcs0_2d + 5) = sinking2(:,:,ik100) * zfact * tmask(:,:,1)585 trc2d(:,: ,jp_pcs0_2d + 6) = sinkfer (:,:,ik100) * zfact * tmask(:,:,1)586 trc2d(:,: ,jp_pcs0_2d + 7) = sinksil (:,:,ik100) * zfact * tmask(:,:,1)587 trc2d(:,: ,jp_pcs0_2d + 8) = sinkcal (:,:,ik100) * zfact * tmask(:,:,1)588 trc3d(:,:,:,jp_pcs0_3d + 11) = sinking (:,:,:) * zfact * tmask(:,:,:)589 trc3d(:,:,:,jp_pcs0_3d + 12) = sinking2(:,:,:) * zfact * tmask(:,:,:)590 trc3d(:,:,:,jp_pcs0_3d + 13) = sinksil (:,:,:) * zfact * tmask(:,:,:)591 trc3d(:,:,:,jp_pcs0_3d + 14) = sinkcal (:,:,:) * zfact * tmask(:,:,:)592 trc3d(:,:,:,jp_pcs0_3d + 15) = znum3d (:,:,:) * tmask(:,:,:)593 trc3d(:,:,:,jp_pcs0_3d + 16) = wsbio3 (:,:,:) * tmask(:,:,:)594 trc3d(:,:,:,jp_pcs0_3d + 17) = wsbio4 (:,:,:) * tmask(:,:,:)595 ENDIF596 ENDIF597 598 !599 IF(ln_ctl) THEN ! print mean trends (used for debugging)600 WRITE(charout, FMT="('sink')")601 CALL prt_ctl_trc_info(charout)602 CALL prt_ctl_trc(tab4d=tra, mask=tmask, clinfo=ctrcnm)603 ENDIF604 !605 CALL wrk_dealloc( jpi, jpj, jpk, znum3d )606 !607 IF( nn_timing == 1 ) CALL timing_stop('p4z_sink')608 !609 END SUBROUTINE p4z_sink610 611 612 SUBROUTINE p4z_sink_init613 !!----------------------------------------------------------------------614 !! *** ROUTINE p4z_sink_init ***615 !!616 !! ** Purpose : Initialization of sinking parameters617 !! Kriest parameterization only618 !!619 !! ** Method : Read the nampiskrs namelist and check the parameters620 !! called at the first timestep621 !!622 !! ** input : Namelist nampiskrs623 !!----------------------------------------------------------------------624 INTEGER :: jk, jn, kiter625 INTEGER :: ios ! Local integer output status for namelist read626 REAL(wp) :: znum, zdiv627 REAL(wp) :: zws, zwr, zwl,wmax, znummax628 REAL(wp) :: zmin, zmax, zl, zr, xacc629 !630 NAMELIST/nampiskrs/ xkr_sfact, xkr_stick , &631 & xkr_nnano, xkr_ndiat, xkr_nmicro, xkr_nmeso, xkr_naggr632 !!----------------------------------------------------------------------633 !634 IF( nn_timing == 1 ) CALL timing_start('p4z_sink_init')635 !636 637 REWIND( numnatp_ref ) ! Namelist nampiskrs in reference namelist : Pisces sinking Kriest638 READ ( numnatp_ref, nampiskrs, IOSTAT = ios, ERR = 901)639 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nampiskrs in reference namelist', lwp )640 641 REWIND( numnatp_cfg ) ! Namelist nampiskrs in configuration namelist : Pisces sinking Kriest642 READ ( numnatp_cfg, nampiskrs, IOSTAT = ios, ERR = 902 )643 902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nampiskrs in configuration namelist', lwp )644 IF(lwm) WRITE ( numonp, nampiskrs )645 646 IF(lwp) THEN647 WRITE(numout,*)648 WRITE(numout,*) ' Namelist : nampiskrs'649 WRITE(numout,*) ' Sinking factor xkr_sfact = ', xkr_sfact650 WRITE(numout,*) ' Stickiness xkr_stick = ', xkr_stick651 WRITE(numout,*) ' Nbr of cell in nano size class xkr_nnano = ', xkr_nnano652 WRITE(numout,*) ' Nbr of cell in diatoms size class xkr_ndiat = ', xkr_ndiat653 WRITE(numout,*) ' Nbr of cell in microzoo size class xkr_nmicro = ', xkr_nmicro654 WRITE(numout,*) ' Nbr of cell in mesozoo size class xkr_nmeso = ', xkr_nmeso655 WRITE(numout,*) ' Nbr of cell in aggregates size class xkr_naggr = ', xkr_naggr656 ENDIF657 658 659 ! max and min vertical particle speed660 xkr_wsbio_min = xkr_sfact * xkr_mass_min**xkr_eta661 xkr_wsbio_max = xkr_sfact * xkr_mass_max**xkr_eta662 IF (lwp) WRITE(numout,*) ' max and min vertical particle speed ', xkr_wsbio_min, xkr_wsbio_max663 664 !665 ! effect of the sizes of the different living pools on particle numbers666 ! nano = 2um-20um -> mean size=6.32 um -> ws=2.596 -> xnum=xnnano=2.337667 ! diat and microzoo = 10um-200um -> 44.7 -> 8.732 -> xnum=xndiat=3.718668 ! mesozoo = 200um-2mm -> 632.45 -> 45.14 -> xnum=xnmeso=7.147669 ! aggregates = 200um-10mm -> 1414 -> 74.34 -> xnum=xnaggr=9.877670 ! doc aggregates = 1um671 ! ----------------------------------------------------------672 673 xkr_dnano = 1. / ( xkr_massp * xkr_nnano )674 xkr_ddiat = 1. / ( xkr_massp * xkr_ndiat )675 xkr_dmicro = 1. / ( xkr_massp * xkr_nmicro )676 xkr_dmeso = 1. / ( xkr_massp * xkr_nmeso )677 xkr_daggr = 1. / ( xkr_massp * xkr_naggr )678 679 !!---------------------------------------------------------------------680 !! 'key_kriest' ???681 !!---------------------------------------------------------------------682 ! COMPUTATION OF THE VERTICAL PROFILE OF MAXIMUM SINKING SPEED683 ! Search of the maximum number of particles in aggregates for each k-level.684 ! Bissection Method685 !--------------------------------------------------------------------686 IF (lwp) THEN687 WRITE(numout,*)688 WRITE(numout,*)' kriest : Compute maximum number of particles in aggregates'689 ENDIF690 691 xacc = 0.001_wp692 kiter = 50693 zmin = 1.10_wp694 zmax = xkr_mass_max / xkr_mass_min695 xkr_frac = zmax696 697 DO jk = 1,jpk698 zl = zmin699 zr = zmax700 wmax = 0.5 * e3t_n(1,1,jk) * rday * float(niter1max) / rfact2701 zdiv = xkr_zeta + xkr_eta - xkr_eta * zl702 znum = zl - 1.703 zwl = xkr_wsbio_min * xkr_zeta / zdiv &704 & - ( xkr_wsbio_max * xkr_eta * znum * &705 & xkr_frac**( -xkr_zeta / znum ) / zdiv ) &706 & - wmax707 708 zdiv = xkr_zeta + xkr_eta - xkr_eta * zr709 znum = zr - 1.710 zwr = xkr_wsbio_min * xkr_zeta / zdiv &711 & - ( xkr_wsbio_max * xkr_eta * znum * &712 & xkr_frac**( -xkr_zeta / znum ) / zdiv ) &713 & - wmax714 iflag: DO jn = 1, kiter715 IF ( zwl == 0._wp ) THEN ; znummax = zl716 ELSEIF( zwr == 0._wp ) THEN ; znummax = zr717 ELSE718 znummax = ( zr + zl ) / 2.719 zdiv = xkr_zeta + xkr_eta - xkr_eta * znummax720 znum = znummax - 1.721 zws = xkr_wsbio_min * xkr_zeta / zdiv &722 & - ( xkr_wsbio_max * xkr_eta * znum * &723 & xkr_frac**( -xkr_zeta / znum ) / zdiv ) &724 & - wmax725 IF( zws * zwl < 0. ) THEN ; zr = znummax726 ELSE ; zl = znummax727 ENDIF728 zdiv = xkr_zeta + xkr_eta - xkr_eta * zl729 znum = zl - 1.730 zwl = xkr_wsbio_min * xkr_zeta / zdiv &731 & - ( xkr_wsbio_max * xkr_eta * znum * &732 & xkr_frac**( -xkr_zeta / znum ) / zdiv ) &733 & - wmax734 735 zdiv = xkr_zeta + xkr_eta - xkr_eta * zr736 znum = zr - 1.737 zwr = xkr_wsbio_min * xkr_zeta / zdiv &738 & - ( xkr_wsbio_max * xkr_eta * znum * &739 & xkr_frac**( -xkr_zeta / znum ) / zdiv ) &740 & - wmax741 !742 IF ( ABS ( zws ) <= xacc ) EXIT iflag743 !744 ENDIF745 !746 END DO iflag747 748 xnumm(jk) = znummax749 IF (lwp) WRITE(numout,*) ' jk = ', jk, ' wmax = ', wmax,' xnum max = ', xnumm(jk)750 !751 END DO752 !753 ik100 = 10 ! last level where depth less than 100 m754 DO jk = jpkm1, 1, -1755 IF( gdept_1d(jk) > 100. ) iksed = jk - 1756 END DO757 IF (lwp) WRITE(numout,*)758 IF (lwp) WRITE(numout,*) ' Level corresponding to 100m depth ', ik100 + 1759 IF (lwp) WRITE(numout,*)760 !761 t_oce_co2_exp = 0._wp762 !763 IF( nn_timing == 1 ) CALL timing_stop('p4z_sink_init')764 !765 END SUBROUTINE p4z_sink_init766 767 #endif768 284 769 285 SUBROUTINE p4z_sink2( pwsink, psinkflx, jp_tra, kiter ) … … 794 310 CALL wrk_alloc( jpi, jpj, jpk, ztraz, zakz, zwsink2, ztrb ) 795 311 796 zstep = rfact2 / FLOAT( kiter) / 2.312 zstep = rfact2 / REAL( kiter, wp ) / 2. 797 313 798 314 ztraz(:,:,:) = 0.e0 … … 804 320 END DO 805 321 zwsink2(:,:,1) = 0.e0 806 IF( lk_degrad ) THEN807 zwsink2(:,:,:) = zwsink2(:,:,:) * facvol(:,:,:)808 ENDIF809 322 810 323 … … 887 400 !! *** ROUTINE p4z_sink_alloc *** 888 401 !!---------------------------------------------------------------------- 889 ALLOCATE( wsbio3 (jpi,jpj,jpk) , wsbio4 (jpi,jpj,jpk) , wscal(jpi,jpj,jpk) , & 890 & sinking(jpi,jpj,jpk) , sinking2(jpi,jpj,jpk) , & 891 & sinkcal(jpi,jpj,jpk) , sinksil (jpi,jpj,jpk) , & 892 #if defined key_kriest 893 & xnumm(jpk) , & 894 #else 895 & sinkfer2(jpi,jpj,jpk) , & 896 #endif 897 & sinkfer(jpi,jpj,jpk) , STAT=p4z_sink_alloc ) 402 INTEGER :: ierr(3) 403 404 ierr(:) = 0 405 ! 406 ALLOCATE( sinking(jpi,jpj,jpk) , sinking2(jpi,jpj,jpk) , & 407 & sinkcal(jpi,jpj,jpk) , sinksil (jpi,jpj,jpk) , & 408 & sinkfer2(jpi,jpj,jpk) , & 409 & sinkfer(jpi,jpj,jpk) , STAT=ierr(1) ) 898 410 ! 411 IF( ln_ligand ) ALLOCATE( sinkfep(jpi,jpj,jpk) , STAT=ierr(2) ) 412 413 IF( ln_p5z ) ALLOCATE( sinkingn(jpi,jpj,jpk), sinking2n(jpi,jpj,jpk) , & 414 & sinkingp(jpi,jpj,jpk), sinking2p(jpi,jpj,jpk) , STAT=ierr(3) ) 415 ! 416 p4z_sink_alloc = MAXVAL( ierr ) 899 417 IF( p4z_sink_alloc /= 0 ) CALL ctl_warn('p4z_sink_alloc : failed to allocate arrays.') 900 418 ! 901 419 END FUNCTION p4z_sink_alloc 902 420 903 #else904 !!======================================================================905 !! Dummy module : No PISCES bio-model906 !!======================================================================907 CONTAINS908 SUBROUTINE p4z_sink ! Empty routine909 END SUBROUTINE p4z_sink910 #endif911 912 421 !!====================================================================== 913 422 END MODULE p4zsink
Note: See TracChangeset
for help on using the changeset viewer.
