- 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/P5Z/p5zsink.F90
r6841 r8003 40 40 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: sinkcal, sinksil !: CaCO3 and BSi sinking fluxes 41 41 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: sinkfer !: Small BFe sinking fluxes 42 #if ! defined key_kriest43 42 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: sinkfer2 !: Big iron sinking fluxes 44 #endif45 43 #if defined key_ligand 46 44 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: sinkfep !: Fep sinking fluxes … … 49 47 50 48 INTEGER :: ik100 51 52 #if defined key_kriest53 REAL(wp) :: xkr_sfact !: Sinking factor54 REAL(wp) :: xkr_stick !: Stickiness55 REAL(wp) :: xkr_nnano !: Nbr of cell in nano size class56 REAL(wp) :: xkr_ndiat !: Nbr of cell in diatoms size class57 REAL(wp) :: xkr_nmicro !: Nbr of cell in microzoo size class58 REAL(wp) :: xkr_nmeso !: Nbr of cell in mesozoo size class59 REAL(wp) :: xkr_naggr !: Nbr of cell in aggregates size class60 61 REAL(wp) :: xkr_frac62 63 REAL(wp), PUBLIC :: xkr_dnano !: Size of particles in nano pool64 REAL(wp), PUBLIC :: xkr_ddiat !: Size of particles in diatoms pool65 REAL(wp), PUBLIC :: xkr_dmicro !: Size of particles in microzoo pool66 REAL(wp), PUBLIC :: xkr_dmeso !: Size of particles in mesozoo pool67 REAL(wp), PUBLIC :: xkr_daggr !: Size of particles in aggregates pool68 REAL(wp), PUBLIC :: xkr_wsbio_min !: min vertical particle speed69 REAL(wp), PUBLIC :: xkr_wsbio_max !: max vertical particle speed70 71 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: xnumm !: maximum number of particles in aggregates72 #endif73 49 74 50 !!* Substitution … … 81 57 CONTAINS 82 58 83 #if ! defined key_kriest84 59 !!---------------------------------------------------------------------- 85 60 !! 'standard sinking parameterisation' ??? … … 308 283 END SUBROUTINE p5z_sink_init 309 284 310 #else311 !!----------------------------------------------------------------------312 !! 'Kriest sinking parameterisation' key_kriest ???313 !!----------------------------------------------------------------------314 315 SUBROUTINE p5z_sink ( kt, knt )316 !!---------------------------------------------------------------------317 !! *** ROUTINE p5z_sink ***318 !!319 !! ** Purpose : Compute vertical flux of particulate matter due to320 !! gravitational sinking - Kriest parameterization321 !!322 !! ** Method : - ???323 !!---------------------------------------------------------------------324 !325 INTEGER, INTENT(in) :: kt, knt326 !327 INTEGER :: ji, jj, jk, jit, niter1, niter2328 REAL(wp) :: znum , zeps, zfm, zgm, zsm, zfactn, zfactp329 REAL(wp) :: zdiv , zdiv1, zdiv2, zdiv3, zdiv4, zdiv5330 REAL(wp) :: zval1, zval2, zval3, zval4331 CHARACTER (len=25) :: charout332 REAL(wp), POINTER, DIMENSION(:,:,:) :: znum3d333 REAL(wp), POINTER, DIMENSION(:,:,:) :: zw3d334 REAL(wp), POINTER, DIMENSION(:,: ) :: zw2d335 !!---------------------------------------------------------------------336 !337 IF( nn_timing == 1 ) CALL timing_start('p5z_sink')338 !339 CALL wrk_alloc( jpi, jpj, jpk, znum3d )340 !341 ! Initialisation of variables used to compute Sinking Speed342 ! ---------------------------------------------------------343 344 znum3d(:,:,:) = 0.e0345 zval1 = 1. + xkr_zeta346 zval2 = 1. + xkr_zeta + xkr_eta347 zval3 = 1. + xkr_eta348 349 ! Computation of the vertical sinking speed : Kriest et Evans, 2000350 ! -----------------------------------------------------------------351 352 DO jk = 1, jpkm1353 DO jj = 1, jpj354 DO ji = 1, jpi355 IF( tmask(ji,jj,jk) /= 0.e0 ) THEN356 znum = trb(ji,jj,jk,jppoc) / ( trb(ji,jj,jk,jpnum) + rtrn ) / xkr_massp357 ! -------------- To avoid sinking speed over 50 m/day -------358 znum = MIN( xnumm(jk), znum )359 znum = MAX( 1.1 , znum )360 znum3d(ji,jj,jk) = znum361 !------------------------------------------------------------362 zeps = ( zval1 * znum - 1. )/ ( znum - 1. )363 zfm = xkr_frac**( 1. - zeps )364 zgm = xkr_frac**( zval1 - zeps )365 zdiv = MAX( 1.e-4, ABS( zeps - zval2 ) ) * SIGN( 1., ( zeps - zval2 ) )366 zdiv1 = zeps - zval3367 wsbio3(ji,jj,jk) = xkr_wsbio_min * ( zeps - zval1 ) / zdiv &368 & - xkr_wsbio_max * zgm * xkr_eta / zdiv369 wsbio4(ji,jj,jk) = xkr_wsbio_min * ( zeps-1. ) / zdiv1 &370 & - xkr_wsbio_max * zfm * xkr_eta / zdiv1371 IF( znum == 1.1) wsbio3(ji,jj,jk) = wsbio4(ji,jj,jk)372 ENDIF373 END DO374 END DO375 END DO376 377 wscal(:,:,:) = MAX( wsbio3(:,:,:), 30._wp )378 #if defined key_ligand379 wsfep (:,:,:) = wfep380 #endif381 382 ! INITIALIZE TO ZERO ALL THE SINKING ARRAYS383 ! -----------------------------------------384 385 sinking (:,:,:) = 0.e0386 sinkingn(:,:,:) = 0.e0387 sinkingp(:,:,:) = 0.e0388 sinking2(:,:,:) = 0.e0389 sinkcal (:,:,:) = 0.e0390 sinkfer (:,:,:) = 0.e0391 sinksil (:,:,:) = 0.e0392 #if defined key_ligand393 sinkfep(:,:,:) = 0.e0394 #endif395 396 ! Compute the sedimentation term using p4zsink2 for all the sinking particles397 ! -----------------------------------------------------398 399 niter1 = niter1max400 niter2 = niter2max401 402 DO jit = 1, niter1403 CALL p4z_sink2( wsbio3, sinking , jppoc, niter1 )404 CALL p4z_sink2( wsbio3, sinkingn, jppon, niter1 )405 CALL p4z_sink2( wsbio3, sinkingp, jppop, niter1 )406 CALL p4z_sink2( wsbio3, sinkfer , jpsfe, niter1 )407 CALL p4z_sink2( wscal , sinksil , jpgsi, niter1 )408 CALL p4z_sink2( wscal , sinkcal , jpcal, niter1 )409 #if defined key_ligand410 CALL p4z_sink2( wsfep , sinkfep , jpfep, niter1 )411 #endif412 END DO413 414 DO jit = 1, niter2415 CALL p4z_sink2( wsbio4, sinking2, jpnum, niter2 )416 END DO417 418 ! Total carbon export per year419 IF( iom_use( "tcexp" ) .OR. ( ln_check_mass .AND. kt == nitend .AND. knt == nrdttrc ) ) &420 & t_oce_co2_exp = glob_sum( ( sinking(:,:,ik100) + sinking2(:,:,ik100) ) * e1e2t(:,:) * tmask(:,:,1) )421 !422 IF( lk_iomput ) THEN423 IF( knt == nrdttrc ) THEN424 CALL wrk_alloc( jpi, jpj, zw2d )425 CALL wrk_alloc( jpi, jpj, jpk, zw3d )426 zfact = 1.e+3 * rfact2r ! conversion from mol/l/kt to mol/m3/s427 !428 IF( iom_use( "EPC100" ) ) THEN429 zw2d(:,:) = ( sinking(:,:,ik100) + sinking2(:,:,ik100) ) * zfact * tmask(:,:,1) ! Export of carbon at 100m430 CALL iom_put( "EPC100" , zw2d )431 ENDIF432 IF( iom_use( "EPFE100" ) ) THEN433 zw2d(:,:) = ( sinkfer(:,:,ik100) + sinkfer2(:,:,ik100) ) * zfact * tmask(:,:,1) ! Export of iron at 100m434 CALL iom_put( "EPFE100" , zw2d )435 ENDIF436 IF( iom_use( "EPCAL100" ) ) THEN437 zw2d(:,:) = sinkcal(:,:,ik100) * zfact * tmask(:,:,1) ! Export of calcite at 100m438 CALL iom_put( "EPCAL100" , zw2d )439 ENDIF440 IF( iom_use( "EPSI100" ) ) THEN441 zw2d(:,:) = sinksil(:,:,ik100) * zfact * tmask(:,:,1) ! Export of bigenic silica at 100m442 CALL iom_put( "EPSI100" , zw2d )443 ENDIF444 IF( iom_use( "EXPC" ) ) THEN445 zw3d(:,:,:) = ( sinking(:,:,:) + sinking2(:,:,:) ) * zfact * tmask(:,:,:) ! Export of carbon in the water column446 CALL iom_put( "EXPC" , zw3d )447 ENDIF448 IF( iom_use( "EXPFE" ) ) THEN449 zw3d(:,:,:) = ( sinkfer(:,:,:) + sinkfer2(:,:,:) ) * zfact * tmask(:,:,:) ! Export of iron450 CALL iom_put( "EXPFE" , zw3d )451 ENDIF452 IF( iom_use( "EXPCAL" ) ) THEN453 zw3d(:,:,:) = sinkcal(:,:,:) * zfact * tmask(:,:,:) ! Export of calcite454 CALL iom_put( "EXPCAL" , zw3d )455 ENDIF456 IF( iom_use( "EXPSI" ) ) THEN457 zw3d(:,:,:) = sinksil(:,:,:) * zfact * tmask(:,:,:) ! Export of bigenic silica458 CALL iom_put( "EXPSI" , zw3d )459 ENDIF460 IF( iom_use( "XNUM" ) ) THEN461 zw3d(:,:,:) = znum3d(:,:,:) * tmask(:,:,:) ! Number of particles on aggregats462 CALL iom_put( "XNUM" , zw3d )463 ENDIF464 IF( iom_use( "WSC" ) ) THEN465 zw3d(:,:,:) = wsbio3(:,:,:) * tmask(:,:,:) ! Sinking speed of carbon particles466 CALL iom_put( "WSC" , zw3d )467 ENDIF468 IF( iom_use( "WSN" ) ) THEN469 zw3d(:,:,:) = wsbio4(:,:,:) * tmask(:,:,:) ! Sinking speed of particles number470 CALL iom_put( "WSN" , zw3d )471 ENDIF472 IF( iom_use( "tcexp" ) ) CALL iom_put( "tcexp" , t_oce_co2_exp * zfact ) ! molC/s473 !474 CALL wrk_dealloc( jpi, jpj, zw2d )475 CALL wrk_dealloc( jpi, jpj, jpk, zw3d )476 ENDIF477 ELSE478 IF( ln_diatrc ) THEN479 zfact = 1.e3 * rfact2r480 trc2d(:,:,jp_pcs0_2d + 4) = sinking (:,:,ik100) * zfact * tmask(:,:,1)481 trc2d(:,:,jp_pcs0_2d + 5) = sinking2(:,:,ik100) * zfact * tmask(:,:,1)482 trc2d(:,:,jp_pcs0_2d + 6) = sinkfer (:,:,ik100) * zfact * tmask(:,:,1)483 trc2d(:,:,jp_pcs0_2d + 7) = sinkfer2(:,:,ik100) * zfact * tmask(:,:,1)484 trc2d(:,:,jp_pcs0_2d + 8) = sinksil (:,:,ik100) * zfact * tmask(:,:,1)485 trc2d(:,:,jp_pcs0_2d + 9) = sinkcal (:,:,ik100) * zfact * tmask(:,:,1)486 ENDIF487 ENDIF488 !489 IF(ln_ctl) THEN ! print mean trends (used for debugging)490 WRITE(charout, FMT="('sink')")491 CALL prt_ctl_trc_info(charout)492 CALL prt_ctl_trc(tab4d=tra, mask=tmask, clinfo=ctrcnm)493 ENDIF494 !495 CALL wrk_dealloc( jpi, jpj, jpk, znum3d )496 !497 IF( nn_timing == 1 ) CALL timing_stop('p5z_sink')498 !499 END SUBROUTINE p5z_sink500 501 502 SUBROUTINE p5z_sink_init503 !!----------------------------------------------------------------------504 !! *** ROUTINE p5z_sink_init ***505 !!506 !! ** Purpose : Initialization of sinking parameters507 !! Kriest parameterization only508 !!509 !! ** Method : Read the nampiskrs namelist and check the parameters510 !! called at the first timestep511 !!512 !! ** input : Namelist nampiskrs513 !!----------------------------------------------------------------------514 INTEGER :: jk, jn, kiter515 INTEGER :: ios ! Local integer output status for namelist read516 REAL(wp) :: znum, zdiv517 REAL(wp) :: zws, zwr, zwl,wmax, znummax518 REAL(wp) :: zmin, zmax, zl, zr, xacc519 !520 NAMELIST/nampiskrs/ xkr_sfact, xkr_stick , &521 & xkr_nnano, xkr_ndiat, xkr_nmicro, xkr_nmeso, xkr_naggr522 !!----------------------------------------------------------------------523 !524 IF( nn_timing == 1 ) CALL timing_start('p5z_sink_init')525 !526 527 REWIND( numnatp_ref ) ! Namelist nampiskrs in reference namelist : Pisces sinking Kriest528 READ ( numnatp_ref, nampiskrs, IOSTAT = ios, ERR = 901)529 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nampiskrs in reference namelist', lwp )530 531 REWIND( numnatp_cfg ) ! Namelist nampiskrs in configuration namelist : Pisces sinking Kriest532 READ ( numnatp_cfg, nampiskrs, IOSTAT = ios, ERR = 902 )533 902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nampiskrs in configuration namelist', lwp )534 IF(lwm) WRITE ( numonp, nampiskrs )535 536 IF(lwp) THEN537 WRITE(numout,*)538 WRITE(numout,*) ' Namelist : nampiskrs'539 WRITE(numout,*) ' Sinking factor xkr_sfact = ', xkr_sfact540 WRITE(numout,*) ' Stickiness xkr_stick = ', xkr_stick541 WRITE(numout,*) ' Nbr of cell in nano size class xkr_nnano = ', xkr_nnano542 WRITE(numout,*) ' Nbr of cell in diatoms size class xkr_ndiat = ', xkr_ndiat543 WRITE(numout,*) ' Nbr of cell in microzoo size class xkr_nmicro = ', xkr_nmicro544 WRITE(numout,*) ' Nbr of cell in mesozoo size class xkr_nmeso = ', xkr_nmeso545 WRITE(numout,*) ' Nbr of cell in aggregates size class xkr_naggr = ', xkr_naggr546 ENDIF547 548 549 ! max and min vertical particle speed550 xkr_wsbio_min = xkr_sfact * xkr_mass_min**xkr_eta551 xkr_wsbio_max = xkr_sfact * xkr_mass_max**xkr_eta552 IF (lwp) WRITE(numout,*) ' max and min vertical particle speed ', xkr_wsbio_min, xkr_wsbio_max553 554 !555 ! effect of the sizes of the different living pools on particle numbers556 ! nano = 2um-20um -> mean size=6.32 um -> ws=2.596 -> xnum=xnnano=2.337557 ! diat and microzoo = 10um-200um -> 44.7 -> 8.732 -> xnum=xndiat=3.718558 ! mesozoo = 200um-2mm -> 632.45 -> 45.14 -> xnum=xnmeso=7.147559 ! aggregates = 200um-10mm -> 1414 -> 74.34 -> xnum=xnaggr=9.877560 ! doc aggregates = 1um561 ! ----------------------------------------------------------562 563 xkr_dnano = 1. / ( xkr_massp * xkr_nnano )564 xkr_ddiat = 1. / ( xkr_massp * xkr_ndiat )565 xkr_dmicro = 1. / ( xkr_massp * xkr_nmicro )566 xkr_dmeso = 1. / ( xkr_massp * xkr_nmeso )567 xkr_daggr = 1. / ( xkr_massp * xkr_naggr )568 569 !!---------------------------------------------------------------------570 !! 'key_kriest' ???571 !!---------------------------------------------------------------------572 ! COMPUTATION OF THE VERTICAL PROFILE OF MAXIMUM SINKING SPEED573 ! Search of the maximum number of particles in aggregates for each k-level.574 ! Bissection Method575 !--------------------------------------------------------------------576 IF (lwp) THEN577 WRITE(numout,*)578 WRITE(numout,*)' kriest : Compute maximum number of particles in aggregates'579 ENDIF580 581 xacc = 0.001_wp582 kiter = 50583 zmin = 1.10_wp584 zmax = xkr_mass_max / xkr_mass_min585 xkr_frac = zmax586 587 DO jk = 1,jpk588 zl = zmin589 zr = zmax590 wmax = 0.5 * fse3t(1,1,jk) * rday * float(niter1max) / rfact2591 zdiv = xkr_zeta + xkr_eta - xkr_eta * zl592 znum = zl - 1.593 zwl = xkr_wsbio_min * xkr_zeta / zdiv &594 & - ( xkr_wsbio_max * xkr_eta * znum * &595 & xkr_frac**( -xkr_zeta / znum ) / zdiv ) &596 & - wmax597 598 zdiv = xkr_zeta + xkr_eta - xkr_eta * zr599 znum = zr - 1.600 zwr = xkr_wsbio_min * xkr_zeta / zdiv &601 & - ( xkr_wsbio_max * xkr_eta * znum * &602 & xkr_frac**( -xkr_zeta / znum ) / zdiv ) &603 & - wmax604 iflag: DO jn = 1, kiter605 IF ( zwl == 0._wp ) THEN ; znummax = zl606 ELSEIF( zwr == 0._wp ) THEN ; znummax = zr607 ELSE608 znummax = ( zr + zl ) / 2.609 zdiv = xkr_zeta + xkr_eta - xkr_eta * znummax610 znum = znummax - 1.611 zws = xkr_wsbio_min * xkr_zeta / zdiv &612 & - ( xkr_wsbio_max * xkr_eta * znum * &613 & xkr_frac**( -xkr_zeta / znum ) / zdiv ) &614 & - wmax615 IF( zws * zwl < 0. ) THEN ; zr = znummax616 ELSE ; zl = znummax617 ENDIF618 zdiv = xkr_zeta + xkr_eta - xkr_eta * zl619 znum = zl - 1.620 zwl = xkr_wsbio_min * xkr_zeta / zdiv &621 & - ( xkr_wsbio_max * xkr_eta * znum * &622 & xkr_frac**( -xkr_zeta / znum ) / zdiv ) &623 & - wmax624 625 zdiv = xkr_zeta + xkr_eta - xkr_eta * zr626 znum = zr - 1.627 zwr = xkr_wsbio_min * xkr_zeta / zdiv &628 & - ( xkr_wsbio_max * xkr_eta * znum * &629 & xkr_frac**( -xkr_zeta / znum ) / zdiv ) &630 & - wmax631 !632 IF ( ABS ( zws ) <= xacc ) EXIT iflag633 !634 ENDIF635 !636 END DO iflag637 638 xnumm(jk) = znummax639 IF (lwp) WRITE(numout,*) ' jk = ', jk, ' wmax = ', wmax,' xnum max = ', xnumm(jk)640 !641 END DO642 !643 ik100 = 10 ! last level where depth less than 100 m644 DO jk = jpkm1, 1, -1645 IF( gdept_1d(jk) > 100. ) ik100 = jk - 1646 END DO647 IF (lwp) WRITE(numout,*)648 IF (lwp) WRITE(numout,*) ' Level corresponding to 100m depth ', ik100 + 1649 IF (lwp) WRITE(numout,*)650 !651 t_oce_co2_exp = 0._wp652 !653 IF( nn_timing == 1 ) CALL timing_stop('p5z_sink_init')654 !655 END SUBROUTINE p5z_sink_init656 657 #endif658 659 285 SUBROUTINE p4z_sink2( pwsink, psinkflx, jp_tra, kiter ) 660 286 !!--------------------------------------------------------------------- … … 782 408 783 409 & sinkcal(jpi,jpj,jpk) , sinksil (jpi,jpj,jpk) , & 784 #if defined key_kriest785 & xnumm(jpk) , &786 #else787 410 & sinkfer2(jpi,jpj,jpk) , & 788 #endif789 411 #if defined key_ligand 790 412 & wsfep(jpi,jpj,jpk) , sinkfep(jpi,jpj,jpk) , &
Note: See TracChangeset
for help on using the changeset viewer.
