- Timestamp:
- 2017-05-09T12:14:45+02:00 (7 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zsink.F90
r7180 r8003 21 21 USE iom ! I/O manager 22 22 USE lib_mpp 23 USE p4zsbc 23 24 24 25 IMPLICIT NONE … … 37 38 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: sinkcal, sinksil !: CaCO3 and BSi sinking fluxes 38 39 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: sinkfer !: Small BFe sinking fluxes 39 #if ! defined key_kriest40 40 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: sinkfer2 !: Big iron sinking fluxes 41 #endif42 41 #if defined key_ligand 43 42 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: sinkfep !: Fep sinking fluxes … … 46 45 47 46 INTEGER :: ik100 48 49 #if defined key_kriest50 REAL(wp) :: xkr_sfact !: Sinking factor51 REAL(wp) :: xkr_stick !: Stickiness52 REAL(wp) :: xkr_nnano !: Nbr of cell in nano size class53 REAL(wp) :: xkr_ndiat !: Nbr of cell in diatoms size class54 REAL(wp) :: xkr_nmicro !: Nbr of cell in microzoo size class55 REAL(wp) :: xkr_nmeso !: Nbr of cell in mesozoo size class56 REAL(wp) :: xkr_naggr !: Nbr of cell in aggregates size class57 58 REAL(wp) :: xkr_frac59 60 REAL(wp), PUBLIC :: xkr_dnano !: Size of particles in nano pool61 REAL(wp), PUBLIC :: xkr_ddiat !: Size of particles in diatoms pool62 REAL(wp), PUBLIC :: xkr_dmicro !: Size of particles in microzoo pool63 REAL(wp), PUBLIC :: xkr_dmeso !: Size of particles in mesozoo pool64 REAL(wp), PUBLIC :: xkr_daggr !: Size of particles in aggregates pool65 REAL(wp), PUBLIC :: xkr_wsbio_min !: min vertical particle speed66 REAL(wp), PUBLIC :: xkr_wsbio_max !: max vertical particle speed67 68 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: xnumm !: maximum number of particles in aggregates69 #endif70 47 71 48 !!* Substitution … … 78 55 CONTAINS 79 56 80 #if ! defined key_kriest81 57 !!---------------------------------------------------------------------- 82 58 !! 'standard sinking parameterisation' ??? … … 255 231 CALL wrk_dealloc( jpi, jpj, jpk, zw3d ) 256 232 ENDIF 257 ELSE258 IF( ln_diatrc ) THEN259 zfact = 1.e3 * rfact2r260 trc2d(:,:,jp_pcs0_2d + 4) = sinking (:,:,ik100) * zfact * tmask(:,:,1)261 trc2d(:,:,jp_pcs0_2d + 5) = sinking2(:,:,ik100) * zfact * tmask(:,:,1)262 trc2d(:,:,jp_pcs0_2d + 6) = sinkfer (:,:,ik100) * zfact * tmask(:,:,1)263 trc2d(:,:,jp_pcs0_2d + 7) = sinkfer2(:,:,ik100) * zfact * tmask(:,:,1)264 trc2d(:,:,jp_pcs0_2d + 8) = sinksil (:,:,ik100) * zfact * tmask(:,:,1)265 trc2d(:,:,jp_pcs0_2d + 9) = sinkcal (:,:,ik100) * zfact * tmask(:,:,1)266 ENDIF267 233 ENDIF 268 234 ! … … 294 260 ! 295 261 END SUBROUTINE p4z_sink_init 296 297 #else298 !!----------------------------------------------------------------------299 !! 'Kriest sinking parameterisation' key_kriest ???300 !!----------------------------------------------------------------------301 302 SUBROUTINE p4z_sink ( kt, knt )303 !!---------------------------------------------------------------------304 !! *** ROUTINE p4z_sink ***305 !!306 !! ** Purpose : Compute vertical flux of particulate matter due to307 !! gravitational sinking - Kriest parameterization308 !!309 !! ** Method : - ???310 !!---------------------------------------------------------------------311 !312 INTEGER, INTENT(in) :: kt, knt313 !314 INTEGER :: ji, jj, jk, jit, niter1, niter2315 REAL(wp) :: znum , zeps, zfm, zgm, zsm316 REAL(wp) :: zdiv , zdiv1, zdiv2, zdiv3, zdiv4, zdiv5317 REAL(wp) :: zval1, zval2, zval3318 REAL(wp) :: zfact319 INTEGER :: ik1320 CHARACTER (len=25) :: charout321 REAL(wp), POINTER, DIMENSION(:,:,:) :: znum3d322 REAL(wp), POINTER, DIMENSION(:,:,:) :: zw3d323 REAL(wp), POINTER, DIMENSION(:,: ) :: zw2d324 !!---------------------------------------------------------------------325 !326 IF( nn_timing == 1 ) CALL timing_start('p4z_sink')327 !328 CALL wrk_alloc( jpi, jpj, jpk, znum3d )329 !330 ! Initialisation of variables used to compute Sinking Speed331 ! ---------------------------------------------------------332 333 znum3d(:,:,:) = 0.e0334 zval1 = 1. + xkr_zeta335 zval2 = 1. + xkr_zeta + xkr_eta336 zval3 = 1. + xkr_eta337 338 ! Computation of the vertical sinking speed : Kriest et Evans, 2000339 ! -----------------------------------------------------------------340 341 DO jk = 1, jpkm1342 DO jj = 1, jpj343 DO ji = 1, jpi344 IF( tmask(ji,jj,jk) /= 0.e0 ) THEN345 znum = trb(ji,jj,jk,jppoc) / ( trb(ji,jj,jk,jpnum) + rtrn ) / xkr_massp346 ! -------------- To avoid sinking speed over 50 m/day -------347 znum = MIN( xnumm(jk), znum )348 znum = MAX( 1.1 , znum )349 znum3d(ji,jj,jk) = znum350 !------------------------------------------------------------351 zeps = ( zval1 * znum - 1. )/ ( znum - 1. )352 zfm = xkr_frac**( 1. - zeps )353 zgm = xkr_frac**( zval1 - zeps )354 zdiv = MAX( 1.e-4, ABS( zeps - zval2 ) ) * SIGN( 1., ( zeps - zval2 ) )355 zdiv1 = zeps - zval3356 wsbio3(ji,jj,jk) = xkr_wsbio_min * ( zeps - zval1 ) / zdiv &357 & - xkr_wsbio_max * zgm * xkr_eta / zdiv358 wsbio4(ji,jj,jk) = xkr_wsbio_min * ( zeps-1. ) / zdiv1 &359 & - xkr_wsbio_max * zfm * xkr_eta / zdiv1360 IF( znum == 1.1) wsbio3(ji,jj,jk) = wsbio4(ji,jj,jk)361 ENDIF362 END DO363 END DO364 END DO365 366 wscal(:,:,:) = MAX( wsbio3(:,:,:), 30._wp )367 #if defined key_ligand368 wsfep (:,:,:) = wfep369 #endif370 371 ! INITIALIZE TO ZERO ALL THE SINKING ARRAYS372 ! -----------------------------------------373 374 sinking (:,:,:) = 0.e0375 sinking2(:,:,:) = 0.e0376 sinkcal (:,:,:) = 0.e0377 sinkfer (:,:,:) = 0.e0378 sinksil (:,:,:) = 0.e0379 #if defined key_ligand380 sinkfep(:,:,:) = 0.e0381 #endif382 383 ! Compute the sedimentation term using p4zsink2 for all the sinking particles384 ! -----------------------------------------------------385 386 niter1 = niter1max387 niter2 = niter2max388 389 DO jit = 1, niter1390 CALL p4z_sink2( wsbio3, sinking , jppoc, niter1 )391 CALL p4z_sink2( wsbio3, sinkfer , jpsfe, niter1 )392 CALL p4z_sink2( wscal , sinksil , jpgsi, niter1 )393 CALL p4z_sink2( wscal , sinkcal , jpcal, niter1 )394 #if defined key_ligand395 CALL p4z_sink2( wsfep, sinkfep , jpfep, iiter1 )396 #endif397 END DO398 399 DO jit = 1, niter2400 CALL p4z_sink2( wsbio4, sinking2, jpnum, niter2 )401 END DO402 403 IF( iom_use( "tcexp" ) .OR. ( ln_check_mass .AND. kt == nitend .AND. knt == nrdttrc ) ) &404 & t_oce_co2_exp = glob_sum( sinking(:,:,ik100) * e1e2t(:,:) * tmask(:,:,1) )405 !406 IF( lk_iomput ) THEN407 IF( knt == nrdttrc ) THEN408 CALL wrk_alloc( jpi, jpj, zw2d )409 CALL wrk_alloc( jpi, jpj, jpk, zw3d )410 zfact = 1.e+3 * rfact2r ! conversion from mol/l/kt to mol/m3/s411 !412 IF( iom_use( "EPC100" ) ) THEN413 zw2d(:,:) = sinking(:,:,ik100) * zfact * tmask(:,:,1) ! Export of carbon at 100m414 CALL iom_put( "EPC100" , zw2d )415 ENDIF416 IF( iom_use( "EPN100" ) ) THEN417 zw2d(:,:) = sinking2(:,:,ik100) * zfact * tmask(:,:,1) ! Export of number of aggregates ?418 CALL iom_put( "EPN100" , zw2d )419 ENDIF420 IF( iom_use( "EPCAL100" ) ) THEN421 zw2d(:,:) = sinkcal(:,:,ik100) * zfact * tmask(:,:,1) ! Export of calcite at 100m422 CALL iom_put( "EPCAL100" , zw2d )423 ENDIF424 IF( iom_use( "EPSI100" ) ) THEN425 zw2d(:,:) = sinksil(:,:,ik100) * zfact * tmask(:,:,1) ! Export of bigenic silica at 100m426 CALL iom_put( "EPSI100" , zw2d )427 ENDIF428 IF( iom_use( "EXPC" ) ) THEN429 zw3d(:,:,:) = sinking(:,:,:) * zfact * tmask(:,:,:) ! Export of carbon in the water column430 CALL iom_put( "EXPC" , zw3d )431 ENDIF432 IF( iom_use( "EXPN" ) ) THEN433 zw3d(:,:,:) = sinking(:,:,:) * zfact * tmask(:,:,:) ! Export of carbon in the water column434 CALL iom_put( "EXPN" , zw3d )435 ENDIF436 IF( iom_use( "EXPCAL" ) ) THEN437 zw3d(:,:,:) = sinkcal(:,:,:) * zfact * tmask(:,:,:) ! Export of calcite438 CALL iom_put( "EXPCAL" , zw3d )439 ENDIF440 IF( iom_use( "EXPSI" ) ) THEN441 zw3d(:,:,:) = sinksil(:,:,:) * zfact * tmask(:,:,:) ! Export of bigenic silica442 CALL iom_put( "EXPSI" , zw3d )443 ENDIF444 IF( iom_use( "XNUM" ) ) THEN445 zw3d(:,:,:) = znum3d(:,:,:) * tmask(:,:,:) ! Number of particles on aggregats446 CALL iom_put( "XNUM" , zw3d )447 ENDIF448 IF( iom_use( "WSC" ) ) THEN449 zw3d(:,:,:) = wsbio3(:,:,:) * tmask(:,:,:) ! Sinking speed of carbon particles450 CALL iom_put( "WSC" , zw3d )451 ENDIF452 IF( iom_use( "WSN" ) ) THEN453 zw3d(:,:,:) = wsbio4(:,:,:) * tmask(:,:,:) ! Sinking speed of particles number454 CALL iom_put( "WSN" , zw3d )455 ENDIF456 !457 CALL wrk_dealloc( jpi, jpj, zw2d )458 CALL wrk_dealloc( jpi, jpj, jpk, zw3d )459 ELSE460 IF( ln_diatrc ) THEN461 zfact = 1.e3 * rfact2r462 trc2d(:,: ,jp_pcs0_2d + 4) = sinking (:,:,ik100) * zfact * tmask(:,:,1)463 trc2d(:,: ,jp_pcs0_2d + 5) = sinking2(:,:,ik100) * zfact * tmask(:,:,1)464 trc2d(:,: ,jp_pcs0_2d + 6) = sinkfer (:,:,ik100) * zfact * tmask(:,:,1)465 trc2d(:,: ,jp_pcs0_2d + 7) = sinksil (:,:,ik100) * zfact * tmask(:,:,1)466 trc2d(:,: ,jp_pcs0_2d + 8) = sinkcal (:,:,ik100) * zfact * tmask(:,:,1)467 trc3d(:,:,:,jp_pcs0_3d + 11) = sinking (:,:,:) * zfact * tmask(:,:,:)468 trc3d(:,:,:,jp_pcs0_3d + 12) = sinking2(:,:,:) * zfact * tmask(:,:,:)469 trc3d(:,:,:,jp_pcs0_3d + 13) = sinksil (:,:,:) * zfact * tmask(:,:,:)470 trc3d(:,:,:,jp_pcs0_3d + 14) = sinkcal (:,:,:) * zfact * tmask(:,:,:)471 trc3d(:,:,:,jp_pcs0_3d + 15) = znum3d (:,:,:) * tmask(:,:,:)472 trc3d(:,:,:,jp_pcs0_3d + 16) = wsbio3 (:,:,:) * tmask(:,:,:)473 trc3d(:,:,:,jp_pcs0_3d + 17) = wsbio4 (:,:,:) * tmask(:,:,:)474 ENDIF475 ENDIF476 477 !478 IF(ln_ctl) THEN ! print mean trends (used for debugging)479 WRITE(charout, FMT="('sink')")480 CALL prt_ctl_trc_info(charout)481 CALL prt_ctl_trc(tab4d=tra, mask=tmask, clinfo=ctrcnm)482 ENDIF483 !484 CALL wrk_dealloc( jpi, jpj, jpk, znum3d )485 !486 IF( nn_timing == 1 ) CALL timing_stop('p4z_sink')487 !488 END SUBROUTINE p4z_sink489 490 491 SUBROUTINE p4z_sink_init492 !!----------------------------------------------------------------------493 !! *** ROUTINE p4z_sink_init ***494 !!495 !! ** Purpose : Initialization of sinking parameters496 !! Kriest parameterization only497 !!498 !! ** Method : Read the nampiskrs namelist and check the parameters499 !! called at the first timestep500 !!501 !! ** input : Namelist nampiskrs502 !!----------------------------------------------------------------------503 INTEGER :: jk, jn, kiter504 INTEGER :: ios ! Local integer output status for namelist read505 REAL(wp) :: znum, zdiv506 REAL(wp) :: zws, zwr, zwl,wmax, znummax507 REAL(wp) :: zmin, zmax, zl, zr, xacc508 !509 NAMELIST/nampiskrs/ xkr_sfact, xkr_stick , &510 & xkr_nnano, xkr_ndiat, xkr_nmicro, xkr_nmeso, xkr_naggr511 !!----------------------------------------------------------------------512 !513 IF( nn_timing == 1 ) CALL timing_start('p4z_sink_init')514 !515 516 REWIND( numnatp_ref ) ! Namelist nampiskrs in reference namelist : Pisces sinking Kriest517 READ ( numnatp_ref, nampiskrs, IOSTAT = ios, ERR = 901)518 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nampiskrs in reference namelist', lwp )519 520 REWIND( numnatp_cfg ) ! Namelist nampiskrs in configuration namelist : Pisces sinking Kriest521 READ ( numnatp_cfg, nampiskrs, IOSTAT = ios, ERR = 902 )522 902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nampiskrs in configuration namelist', lwp )523 IF(lwm) WRITE ( numonp, nampiskrs )524 525 IF(lwp) THEN526 WRITE(numout,*)527 WRITE(numout,*) ' Namelist : nampiskrs'528 WRITE(numout,*) ' Sinking factor xkr_sfact = ', xkr_sfact529 WRITE(numout,*) ' Stickiness xkr_stick = ', xkr_stick530 WRITE(numout,*) ' Nbr of cell in nano size class xkr_nnano = ', xkr_nnano531 WRITE(numout,*) ' Nbr of cell in diatoms size class xkr_ndiat = ', xkr_ndiat532 WRITE(numout,*) ' Nbr of cell in microzoo size class xkr_nmicro = ', xkr_nmicro533 WRITE(numout,*) ' Nbr of cell in mesozoo size class xkr_nmeso = ', xkr_nmeso534 WRITE(numout,*) ' Nbr of cell in aggregates size class xkr_naggr = ', xkr_naggr535 ENDIF536 537 538 ! max and min vertical particle speed539 xkr_wsbio_min = xkr_sfact * xkr_mass_min**xkr_eta540 xkr_wsbio_max = xkr_sfact * xkr_mass_max**xkr_eta541 IF (lwp) WRITE(numout,*) ' max and min vertical particle speed ', xkr_wsbio_min, xkr_wsbio_max542 543 !544 ! effect of the sizes of the different living pools on particle numbers545 ! nano = 2um-20um -> mean size=6.32 um -> ws=2.596 -> xnum=xnnano=2.337546 ! diat and microzoo = 10um-200um -> 44.7 -> 8.732 -> xnum=xndiat=3.718547 ! mesozoo = 200um-2mm -> 632.45 -> 45.14 -> xnum=xnmeso=7.147548 ! aggregates = 200um-10mm -> 1414 -> 74.34 -> xnum=xnaggr=9.877549 ! doc aggregates = 1um550 ! ----------------------------------------------------------551 552 xkr_dnano = 1. / ( xkr_massp * xkr_nnano )553 xkr_ddiat = 1. / ( xkr_massp * xkr_ndiat )554 xkr_dmicro = 1. / ( xkr_massp * xkr_nmicro )555 xkr_dmeso = 1. / ( xkr_massp * xkr_nmeso )556 xkr_daggr = 1. / ( xkr_massp * xkr_naggr )557 558 !!---------------------------------------------------------------------559 !! 'key_kriest' ???560 !!---------------------------------------------------------------------561 ! COMPUTATION OF THE VERTICAL PROFILE OF MAXIMUM SINKING SPEED562 ! Search of the maximum number of particles in aggregates for each k-level.563 ! Bissection Method564 !--------------------------------------------------------------------565 IF (lwp) THEN566 WRITE(numout,*)567 WRITE(numout,*)' kriest : Compute maximum number of particles in aggregates'568 ENDIF569 570 xacc = 0.001_wp571 kiter = 50572 zmin = 1.10_wp573 zmax = xkr_mass_max / xkr_mass_min574 xkr_frac = zmax575 576 DO jk = 1,jpk577 zl = zmin578 zr = zmax579 wmax = 0.5 * fse3t(1,1,jk) * rday * float(niter1max) / rfact2580 zdiv = xkr_zeta + xkr_eta - xkr_eta * zl581 znum = zl - 1.582 zwl = xkr_wsbio_min * xkr_zeta / zdiv &583 & - ( xkr_wsbio_max * xkr_eta * znum * &584 & xkr_frac**( -xkr_zeta / znum ) / zdiv ) &585 & - wmax586 587 zdiv = xkr_zeta + xkr_eta - xkr_eta * zr588 znum = zr - 1.589 zwr = xkr_wsbio_min * xkr_zeta / zdiv &590 & - ( xkr_wsbio_max * xkr_eta * znum * &591 & xkr_frac**( -xkr_zeta / znum ) / zdiv ) &592 & - wmax593 iflag: DO jn = 1, kiter594 IF ( zwl == 0._wp ) THEN ; znummax = zl595 ELSEIF( zwr == 0._wp ) THEN ; znummax = zr596 ELSE597 znummax = ( zr + zl ) / 2.598 zdiv = xkr_zeta + xkr_eta - xkr_eta * znummax599 znum = znummax - 1.600 zws = xkr_wsbio_min * xkr_zeta / zdiv &601 & - ( xkr_wsbio_max * xkr_eta * znum * &602 & xkr_frac**( -xkr_zeta / znum ) / zdiv ) &603 & - wmax604 IF( zws * zwl < 0. ) THEN ; zr = znummax605 ELSE ; zl = znummax606 ENDIF607 zdiv = xkr_zeta + xkr_eta - xkr_eta * zl608 znum = zl - 1.609 zwl = xkr_wsbio_min * xkr_zeta / zdiv &610 & - ( xkr_wsbio_max * xkr_eta * znum * &611 & xkr_frac**( -xkr_zeta / znum ) / zdiv ) &612 & - wmax613 614 zdiv = xkr_zeta + xkr_eta - xkr_eta * zr615 znum = zr - 1.616 zwr = xkr_wsbio_min * xkr_zeta / zdiv &617 & - ( xkr_wsbio_max * xkr_eta * znum * &618 & xkr_frac**( -xkr_zeta / znum ) / zdiv ) &619 & - wmax620 !621 IF ( ABS ( zws ) <= xacc ) EXIT iflag622 !623 ENDIF624 !625 END DO iflag626 627 xnumm(jk) = znummax628 IF (lwp) WRITE(numout,*) ' jk = ', jk, ' wmax = ', wmax,' xnum max = ', xnumm(jk)629 !630 END DO631 !632 ik100 = 10 ! last level where depth less than 100 m633 DO jk = jpkm1, 1, -1634 IF( gdept_1d(jk) > 100. ) ik100 = jk - 1635 END DO636 IF (lwp) WRITE(numout,*)637 IF (lwp) WRITE(numout,*) ' Level corresponding to 100m depth ', ik100 + 1638 IF (lwp) WRITE(numout,*)639 !640 t_oce_co2_exp = 0._wp641 !642 IF( nn_timing == 1 ) CALL timing_stop('p4z_sink_init')643 !644 END SUBROUTINE p4z_sink_init645 646 #endif647 262 648 263 SUBROUTINE p4z_sink2( pwsink, psinkflx, jp_tra, kiter ) … … 769 384 & sinking(jpi,jpj,jpk) , sinking2(jpi,jpj,jpk) , & 770 385 & sinkcal(jpi,jpj,jpk) , sinksil (jpi,jpj,jpk) , & 771 #if defined key_kriest772 & xnumm(jpk) , &773 #else774 386 & sinkfer2(jpi,jpj,jpk) , & 775 #endif776 387 #if defined key_ligand 777 388 & wsfep(jpi,jpj,jpk) , sinkfep(jpi,jpj,jpk) , &
Note: See TracChangeset
for help on using the changeset viewer.
