Skip to content
Snippets Groups Projects

Compare revisions

Changes are shown as if the source revision was being merged into the target revision. Learn more about comparing revisions.

Source

Select target project
No results found

Target

Select target project
  • nemo/nemo
  • sparonuz/nemo
  • hatfield/nemo
  • extdevs/nemo
4 results
Show changes
Showing
with 769 additions and 387 deletions
......@@ -65,7 +65,7 @@ CONTAINS
!!
INTEGER :: ji, jj, jk, jl, ikt
REAL(wp) :: zgeolpoc, zfact, zwork, ze3t, zsedpocd, zmaskt
REAL(wp), DIMENSION(jpi,jpj) :: zsedpoca
REAL(wp), DIMENSION(A2D(0)) :: zsedpoca
CHARACTER (len=25) :: charout
!!---------------------------------------------------------------------
!
......@@ -106,7 +106,7 @@ CONTAINS
tr(ji,jj,1,jpno3,Krhs) = tr(ji,jj,1,jpno3,Krhs) + zgeolpoc * cmask(ji,jj) / areacot / e3t(ji,jj,1,Kmm)
END_2D
CALL lbc_lnk( 'p2zexp', sedpocn, 'T', 1.0_wp )
! CALL lbc_lnk( 'p2zexp', sedpocn, 'T', 1.0_wp )
! Oa & Ek: diagnostics depending on jpdia2d ! left as example
IF( lk_iomput ) CALL iom_put( "SEDPOC" , sedpocn )
......@@ -120,7 +120,7 @@ CONTAINS
!
ELSE
!
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
DO_2D( 0, 0, 0, 0 )
zsedpocd = zsedpoca(ji,jj) - 2. * sedpocn(ji,jj) + sedpocb(ji,jj) ! time laplacian on tracers
sedpocb(ji,jj) = sedpocn(ji,jj) + rn_atfp * zsedpocd ! sedpocb <-- filtered sedpocn
sedpocn(ji,jj) = zsedpoca(ji,jj) ! sedpocn <-- sedpoca
......@@ -156,8 +156,8 @@ CONTAINS
INTEGER, INTENT(in) :: Kmm ! time level index
INTEGER :: ji, jj, jk
REAL(wp) :: zmaskt, zfluo, zfluu
REAL(wp), DIMENSION(jpi,jpj ) :: zrro
REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdm0
REAL(wp), DIMENSION(A2D(0) ) :: zrro, zarea
REAL(wp), DIMENSION(A2D(0),jpk) :: zdm0
!!---------------------------------------------------------------------
!
IF(lwp) THEN
......@@ -171,9 +171,9 @@ CONTAINS
! Calculate vertical distribution of newly formed biogenic poc
! in the water column in the case of max. possible bottom depth
! ------------------------------------------------------------
zdm0 = 0._wp
zrro = 1._wp
DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, jpkb, jpkm1 )
zdm0(:,:,:) = 0._wp
zrro(:,:) = 1._wp
DO_3D( 0, 0, 0, 0, jpkb, jpkm1 )
zfluo = ( gdepw(ji,jj,jk ,Kmm) / gdepw(ji,jj,jpkb,Kmm) )**xhr
zfluu = ( gdepw(ji,jj,jk+1,Kmm) / gdepw(ji,jj,jpkb,Kmm) )**xhr
IF( zfluo.GT.1. ) zfluo = 1._wp
......@@ -190,14 +190,14 @@ CONTAINS
! ----------------------------------------------------------------------
dminl(:,:) = 0._wp
dmin3(:,:,:) = zdm0
DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpk )
DO_3D( 0, 0, 0, 0, 1, jpk )
IF( tmask(ji,jj,jk) == 0._wp ) THEN
dminl(ji,jj) = dminl(ji,jj) + dmin3(ji,jj,jk)
dmin3(ji,jj,jk) = 0._wp
ENDIF
END_3D
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
DO_2D( 0, 0, 0, 0 )
IF( tmask(ji,jj,1) == 0 ) dmin3(ji,jj,1) = 0._wp
END_2D
......@@ -209,8 +209,11 @@ CONTAINS
IF( zmaskt == 0. ) cmask(ji,jj) = 1._wp
END IF
END_2D
CALL lbc_lnk( 'p2zexp', cmask , 'T', 1.0_wp ) ! lateral boundary conditions on cmask (sign unchanged)
areacot = glob_sum( 'p2zexp', e1e2t(:,:) * cmask(:,:) )
! CALL lbc_lnk( 'p2zexp', cmask , 'T', 1.0_wp ) ! lateral boundary conditions on cmask (sign unchanged)
DO_2D( 0, 0, 0, 0 )
zarea(ji,jj) = e1e2t(ji,jj) * cmask(ji,jj)
END_2D
areacot = glob_sum( 'p2zexp', zarea(:,:) )
!
IF( ln_rsttr ) THEN
CALL iom_get( numrtr, jpdom_auto, 'SEDB'//ctrcnm(jpdet), sedpocb(:,:) )
......@@ -226,8 +229,8 @@ CONTAINS
!!----------------------------------------------------------------------
!! *** ROUTINE p2z_exp_alloc ***
!!----------------------------------------------------------------------
ALLOCATE( cmask(jpi,jpj) , dminl(jpi,jpj) , dmin3(jpi,jpj,jpk), &
& sedpocb(jpi,jpj) , sedpocn(jpi,jpj), STAT=p2z_exp_alloc )
ALLOCATE( cmask(A2D(0)) , dminl(A2D(0)) , dmin3(A2D(0),jpk), &
& sedpocb(A2D(0)) , sedpocn(A2D(0)), STAT=p2z_exp_alloc )
IF( p2z_exp_alloc /= 0 ) CALL ctl_stop( 'STOP', 'p2z_exp_alloc : failed to allocate arrays.' )
!
END FUNCTION p2z_exp_alloc
......
......@@ -70,8 +70,8 @@ CONTAINS
REAL(wp) :: zpig ! log of the total pigment
REAL(wp) :: zkr, zkg ! total absorption coefficient in red and green
REAL(wp) :: zcoef ! temporary scalar
REAL(wp), DIMENSION(jpi,jpj ) :: zpar100, zpar0m
REAL(wp), DIMENSION(jpi,jpj,jpk) :: zparr, zparg
REAL(wp), DIMENSION(A2D(0) ) :: zpar100, zpar0m
REAL(wp), DIMENSION(A2D(0),jpk) :: zparr, zparg
!!---------------------------------------------------------------------
!
IF( ln_timing ) CALL timing_start('p2z_opt')
......@@ -85,8 +85,14 @@ CONTAINS
! ! surface irradiance
! ! ------------------
IF( ln_dm2dc ) THEN ; zpar0m(:,:) = qsr_mean(:,:) * 0.43
ELSE ; zpar0m(:,:) = qsr (:,:) * 0.43
IF( ln_dm2dc ) THEN
DO_2D( 0, 0, 0, 0 )
zpar0m(ji,jj) = qsr_mean(ji,jj) * 0.43
END_2D
ELSE
DO_2D( 0, 0, 0, 0 )
zpar0m(ji,jj) = qsr(ji,jj) * 0.43
END_2D
ENDIF
zpar100(:,:) = zpar0m(:,:) * 0.01
zparr (:,:,1) = zpar0m(:,:) * 0.5
......@@ -94,14 +100,14 @@ CONTAINS
! ! Photosynthetically Available Radiation (PAR)
zcoef = 12 * redf / rcchl / rpig ! --------------------------------------
DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 2, jpk )
DO_3D( 0, 0, 0, 0, 2, jpk )
zpig = LOG( MAX( TINY(0.), tr(ji,jj,jk-1,jpphy,Kmm) ) * zcoef )
zkr = xkr0 + xkrp * EXP( xlr * zpig )
zkg = xkg0 + xkgp * EXP( xlg * zpig )
zparr(ji,jj,jk) = zparr(ji,jj,jk-1) * EXP( -zkr * e3t(ji,jj,jk-1,Kmm) )
zparg(ji,jj,jk) = zparg(ji,jj,jk-1) * EXP( -zkg * e3t(ji,jj,jk-1,Kmm) )
END_3D
DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1 ) ! mean par at t-levels
DO_3D( 0, 0, 0, 0, 1, jpkm1 ) ! mean par at t-levels
zpig = LOG( MAX( TINY(0.), tr(ji,jj,jk,jpphy,Kmm) ) * zcoef )
zkr = xkr0 + xkrp * EXP( xlr * zpig )
zkg = xkg0 + xkgp * EXP( xlg * zpig )
......@@ -113,11 +119,11 @@ CONTAINS
! ! Euphotic layer
! ! --------------
neln(:,:) = 1 ! euphotic layer level
DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1 ) ! (i.e. 1rst T-level strictly below EL bottom)
DO_3D( 0, 0, 0, 0, 1, jpkm1 ) ! (i.e. 1rst T-level strictly below EL bottom)
IF( etot(ji,jj,jk) >= zpar100(ji,jj) ) neln(ji,jj) = jk + 1
END_3D
! ! Euphotic layer depth
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
DO_2D( 0, 0, 0, 0 )
heup(ji,jj) = gdepw(ji,jj,neln(ji,jj),Kmm)
END_2D
......
......@@ -61,10 +61,11 @@ CONTAINS
INTEGER, INTENT( in ) :: kt ! ocean time-step index
INTEGER, INTENT( in ) :: Kmm, Krhs ! time level indices
!
INTEGER :: ji, jj, jk, jl, ierr
INTEGER :: ji, jj, jk
CHARACTER (len=25) :: charout
REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zw2d
REAL(wp), DIMENSION(jpi,jpj,jpk) :: zwork, ztra
REAL(wp), DIMENSION(A2D(0),jpk) :: zwork
REAL(wp) :: ztra
!!---------------------------------------------------------------------
!
IF( ln_timing ) CALL timing_start('p2z_sed')
......@@ -83,22 +84,26 @@ CONTAINS
zwork(:,:,jpk) = 0.e0 ! bottom value set to zero
! tracer flux at w-point: we use -vsed (downward flux) with simplification : no e1*e2
DO jk = 2, jpkm1
zwork(:,:,jk) = -vsed * tr(:,:,jk-1,jpdet,Kmm)
END DO
DO_3D( 0, 0, 0, 0, 2, jpkm1 )
zwork(ji,jj,jk) = -vsed * tr(ji,jj,jk-1,jpdet,Kmm)
END_3D
! tracer flux divergence at t-point added to the general trend
DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1 )
ztra(ji,jj,jk) = - ( zwork(ji,jj,jk) - zwork(ji,jj,jk+1) ) / e3t(ji,jj,jk,Kmm)
tr(ji,jj,jk,jpdet,Krhs) = tr(ji,jj,jk,jpdet,Krhs) + ztra(ji,jj,jk)
DO_3D( 0, 0, 0, 0, 1, jpkm1 )
ztra = - ( zwork(ji,jj,jk) - zwork(ji,jj,jk+1) ) / e3t(ji,jj,jk,Kmm)
tr(ji,jj,jk,jpdet,Krhs) = tr(ji,jj,jk,jpdet,Krhs) + ztra
END_3D
IF( lk_iomput ) THEN
IF( iom_use( "TDETSED" ) ) THEN
ALLOCATE( zw2d(jpi,jpj) )
zw2d(:,:) = ztra(:,:,1) * e3t(:,:,1,Kmm) * 86400._wp
ALLOCATE( zw2d(A2D(0)) )
DO_2D( 0, 0, 0, 0 )
zw2d(ji,jj) = - ( zwork(ji,jj,1) - zwork(ji,jj,2) ) * 86400._wp
END_2D
DO jk = 2, jpkm1
zw2d(:,:) = zw2d(:,:) + ztra(:,:,jk) * e3t(:,:,jk,Kmm) * 86400._wp
DO_2D( 0, 0, 0, 0 )
zw2d(ji,jj) = zw2d(ji,jj) - ( zwork(ji,jj,jk) - zwork(ji,jj,jk+1) ) * 86400._wp
END_2D
END DO
CALL iom_put( "TDETSED", zw2d )
DEALLOCATE( zw2d )
......
......@@ -70,7 +70,7 @@ CONTAINS
! PISCES part
IF( ln_p4z ) THEN
!
DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1 )
DO_3D( 0, 0, 0, 0, 1, jpkm1 )
!
zfact = xstep * xdiss(ji,jj,jk)
! Part I : Coagulation dependent on turbulence
......@@ -117,7 +117,7 @@ CONTAINS
ELSE ! ln_p5z
! PISCES-QUOTA part
!
DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1 )
DO_3D( 0, 0, 0, 0, 1, jpkm1 )
!
zfact = xstep * xdiss(ji,jj,jk)
! Part I : Coagulation dependent on turbulence
......
......@@ -13,6 +13,7 @@ MODULE p4zbc
USE sms_pisces ! PISCES Source Minus Sink variables
USE iom ! I/O manager
USE fldread ! time interpolation
USE prtctl ! print control for debugging
USE trcbc
IMPLICIT NONE
......@@ -36,7 +37,6 @@ MODULE p4zbc
LOGICAL , PUBLIC :: ll_river !: boolean for river input of nutrients
LOGICAL , PUBLIC :: ll_ndepo !: boolean for atmospheric deposition of N
TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: sf_dust ! structure of input dust
TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: sf_ironsed ! structure of input iron from sediment
TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: sf_hydrofe ! structure of input iron from sediment
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: dust !: dust fields
......@@ -74,6 +74,8 @@ CONTAINS
REAL(wp) :: zcoef, zwdust, zrivdin, zdustdep, zndep
!
CHARACTER (len=25) :: charout
REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: zw3d
REAL(wp), ALLOCATABLE, DIMENSION(:,: ) :: zw2d
!!---------------------------------------------------------------------
!
IF( ln_timing ) CALL timing_start('p4z_bc')
......@@ -84,7 +86,9 @@ CONTAINS
IF( ll_dust ) THEN
!
CALL fld_read( kt, 1, sf_dust )
dust(:,:) = MAX( rtrn, sf_dust(1)%fnow(:,:,1) )
DO_2D( 0, 0, 0, 0 )
dust(ji,jj) = MAX( rtrn, sf_dust(1)%fnow(ji,jj,1) )
END_2D
!
! Iron solubilization of particles in the water column
! dust in kg/m2/s ---> 1/55.85 to put in mol/Fe ; wdust in m/d
......@@ -99,7 +103,7 @@ CONTAINS
! Atmospheric input of Iron dissolves in the water column
IF ( ln_trc_sbc(jpfer) ) THEN
DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 2, jpkm1 )
DO_3D( 0, 0, 0, 0, 2, jpkm1 )
zdustdep = dust(ji,jj) * zwdust * rfact * EXP( -gdept(ji,jj,jk,Kmm) /( 250. * wdust ) )
!
tr(ji,jj,jk,jpfer,Krhs) = tr(ji,jj,jk,jpfer,Krhs) + zdustdep * mfrac / mMass_Fe
......@@ -107,16 +111,18 @@ CONTAINS
IF( lk_iomput ) THEN
! surface downward dust depo of iron
ALLOCATE( zw2d(A2D(0)) )
jl = n_trc_indsbc(jpfer)
CALL iom_put( "Irondep", ( rf_trsfac(jl) * sf_trcsbc(jl)%fnow(:,:,1) / rn_sbc_time ) * 1.e+3 * tmask(:,:,1) )
zw2d(A2D(0)) = rf_trsfac(jl) * ( sf_trcsbc(jl)%fnow(A2D(0),1) / rn_sbc_time ) * 1.e+3 * tmask(A2D(0),1)
CALL iom_put( "Irondep", zw2d )
DEALLOCATE( zw2d )
ENDIF
ENDIF
! Atmospheric input of PO4 dissolves in the water column
IF ( ln_trc_sbc(jppo4) ) THEN
DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 2, jpkm1 )
DO_3D( 0, 0, 0, 0, 2, jpkm1 )
zdustdep = dust(ji,jj) * zwdust * rfact * EXP( -gdept(ji,jj,jk,Kmm) /( 250. * wdust ) )
!
tr(ji,jj,jk,jppo4,Krhs) = tr(ji,jj,jk,jppo4,Krhs) + zdustdep * 1.e-3 / mMass_P
......@@ -125,7 +131,7 @@ CONTAINS
! Atmospheric input of Si dissolves in the water column
IF ( ln_trc_sbc(jpsil) ) THEN
DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 2, jpkm1 )
DO_3D( 0, 0, 0, 0, 2, jpkm1 )
zdustdep = dust(ji,jj) * zwdust * rfact * EXP( -gdept(ji,jj,jk,Kmm) /( 250. * wdust ) )
!
tr(ji,jj,jk,jpsil,Krhs) = tr(ji,jj,jk,jpsil,Krhs) + zdustdep * 0.269 / mMass_Si
......@@ -135,7 +141,10 @@ CONTAINS
!
IF( lk_iomput ) THEN
! dust concentration at surface
CALL iom_put( "pdust" , dust(:,:) / ( wdust / rday ) * tmask(:,:,1) ) ! dust concentration at surface
ALLOCATE( zw2d(A2D(0)) )
zw2d(A2D(0)) = dust(A2D(0)) / ( wdust / rday ) * tmask(A2D(0),1)
CALL iom_put( "pdust", zw2d )
DEALLOCATE( zw2d )
ENDIF
ENDIF
......@@ -144,7 +153,7 @@ CONTAINS
! ----------------------------------------------------------
IF( ll_river ) THEN
jl = n_trc_indcbc(jpno3)
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
DO_2D( 0, 0, 0, 0 )
DO jk = 1, nk_rnf(ji,jj)
zcoef = rn_rfact / ( e1e2t(ji,jj) * h_rnf(ji,jj) * rn_cbc_time ) * tmask(ji,jj,1)
zrivdin = rf_trcfac(jl) * sf_trccbc(jl)%fnow(ji,jj,1) * zcoef
......@@ -158,14 +167,14 @@ CONTAINS
IF( ll_ndepo ) THEN
IF( ln_trc_sbc(jpno3) ) THEN
jl = n_trc_indsbc(jpno3)
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
DO_2D( 0, 0, 0, 0 )
zndep = rf_trsfac(jl) * sf_trcsbc(jl)%fnow(ji,jj,1) / e3t(ji,jj,1,Kmm) / rn_sbc_time
tr(ji,jj,1,jptal,Krhs) = tr(ji,jj,1,jptal,Krhs) - rno3 * zndep * rfact
END_2D
ENDIF
IF( ln_trc_sbc(jpnh4) ) THEN
jl = n_trc_indsbc(jpnh4)
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
DO_2D( 0, 0, 0, 0 )
zndep = rf_trsfac(jl) * sf_trcsbc(jl)%fnow(ji,jj,1) / e3t(ji,jj,1,Kmm) / rn_sbc_time
tr(ji,jj,1,jptal,Krhs) = tr(ji,jj,1,jptal,Krhs) + rno3 * zndep * rfact
END_2D
......@@ -183,41 +192,71 @@ CONTAINS
! Simple parameterization assuming a fixed constant concentration in
! sea-ice (icefeinput)
! ------------------------------------------------------------------
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
DO_2D( 0, 0, 0, 0 )
zdep = rfact / e3t(ji,jj,1,Kmm)
zwflux = fmmflx(ji,jj) / 1000._wp
zironice = MAX( -0.99 * tr(ji,jj,1,jpfer,Kbb), -zwflux * icefeinput * zdep )
tr(ji,jj,1,jpfer,Krhs) = tr(ji,jj,1,jpfer,Krhs) + zironice
END_2D
!
! iron flux from ice
IF( lk_iomput ) &
& CALL iom_put( "Ironice", MAX( -0.99 * tr(:,:,1,jpfer,Kbb), (-1.*fmmflx(:,:)/1000._wp )*icefeinput*1.e+3*tmask(:,:,1)) )
IF( lk_iomput ) THEN
! iron flux from ice
ALLOCATE( zw2d(A2D(0)) )
zw2d(A2D(0)) = MAX( -0.99 * tr(A2D(0),1,jpfer,Kbb), (-1.*fmmflx(A2D(0))/1000._wp )*icefeinput*1.e+3*tmask(A2D(0),1))
CALL iom_put( "Ironice", zw2d )
DEALLOCATE( zw2d )
ENDIF
!
ENDIF
! Add the external input of iron from sediment mobilization
! ------------------------------------------------------
IF( ln_ironsed .AND. .NOT.lk_sed ) THEN
tr(:,:,:,jpfer,Krhs) = tr(:,:,:,jpfer,Krhs) + ironsed(:,:,:) * rfact
!
IF( lk_iomput ) CALL iom_put( "Ironsed", ironsed(:,:,:) * 1.e+3 * tmask(:,:,:) )
DO_3D( 0, 0, 0, 0, 1, jpkm1 )
tr(ji,jj,jk,jpfer,Krhs) = tr(ji,jj,jk,jpfer,Krhs) + ironsed(ji,jj,jk) * rfact
END_3D
!
IF( lk_iomput ) THEN
ALLOCATE( zw3d(A2D(0),jpk) ) ; zw3d(A2D(0),jpk) = 0._wp
zw3d(A2D(0),1:jpkm1) = ironsed(A2D(0),1:jpkm1) * 1.e+3 * tmask(A2D(0),1:jpkm1)
CALL iom_put( "Ironsed", zw3d )
DEALLOCATE( zw3d )
ENDIF
ENDIF
! Add the external input of iron from hydrothermal vents
! ------------------------------------------------------
IF( ln_hydrofe ) THEN
CALL fld_read( kt, 1, sf_hydrofe )
DO jk = 1, jpk
hydrofe(:,:,jk) = ( MAX( rtrn, sf_hydrofe(1)%fnow(:,:,jk) ) * hratio ) &
& / ( e1e2t(:,:) * e3t(:,:,jk,Kmm) * ryyss + rtrn ) / 1000._wp &
& * tmask(:,:,jk)
ENDDO
tr(:,:,:,jpfer,Krhs) = tr(:,:,:,jpfer,Krhs) + hydrofe(:,:,:) * rfact
IF( ln_ligand ) tr(:,:,:,jplgw,Krhs) = tr(:,:,:,jplgw,Krhs) + ( hydrofe(:,:,:) * lgw_rath ) * rfact
DO_3D( 0, 0, 0, 0, 1, jpkm1 )
hydrofe(ji,jj,jk) = ( MAX( rtrn, sf_hydrofe(1)%fnow(ji,jj,jk) ) * hratio ) &
& / ( e1e2t(ji,jj) * e3t(ji,jj,jk,Kmm) * ryyss + rtrn ) / 1000._wp &
& * tmask(ji,jj,jk)
END_3D
DO_3D( 0, 0, 0, 0, 1, jpkm1 )
tr(ji,jj,jk,jpfer,Krhs) = tr(ji,jj,jk,jpfer,Krhs) + hydrofe(ji,jj,jk) * rfact
END_3D
IF( ln_ligand ) THEN
DO_3D( 0, 0, 0, 0, 1, jpkm1 )
tr(ji,jj,jk,jplgw,Krhs) = tr(ji,jj,jk,jplgw,Krhs) + ( hydrofe(ji,jj,jk) * lgw_rath ) * rfact
END_3D
ENDIF
!
IF( lk_iomput ) CALL iom_put( "HYDR", hydrofe(:,:,:) * 1.e+3 * tmask(:,:,:) ) ! hydrothermal iron input
IF( lk_iomput ) THEN
! hydrothermal iron input
ALLOCATE( zw3d(A2D(0),jpk) ) ; zw3d(A2D(0),jpk) = 0._wp
zw3d(A2D(0),1:jpkm1) = hydrofe(A2D(0),1:jpkm1) * 1.e+3 * tmask(A2D(0),1:jpkm1)
CALL iom_put( "HYDR", zw3d )
DEALLOCATE( zw3d )
ENDIF
ENDIF
!
IF(sn_cfctl%l_prttrc) THEN ! print mean trends (used for debugging)
WRITE(charout, FMT="('bc')")
CALL prt_ctl_info( charout, cdcomp = 'top' )
CALL prt_ctl(tab4d_1=tr(:,:,:,:,Krhs), mask1=tmask, clinfo=ctrcnm)
ENDIF
!
IF( ln_timing ) CALL timing_stop('p4z_bc')
!
END SUBROUTINE p4z_bc
......@@ -303,7 +342,7 @@ CONTAINS
IF(lwp) WRITE(numout,*) ' initialize dust input from atmosphere '
IF(lwp) WRITE(numout,*) ' ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ '
!
ALLOCATE( dust(jpi,jpj) )
ALLOCATE( dust(A2D(0)) )
!
ALLOCATE( sf_dust(1), STAT=ierr ) !* allocate and fill sf_sst (forcing structure) with sn_sst
IF( ierr > 0 ) CALL ctl_stop( 'STOP', 'p4z_bc_init: unable to allocate sf_dust structure' )
......@@ -321,7 +360,7 @@ CONTAINS
IF(lwp) WRITE(numout,*)
IF(lwp) WRITE(numout,*) ' ==>>> ln_ironsed=T , computation of an island mask to enhance coastal supply of iron'
!
ALLOCATE( ironsed(jpi,jpj,jpk) ) ! allocation
ALLOCATE( ironsed(A2D(0),jpk) ) ! allocation
!
CALL iom_open ( TRIM( sn_ironsed%clname ), numiron )
ALLOCATE( zcmask(jpi,jpj,jpk) )
......@@ -350,7 +389,7 @@ CONTAINS
!
CALL lbc_lnk( 'p4zbc', zcmask , 'T', 1.0_wp ) ! lateral boundary conditions on cmask (sign unchanged)
!
DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpk )
DO_3D( 0, 0, 0, 0, 1, jpkm1 )
zexpide = MIN( 8.,( gdept(ji,jj,jk,Kmm) / 500. )**(-1.5) )
zdenitide = -0.9543 + 0.7662 * LOG( zexpide ) - 0.235 * LOG( zexpide )**2
zcmask(ji,jj,jk) = zcmask(ji,jj,jk) * MIN( 1., EXP( zdenitide ) / 0.5 )
......@@ -358,9 +397,9 @@ CONTAINS
! Coastal supply of iron
! -------------------------
ironsed(:,:,jpk) = 0._wp
DO jk = 1, jpkm1
ironsed(:,:,jk) = sedfeinput * zcmask(:,:,jk) / ( e3t_0(:,:,jk) * rday )
END DO
DO_3D( 0, 0, 0, 0, 1, jpkm1 )
ironsed(ji,jj,jk) = sedfeinput * zcmask(ji,jj,jk) / ( e3t_0(ji,jj,jk) * rday )
END_3D
DEALLOCATE( zcmask)
ENDIF
!
......@@ -371,7 +410,7 @@ CONTAINS
IF(lwp) WRITE(numout,*)
IF(lwp) WRITE(numout,*) ' ==>>> ln_hydrofe=T , Input of iron from hydrothermal vents'
!
ALLOCATE( hydrofe(jpi,jpj,jpk) ) ! allocation
ALLOCATE( hydrofe(A2D(0),jpk) ) ! allocation
!
ALLOCATE( sf_hydrofe(1), STAT=ierr ) !* allocate and fill sf_sst (forcing structure) with sn_sst
IF( ierr > 0 ) CALL ctl_stop( 'STOP', 'p4z_bc_init: unable to allocate sf_hydro structure' )
......
......@@ -72,7 +72,7 @@ CONTAINS
! OF PHYTOPLANKTON AND DETRITUS. Shear rate is supposed to equal 1
! in the mixed layer and 0.1 below the mixed layer.
xdiss(:,:,:) = 1.
DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 2, jpkm1 )
DO_3D( 0, 0, 0, 0, 2, jpkm1 )
IF( gdepw(ji,jj,jk+1,Kmm) > hmld(ji,jj) ) xdiss(ji,jj,jk) = 0.01
END_3D
......
......@@ -168,9 +168,13 @@ CONTAINS
! practical salinity
! -------------------------------------------------------------
IF (neos == -1) THEN
salinprac(:,:,:) = ts(:,:,:,jp_sal,Kmm) * 35.0 / 35.16504
DO_3D( 0, 0, 0, 0, 1, jpk )
salinprac(ji,jj,jk) = ts(ji,jj,jk,jp_sal,Kmm) * 35.0 / 35.16504
END_3D
ELSE
salinprac(:,:,:) = ts(:,:,:,jp_sal,Kmm)
DO_3D( 0, 0, 0, 0, 1, jpk )
salinprac(ji,jj,jk) = ts(ji,jj,jk,jp_sal,Kmm)
END_3D
ENDIF
!
......@@ -179,7 +183,7 @@ CONTAINS
! potential temperature to in situ temperature. The errors is less than
! 0.04°C relative to an exact computation
! ---------------------------------------------------------------------
DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpk )
DO_3D( 0, 0, 0, 0, 1, jpk )
zpres = gdept(ji,jj,jk,Kmm) / 1000.
za1 = 0.04 * ( 1.0 + 0.185 * ts(ji,jj,jk,jp_tem,Kmm) + 0.035 * (salinprac(ji,jj,jk) - 35.0) )
za2 = 0.0075 * ( 1.0 - ts(ji,jj,jk,jp_tem,Kmm) / 30.0 )
......@@ -188,7 +192,7 @@ CONTAINS
!
! CHEMICAL CONSTANTS - SURFACE LAYER
! ----------------------------------
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
DO_2D( 0, 0, 0, 0 )
! ! SET ABSOLUTE TEMPERATURE
ztkel = tempis(ji,jj,1) + 273.15
zt = ztkel * 0.01
......@@ -216,7 +220,7 @@ CONTAINS
! OXYGEN SOLUBILITY - DEEP OCEAN
! -------------------------------
DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpk )
DO_3D( 0, 0, 0, 0, 1, jpk )
ztkel = tempis(ji,jj,jk) + 273.15
zsal = salinprac(ji,jj,jk) + ( 1.- tmask(ji,jj,jk) ) * 35.
zsal2 = zsal * zsal
......@@ -235,7 +239,7 @@ CONTAINS
! CHEMICAL CONSTANTS - DEEP OCEAN
! -------------------------------
DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpk )
DO_3D( 0, 0, 0, 0, 1, jpk )
! SET PRESSION ACCORDING TO SAUNDER (1980)
zplat = SIN ( ABS(gphit(ji,jj)*3.141592654/180.) )
zc1 = 5.92E-3 + zplat**2 * 5.25E-3
......@@ -452,7 +456,7 @@ CONTAINS
!! and the 2nd order approximation does not have
!! a solution
!!---------------------------------------------------------------------
REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(OUT) :: p_hini
REAL(wp), DIMENSION(A2D(0),jpk), INTENT(OUT) :: p_hini
INTEGER, INTENT(in) :: Kbb ! time level indices
INTEGER :: ji, jj, jk
REAL(wp) :: zca1, zba1
......@@ -463,7 +467,7 @@ CONTAINS
IF( ln_timing ) CALL timing_start('ahini_for_at')
!
DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpk )
DO_3D( 0, 0, 0, 0, 1, jpk )
zrhd = 1._wp / ( rhd(ji,jj,jk) + 1. )
p_alkcb = tr(ji,jj,jk,jptal,Kbb) * zrhd
p_dictot = tr(ji,jj,jk,jpdic,Kbb) * zrhd
......@@ -512,13 +516,13 @@ CONTAINS
! inf(TA - [OH-] + [H+]) and sup(TA - [OH-] + [H+])
! Argument variables
REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(OUT) :: p_alknw_inf
REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(OUT) :: p_alknw_sup
REAL(wp), DIMENSION(A2D(0),jpk), INTENT(OUT) :: p_alknw_inf
REAL(wp), DIMENSION(A2D(0),jpk), INTENT(OUT) :: p_alknw_sup
INTEGER, INTENT(in) :: Kbb ! time level indices
INTEGER :: ji, jj, jk
REAL(wp) :: zrhd
DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpk )
DO_3D( 0, 0, 0, 0, 1, jpk )
zrhd = 1._wp / ( rhd(ji,jj,jk) + 1. )
p_alknw_inf(ji,jj,jk) = -tr(ji,jj,jk,jppo4,Kbb) * zrhd - sulfat(ji,jj,jk) &
& - fluorid(ji,jj,jk)
......@@ -536,8 +540,8 @@ CONTAINS
! Argument variables
!--------------------
REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(IN) :: p_hini
REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(OUT) :: zhi
REAL(wp), DIMENSION(A2D(0),jpk), INTENT(IN) :: p_hini
REAL(wp), DIMENSION(A2D(0),jpk), INTENT(OUT) :: zhi
INTEGER, INTENT(in) :: Kbb ! time level indices
! Local variables
......@@ -557,17 +561,17 @@ CONTAINS
REAL(wp) :: zrhd, p_alktot, zdic, zbot, zpt, zst, zft, zsit
LOGICAL :: l_exitnow
REAL(wp), PARAMETER :: pz_exp_threshold = 1.0
REAL(wp), DIMENSION(jpi,jpj,jpk) :: zalknw_inf, zalknw_sup, rmask, zh_min, zh_max, zeqn_absmin
REAL(wp), DIMENSION(A2D(0),jpk) :: zalknw_inf, zalknw_sup, rmask, zh_min, zh_max, zeqn_absmin
IF( ln_timing ) CALL timing_start('solve_at_general')
CALL anw_infsup( zalknw_inf, zalknw_sup, Kbb )
rmask(:,:,:) = tmask(:,:,:)
rmask(A2D(0),1:jpk) = tmask(A2D(0),1:jpk)
zhi(:,:,:) = 0.
! TOTAL H+ scale: conversion factor for Htot = aphscale * Hfree
DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpk )
DO_3D( 0, 0, 0, 0, 1, jpk )
IF (rmask(ji,jj,jk) == 1.) THEN
zrhd = 1._wp / ( rhd(ji,jj,jk) + 1. )
p_alktot = tr(ji,jj,jk,jptal,Kbb) * zrhd
......@@ -597,7 +601,7 @@ CONTAINS
zeqn_absmin(:,:,:) = HUGE(1._wp)
DO jn = 1, jp_maxniter_atgen
DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpk )
DO_3D( 0, 0, 0, 0, 1, jpk )
IF (rmask(ji,jj,jk) == 1.) THEN
zrhd = 1._wp / ( rhd(ji,jj,jk) + 1. )
p_alktot = tr(ji,jj,jk,jptal,Kbb) * zrhd
......@@ -797,17 +801,17 @@ CONTAINS
ierr(:) = 0
ALLOCATE( sio3eq(jpi,jpj,jpk), fekeq(jpi,jpj,jpk), chemc(jpi,jpj,3), chemo2(jpi,jpj,jpk), STAT=ierr(1) )
ALLOCATE( sio3eq(A2D(0),jpk), fekeq(A2D(0),jpk), chemc(A2D(0),3), chemo2(A2D(0),jpk), STAT=ierr(1) )
ALLOCATE( akb3(jpi,jpj,jpk) , tempis(jpi, jpj, jpk), &
& akw3(jpi,jpj,jpk) , borat (jpi,jpj,jpk) , &
& aks3(jpi,jpj,jpk) , akf3(jpi,jpj,jpk) , &
& ak1p3(jpi,jpj,jpk) , ak2p3(jpi,jpj,jpk) , &
& ak3p3(jpi,jpj,jpk) , aksi3(jpi,jpj,jpk) , &
& fluorid(jpi,jpj,jpk) , sulfat(jpi,jpj,jpk) , &
& salinprac(jpi,jpj,jpk), STAT=ierr(2) )
ALLOCATE( akb3(A2D(0),jpk) , tempis(A2D(0),jpk), &
& akw3(A2D(0),jpk) , borat (A2D(0),jpk) , &
& aks3(A2D(0),jpk) , akf3(A2D(0),jpk) , &
& ak1p3(A2D(0),jpk) , ak2p3(A2D(0),jpk) , &
& ak3p3(A2D(0),jpk) , aksi3(A2D(0),jpk) , &
& fluorid(A2D(0),jpk) , sulfat(A2D(0),jpk) , &
& salinprac(A2D(0),jpk), STAT=ierr(2) )
ALLOCATE( fesol(jpi,jpj,jpk,5), STAT=ierr(3) )
ALLOCATE( fesol(A2D(0),jpk,5), STAT=ierr(3) )
!* Variable for chemistry of the CO2 cycle
p4z_che_alloc = MAXVAL( ierr )
......
......@@ -61,33 +61,42 @@ CONTAINS
REAL(wp) :: zprecip, zprecipno3, zconsfe, za1
REAL(wp) :: zrfact2
CHARACTER (len=25) :: charout
REAL(wp), DIMENSION(jpi,jpj,jpk) :: zTL1, zFe3, ztotlig, zfeprecip, zFeL1, zfecoll
REAL(wp), DIMENSION(jpi,jpj,jpk) :: zcoll3d, zscav3d, zlcoll3d
REAL(wp), DIMENSION(A2D(0),jpk) :: zTL1, zFe3, ztotlig, zfeprecip, zFeL1, zfecoll
REAL(wp), DIMENSION(A2D(0),jpk) :: zcoll3d, zscav3d, zlcoll3d
REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: zw3d
!!---------------------------------------------------------------------
!
IF( ln_timing ) CALL timing_start('p4z_fechem')
!
zFe3 (:,:,jpk) = 0.
zFeL1 (:,:,jpk) = 0.
zTL1 (:,:,jpk) = 0.
zfeprecip(:,:,jpk) = 0.
zcoll3d (:,:,jpk) = 0.
zscav3d (:,:,jpk) = 0.
zlcoll3d (:,:,jpk) = 0.
zfecoll (:,:,jpk) = 0.
xfecolagg(:,:,jpk) = 0.
xcoagfe (:,:,jpk) = 0.
! zFe3 (:,:,jpk) = 0.
! zFeL1 (:,:,jpk) = 0.
! zTL1 (:,:,jpk) = 0.
! zfeprecip(:,:,jpk) = 0.
! zcoll3d (:,:,jpk) = 0.
! zscav3d (:,:,jpk) = 0.
! zlcoll3d (:,:,jpk) = 0.
! zfecoll (:,:,jpk) = 0.
! xfecolagg(:,:,jpk) = 0.
! xcoagfe (:,:,jpk) = 0.
!
! Total ligand concentration : Ligands can be chosen to be constant or variable
! Parameterization from Pham and Ito (2018)
! -------------------------------------------------
xfecolagg(:,:,:) = ligand * 1E9 + MAX(0., chemo2(:,:,:) - tr(:,:,:,jpoxy,Kbb) ) / 400.E-6
DO_3D( 0, 0, 0, 0, 1, jpkm1)
xfecolagg(ji,jj,jk) = ligand * 1E9 + MAX(0., chemo2(ji,jj,jk) - tr(ji,jj,jk,jpoxy,Kbb) ) / 400.E-6
END_3D
IF( ln_ligvar ) THEN
ztotlig(:,:,:) = 0.09 * 0.667 * tr(:,:,:,jpdoc,Kbb) * 1E6 + xfecolagg(:,:,:)
ztotlig(:,:,:) = MIN( ztotlig(:,:,:), 10. )
DO_3D( 0, 0, 0, 0, 1, jpkm1)
ztotlig(ji,jj,jk) = 0.09 * 0.667 * tr(ji,jj,jk,jpdoc,Kbb) * 1E6 + xfecolagg(ji,jj,jk)
ztotlig(ji,jj,jk) = MIN( ztotlig(ji,jj,jk), 10. )
END_3D
ELSE
IF( ln_ligand ) THEN ; ztotlig(:,:,:) = tr(:,:,:,jplgw,Kbb) * 1E9
ELSE ; ztotlig(:,:,:) = ligand * 1E9
IF( ln_ligand ) THEN
DO_3D( 0, 0, 0, 0, 1, jpkm1)
ztotlig(ji,jj,jk) = tr(ji,jj,jk,jplgw,Kbb) * 1E9
END_3D
ELSE
ztotlig(:,:,:) = ligand * 1E9
ENDIF
ENDIF
......@@ -96,7 +105,7 @@ CONTAINS
! This model is based on one ligand, Fe2+ and Fe3+
! Chemistry is supposed to be fast enough to be at equilibrium
! ------------------------------------------------------------
DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1)
DO_3D( 0, 0, 0, 0, 1, jpkm1)
zTL1(ji,jj,jk) = ztotlig(ji,jj,jk)
zkeq = fekeq(ji,jj,jk)
zklight = 4.77E-7 * etot(ji,jj,jk) * 0.5 / ( 10**(-6.3) )
......@@ -109,7 +118,9 @@ CONTAINS
zFeL1(ji,jj,jk) = MAX( 0., tr(ji,jj,jk,jpfer,Kbb) - zFe3(ji,jj,jk) )
END_3D
!
plig(:,:,:) = MAX( 0., ( zFeL1(:,:,:) / ( tr(:,:,:,jpfer,Kbb) + rtrn ) ) )
DO_3D( 0, 0, 0, 0, 1, jpkm1)
plig(ji,jj,jk) = MAX( 0., ( zFeL1(ji,jj,jk) / ( tr(ji,jj,jk,jpfer,Kbb) + rtrn ) ) )
END_3D
!
zdust = 0. ! if no dust available
......@@ -125,16 +136,25 @@ CONTAINS
! to coagulate
! ----------------------------------------------------------------------
IF (ln_ligand) THEN
zfecoll(:,:,:) = 0.5 * zFeL1(:,:,:) * MAX(0., tr(:,:,:,jplgw,Kbb) - xfecolagg(:,:,:) * 1.0E-9 ) / ( tr(:,:,:,jplgw,Kbb) + rtrn )
DO_3D( 0, 0, 0, 0, 1, jpkm1)
zfecoll(ji,jj,jk) = 0.5 * zFeL1(ji,jj,jk) * MAX(0., tr(ji,jj,jk,jplgw,Kbb) - xfecolagg(ji,jj,jk) * 1.0E-9 ) &
& / ( tr(ji,jj,jk,jplgw,Kbb) + rtrn )
END_3D
ELSE
IF (ln_ligvar) THEN
zfecoll(:,:,:) = 0.5 * zFeL1(:,:,:) * MAX(0., tr(:,:,:,jplgw,Kbb) - xfecolagg(:,:,:) * 1.0E-9 ) / ( tr(:,:,:,jplgw,Kbb) + rtrn )
DO_3D( 0, 0, 0, 0, 1, jpkm1)
!Ch Buguuuugggggggg : tr(ji,jj,jk,jplgw ?????????
zfecoll(ji,jj,jk) = 0.5 * zFeL1(ji,jj,jk) * MAX(0., tr(ji,jj,jk,jplgw,Kbb) - xfecolagg(ji,jj,jk) * 1.0E-9 ) &
& / ( tr(ji,jj,jk,jplgw,Kbb) + rtrn )
END_3D
ELSE
zfecoll(:,:,:) = 0.5 * zFeL1(:,:,:)
DO_3D( 0, 0, 0, 0, 1, jpkm1)
zfecoll(ji,jj,jk) = 0.5 * zFeL1(ji,jj,jk)
END_3D
ENDIF
ENDIF
DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1)
DO_3D( 0, 0, 0, 0, 1, jpkm1)
! Scavenging rate of iron. This scavenging rate depends on the load of particles of sea water.
! This parameterization assumes a simple second order kinetics (k[Particles][Fe]).
! Scavenging onto dust is also included as evidenced from the DUNE experiments.
......@@ -199,20 +219,59 @@ CONTAINS
!
! Define the bioavailable fraction of iron
! ----------------------------------------
biron(:,:,:) = tr(:,:,:,jpfer,Kbb)
DO_3D( 0, 0, 0, 0, 1, jpkm1)
biron(ji,jj,jk) = tr(ji,jj,jk,jpfer,Kbb)
END_3D
!
! Output of some diagnostics variables
! ---------------------------------
IF( lk_iomput .AND. knt == nrdttrc ) THEN
zrfact2 = 1.e3 * rfact2r ! conversion from mol/L/timestep into mol/m3/s
IF( iom_use("Fe3") ) CALL iom_put("Fe3" , zFe3 (:,:,:) * tmask(:,:,:) ) ! Fe3+
IF( iom_use("FeL1") ) CALL iom_put("FeL1" , zFeL1 (:,:,:) * tmask(:,:,:) ) ! FeL1
IF( iom_use("TL1") ) CALL iom_put("TL1" , zTL1 (:,:,:) * tmask(:,:,:) ) ! TL1
IF( iom_use("Totlig") ) CALL iom_put("Totlig" , ztotlig(:,:,:) * tmask(:,:,:) ) ! TL
IF( iom_use("Biron") ) CALL iom_put("Biron" , biron (:,:,:) * 1e9 * tmask(:,:,:) ) ! biron
IF( iom_use("FESCAV") ) CALL iom_put("FESCAV" , zscav3d(:,:,:) * 1e9 * tmask(:,:,:) * zrfact2 )
IF( iom_use("FECOLL") ) CALL iom_put("FECOLL" , zcoll3d(:,:,:) * 1e9 * tmask(:,:,:) * zrfact2 )
IF( iom_use("FEPREC") ) CALL iom_put("FEPREC" , zfeprecip(:,:,:) *1e9*tmask(:,:,:)*zrfact2 )
IF( lk_iomput .AND. knt == nrdttrc ) THEN
!
ALLOCATE( zw3d(A2D(0),jpk) ) ; zw3d(A2D(0),jpk) = 0._wp
zrfact2 = 1.e3 * rfact2r ! conversion from mol/L/timestep into mol/m3/s
!
IF( iom_use ( "Fe3" ) ) THEN ! Fe3+
zw3d(A2D(0),1:jpkm1) = zFe3(A2D(0),1:jpkm1) * tmask(A2D(0),1:jpkm1)
CALL iom_put( "Fe3", zw3d )
ENDIF
!
IF( iom_use ( "FeL1" ) ) THEN ! FeL1
zw3d(A2D(0),1:jpkm1) = zFeL1(A2D(0),1:jpkm1) * tmask(A2D(0),1:jpkm1)
CALL iom_put( "FeL1", zw3d )
ENDIF
!
IF( iom_use ( "TL1" ) ) THEN ! TL1
zw3d(A2D(0),1:jpkm1) = zTL1(A2D(0),1:jpkm1) * tmask(A2D(0),1:jpkm1)
CALL iom_put( "TL1", zw3d )
ENDIF
!
IF( iom_use ( "Totlig" ) ) THEN ! Totlig
zw3d(A2D(0),1:jpkm1) = ztotlig(A2D(0),1:jpkm1) * tmask(A2D(0),1:jpkm1)
CALL iom_put( "Totlig", zw3d )
ENDIF
!
IF( iom_use ( "Biron" ) ) THEN ! biron
zw3d(A2D(0),1:jpkm1) = biron(A2D(0),1:jpkm1) * tmask(A2D(0),1:jpkm1)
CALL iom_put( "Biron", zw3d )
ENDIF
!
IF( iom_use ( "FESCAV" ) ) THEN ! FESCAV
zw3d(A2D(0),1:jpkm1) = zscav3d(A2D(0),1:jpkm1) * tmask(A2D(0),1:jpkm1) * zrfact2
CALL iom_put( "FESCAV", zw3d )
ENDIF
!
IF( iom_use ( "FECOLL" ) ) THEN ! FECOLL
zw3d(A2D(0),1:jpkm1) = zcoll3d(A2D(0),1:jpkm1) * tmask(A2D(0),1:jpkm1) * zrfact2
CALL iom_put( "FECOLL", zw3d )
ENDIF
!
IF( iom_use ( "FEPREC" ) ) THEN ! FEPREC
zw3d(A2D(0),1:jpkm1) = zfeprecip(A2D(0),1:jpkm1) * tmask(A2D(0),1:jpkm1) * zrfact2
CALL iom_put( "FEPREC", zw3d )
ENDIF
!
DEALLOCATE( zw3d )
!
ENDIF
IF(sn_cfctl%l_prttrc) THEN ! print mean trends (used for debugging)
......
......@@ -83,7 +83,8 @@ CONTAINS
REAL(wp) :: zph, zdic, zsch_o2, zsch_co2
REAL(wp) :: zyr_dec, zdco2dt
CHARACTER (len=25) :: charout
REAL(wp), DIMENSION(jpi,jpj) :: zkgco2, zkgo2, zh2co3, zoflx, zpco2atm, zpco2oce
REAL(wp), DIMENSION(A2D(0)) :: zkgco2, zkgo2, zh2co3, zoflx, zpco2atm, zpco2oce
REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zw2d
!!---------------------------------------------------------------------
!
IF( ln_timing ) CALL timing_start('p4z_flx')
......@@ -108,9 +109,13 @@ CONTAINS
satmco2(:,:) = atcco2
ENDIF
IF( l_co2cpl ) satmco2(:,:) = atm_co2(:,:)
IF( l_co2cpl ) THEN
DO_2D( 0, 0, 0, 0 )
satmco2(ji,jj) = atm_co2(ji,jj)
END_2D
ENDIF
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
DO_2D( 0, 0, 0, 0 )
! DUMMY VARIABLES FOR DIC, H+, AND BORATE
zrhd = rhd(ji,jj,1) + 1._wp
zdic = tr(ji,jj,1,jpdic,Kbb)
......@@ -126,7 +131,7 @@ CONTAINS
! FIRST COMPUTE GAS EXCHANGE COEFFICIENTS
! -------------------------------------------
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
DO_2D( 0, 0, 0, 0 )
ztc = MIN( 35., ts(ji,jj,1,jp_tem,Kmm) )
ztc2 = ztc * ztc
ztc3 = ztc * ztc2
......@@ -145,7 +150,7 @@ CONTAINS
END_2D
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
DO_2D( 0, 0, 0, 0 )
ztkel = tempis(ji,jj,1) + 273.15
zsal = salinprac(ji,jj,1) + ( 1.- tmask(ji,jj,1) ) * 35.
zvapsw = EXP(24.4543 - 67.4509*(100.0/ztkel) - 4.8489*LOG(ztkel/100) - 0.000544*zsal)
......@@ -171,11 +176,15 @@ CONTAINS
END_2D
IF( iom_use("tcflx") .OR. iom_use("tcflxcum") .OR. kt == nitrst &
& .OR. (ln_check_mass .AND. kt == nitend) ) &
t_oce_co2_flx = glob_sum( 'p4zflx', oce_co2(:,:) * e1e2t(:,:) * 1000. ) ! Total Flux of Carbon
t_oce_co2_flx_cum = t_oce_co2_flx_cum + t_oce_co2_flx ! Cumulative Total Flux of Carbon
! t_atm_co2_flx = glob_sum( 'p4zflx', satmco2(:,:) * e1e2t(:,:) ) ! Total atmospheric pCO2
t_atm_co2_flx = atcco2 ! Total atmospheric pCO2
& .OR. (ln_check_mass .AND. kt == nitend) ) THEN
ALLOCATE( zw2d(A2D(0)) )
zw2d(A2D(0)) = oce_co2(A2D(0)) * e1e2t(A2D(0)) * 1000._wp
t_oce_co2_flx = glob_sum( 'p4zflx', zw2d(:,:) ) ! Total Flux of Carbon
t_oce_co2_flx_cum = t_oce_co2_flx_cum + t_oce_co2_flx ! Cumulative Total Flux of Carbon
! t_atm_co2_flx = glob_sum( 'p4zflx', satmco2(:,:) * e1e2t(:,:) ) ! Total atmospheric pCO2
t_atm_co2_flx = atcco2 ! Total atmospheric pCO2
DEALLOCATE( zw2d )
ENDIF
IF(sn_cfctl%l_prttrc) THEN ! print mean trends (used for debugging)
WRITE(charout, FMT="('flx ')")
......@@ -184,15 +193,53 @@ CONTAINS
ENDIF
IF( lk_iomput .AND. knt == nrdttrc ) THEN
CALL iom_put( "AtmCo2" , satmco2(:,:) * tmask(:,:,1) ) ! Atmospheric CO2 concentration
CALL iom_put( "Cflx" , oce_co2(:,:) * 1000. )
CALL iom_put( "Oflx" , zoflx(:,:) * 1000. )
CALL iom_put( "Kg" , zkgco2(:,:) * tmask(:,:,1) )
CALL iom_put( "Dpco2" , ( zpco2atm(:,:) - zpco2oce(:,:) ) * tmask(:,:,1) )
CALL iom_put( "pCO2sea" , zpco2oce(:,:) * tmask(:,:,1) )
CALL iom_put( "Dpo2" , ( atcox * patm(:,:) - atcox * tr(:,:,1,jpoxy,Kbb) / ( chemo2(:,:,1) + rtrn ) ) * tmask(:,:,1) )
CALL iom_put( "tcflx" , t_oce_co2_flx ) ! molC/s
CALL iom_put( "tcflxcum", t_oce_co2_flx_cum ) ! molC
!
ALLOCATE( zw2d(A2D(0)) )
!
IF( iom_use ( "AtmCo2" ) ) THEN ! Atmospheric CO2 concentration
zw2d(A2D(0)) = satmco2(A2D(0)) * tmask(A2D(0),1)
CALL iom_put( "AtmCo2", zw2d )
ENDIF
!
IF( iom_use ( "Cflx" ) ) THEN
zw2d(A2D(0)) = oce_co2(A2D(0)) * 1000._wp
CALL iom_put( "Cflx", zw2d )
ENDIF
!
IF( iom_use ( "Dpo2" ) ) THEN
zw2d(A2D(0)) = ( atcox * patm(A2D(0)) - atcox * tr(A2D(0),1,jpoxy,Kbb) / ( chemo2(A2D(0),1) + rtrn ) ) * tmask(A2D(0),1)
CALL iom_put( "Dpo2", zw2d )
ENDIF
!
IF( iom_use ( "tcflx" ) ) THEN ! molC/s
CALL iom_put( "tcflx" , t_oce_co2_flx )
ENDIF
!
IF( iom_use ( "tcflxcum" ) ) THEN ! molC
CALL iom_put( "tcflxcum", t_oce_co2_flx_cum ) ! molC
ENDIF
!
IF( iom_use ( "Oflx" ) ) THEN
zw2d(A2D(0)) = zoflx(A2D(0)) * 1000.
CALL iom_put( "Oflx", zw2d )
ENDIF
!
IF( iom_use ( "Kg" ) ) THEN
zw2d(A2D(0)) = zkgco2(A2D(0)) * tmask(A2D(0),1)
CALL iom_put( "Kg", zw2d )
ENDIF
!
IF( iom_use ( "Dpco2" ) ) THEN
zw2d(A2D(0)) = ( zpco2atm(A2D(0)) - zpco2oce(A2D(0)) ) * tmask(A2D(0),1)
CALL iom_put( "Dpco2", zw2d )
ENDIF
!
IF( iom_use ( "pCO2sea" ) ) THEN
zw2d(A2D(0)) = zpco2oce(A2D(0)) * tmask(A2D(0),1)
CALL iom_put( "pCO2sea", zw2d )
ENDIF
!
DEALLOCATE( zw2d )
ENDIF
!
IF( ln_timing ) CALL timing_stop('p4z_flx')
......@@ -267,7 +314,7 @@ CONTAINS
IF(lwp) WRITE(numout,*) ' Spatialized Atmospheric pCO2 from an external file'
ENDIF
!
oce_co2(:,:) = 0._wp ! Initialization of Flux of Carbon
! oce_co2(:,:) = 0._wp ! Initialization of Flux of Carbon
t_oce_co2_flx = 0._wp
t_atm_co2_flx = 0._wp
!
......@@ -288,6 +335,7 @@ CONTAINS
CHARACTER(len=100) :: cn_dir ! Root directory for location of ssr files
TYPE(FLD_N) :: sn_patm ! informations about the fields to be read
TYPE(FLD_N) :: sn_atmco2 ! informations about the fields to be read
INTEGER :: ji, jj
!!
NAMELIST/nampisatm/ ln_presatm, ln_presatmco2, sn_patm, sn_atmco2, cn_dir
!!----------------------------------------------------------------------
......@@ -337,12 +385,16 @@ CONTAINS
!
IF( ln_presatm ) THEN
CALL fld_read( kt, 1, sf_patm ) !* input Patm provided at kt + 1/2
patm(:,:) = sf_patm(1)%fnow(:,:,1)/101325.0 ! atmospheric pressure
DO_2D( 0, 0, 0, 0 )
patm(ji,jj) = sf_patm(1)%fnow(ji,jj,1)/101325.0 ! atmospheric pressure
END_2D
ENDIF
!
IF( ln_presatmco2 ) THEN
CALL fld_read( kt, 1, sf_atmco2 ) !* input atmco2 provided at kt + 1/2
satmco2(:,:) = sf_atmco2(1)%fnow(:,:,1) ! atmospheric pressure
DO_2D( 0, 0, 0, 0 )
satmco2(ji,jj) = sf_atmco2(1)%fnow(ji,jj,1) ! atmospheric pressure
END_2D
ELSE
satmco2(:,:) = atcco2 ! Initialize atmco2 if no reading from a file
ENDIF
......@@ -354,7 +406,7 @@ CONTAINS
!!----------------------------------------------------------------------
!! *** ROUTINE p4z_flx_alloc ***
!!----------------------------------------------------------------------
ALLOCATE( satmco2(jpi,jpj), patm(jpi,jpj), STAT=p4z_flx_alloc )
ALLOCATE( satmco2(A2D(0)), patm(A2D(0)), STAT=p4z_flx_alloc )
!
IF( p4z_flx_alloc /= 0 ) CALL ctl_stop( 'STOP', 'p4z_flx_alloc : failed to allocate arrays' )
!
......
......@@ -44,16 +44,18 @@ CONTAINS
!
! Computation of phyto and zoo metabolic rate
! -------------------------------------------
! Generic temperature dependence (Eppley, 1972)
tgfunc (:,:,:) = EXP( 0.0631 * ts(:,:,:,jp_tem,Kmm) )
! Temperature dependence of mesozooplankton (Buitenhuis et al. (2005))
tgfunc2(:,:,:) = EXP( 0.0761 * ts(:,:,:,jp_tem,Kmm) )
DO_3D( 0, 0, 0, 0, 1, jpk )
! Generic temperature dependence (Eppley, 1972)
tgfunc (ji,jj,jk) = EXP( 0.0631 * ts(ji,jj,jk,jp_tem,Kmm) )
! Temperature dependence of mesozooplankton (Buitenhuis et al. (2005))
tgfunc2(ji,jj,jk) = EXP( 0.0761 * ts(ji,jj,jk,jp_tem,Kmm) )
END_3D
! Computation of the silicon dependant half saturation constant for silica uptake
! This is based on an old study by Pondaven et al. (1998)
! --------------------------------------------------------------------------------
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
DO_2D( 0, 0, 0, 0 )
zvar = tr(ji,jj,1,jpsil,Kbb) * tr(ji,jj,1,jpsil,Kbb)
xksimax(ji,jj) = MAX( xksimax(ji,jj), ( 1.+ 7.* zvar / ( xksilim * xksilim + zvar ) ) * 1e-6 )
END_2D
......@@ -73,14 +75,14 @@ CONTAINS
zcodel = ASIN( SIN( zrum * rpi * 2._wp ) * SIN( rad * 23.5_wp ) )
! day length in hours
strn(:,:) = 0.
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
! strn(:,:) = 0.
DO_2D( 0, 0, 0, 0 )
zargu = TAN( zcodel ) * TAN( gphit(ji,jj) * rad )
zargu = MAX( -1., MIN( 1., zargu ) )
strn(ji,jj) = MAX( 0.0, 24. - 2. * ACOS( zargu ) / rad / 15. )
END_2D
!
DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1 )
DO_3D( 0, 0, 0, 0, 1, jpkm1 )
! denitrification factor computed from O2 levels
! This factor diagnoses below which level of O2 denitrification
! is active
......
......@@ -47,13 +47,13 @@ CONTAINS
INTEGER :: ji, jj, jk
REAL(wp) :: zlgwp, zlgwpr, zlgwr, zlablgw
REAL(wp) :: zlam1a, zlam1b, zaggliga, zligco
REAL(wp), DIMENSION(jpi,jpj,jpk) :: zligrem, zligpr, zligprod, zlcoll3d
REAL(wp), DIMENSION(A2D(0),jpk) :: zligrem, zligpr, zligprod, zlcoll3d
CHARACTER (len=25) :: charout
!!---------------------------------------------------------------------
!
IF( ln_timing ) CALL timing_start('p4z_ligand')
!
DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1)
DO_3D( 0, 0, 0, 0, 1, jpkm1)
!
! ------------------------------------------------------------------
! Remineralization of iron ligands
......@@ -89,16 +89,16 @@ CONTAINS
! ---------------------------------
IF( lk_iomput .AND. knt == nrdttrc ) THEN
IF( iom_use( "LIGREM" ) ) THEN
zligrem(:,:,jpk) = 0. ; CALL iom_put( "LIGREM", zligrem(:,:,:) * 1e9 * 1.e+3 * rfact2r * tmask(:,:,:) )
zligrem(A2D(0),jpk) = 0. ; CALL iom_put( "LIGREM", zligrem(A2D(0),:) * 1e9 * 1.e+3 * rfact2r * tmask(A2D(0),:) )
ENDIF
IF( iom_use( "LIGPR" ) ) THEN
zligpr(:,:,jpk) = 0. ; CALL iom_put( "LIGPR" , zligpr(:,:,:) * 1e9 * 1.e+3 * rfact2r * tmask(:,:,:) )
zligpr(A2D(0),jpk) = 0. ; CALL iom_put( "LIGPR" , zligpr(A2D(0),:) * 1e9 * 1.e+3 * rfact2r * tmask(A2D(0),:) )
ENDIF
IF( iom_use( "LPRODR" ) ) THEN
zligprod(:,:,jpk) = 0. ; CALL iom_put( "LPRODR", zligprod(:,:,:) * 1e9 * 1.e+3 * rfact2r * tmask(:,:,:) )
zligprod(A2D(0),jpk) = 0. ; CALL iom_put( "LPRODR", zligprod(A2D(0),:) * 1e9 * 1.e+3 * rfact2r * tmask(A2D(0),:) )
ENDIF
IF( iom_use( "LGWCOLL" ) ) THEN
zlcoll3d(:,:,jpk) = 0. ; CALL iom_put( "LGWCOLL", zlcoll3d(:,:,:) * 1.e9 * 1.e+3 * rfact2r * tmask(:,:,:) )
zlcoll3d(A2D(0),jpk) = 0. ; CALL iom_put( "LGWCOLL", zlcoll3d(A2D(0),:) * 1.e9 * 1.e+3 * rfact2r * tmask(A2D(0),:) )
ENDIF
ENDIF
!
......
......@@ -98,13 +98,14 @@ CONTAINS
REAL(wp) :: zconc1d, zconc1dnh4, zconc0n, zconc0nnh4
REAL(wp) :: fananof, fadiatf, znutlim, zfalim
REAL(wp) :: znutlimtot, zlimno3, zlimnh4, zbiron
REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: zw3d
!!---------------------------------------------------------------------
!
IF( ln_timing ) CALL timing_start('p4z_lim')
!
sizena(:,:,:) = 1.0 ; sizeda(:,:,:) = 1.0
!
DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1)
DO_3D( 0, 0, 0, 0, 1, jpkm1)
! Computation of a variable Ks for iron on diatoms taking into account
! that increasing biomass is made of generally bigger cells
......@@ -219,7 +220,7 @@ CONTAINS
! Size estimation of phytoplankton based on total biomass
! Assumes that larger biomass implies addition of larger cells
! ------------------------------------------------------------
DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1)
DO_3D( 0, 0, 0, 0, 1, jpkm1)
zcoef = tr(ji,jj,jk,jpphy,Kbb) - MIN(xsizephy, tr(ji,jj,jk,jpphy,Kbb) )
sizena(ji,jj,jk) = 1. + ( xsizern -1.0 ) * zcoef / ( xsizephy + zcoef )
zcoef = tr(ji,jj,jk,jpdia,Kbb) - MIN(xsizedia, tr(ji,jj,jk,jpdia,Kbb) )
......@@ -231,7 +232,7 @@ CONTAINS
! This is a purely adhoc formulation described in Aumont et al. (2015)
! This fraction depends on nutrient limitation, light, temperature
! --------------------------------------------------------------------
DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1)
DO_3D( 0, 0, 0, 0, 1, jpkm1)
zlim1 = xnanonh4(ji,jj,jk) + xnanono3(ji,jj,jk)
zlim2 = tr(ji,jj,jk,jppo4,Kbb) / ( tr(ji,jj,jk,jppo4,Kbb) + concnnh4 )
zlim3 = tr(ji,jj,jk,jpfer,Kbb) / ( tr(ji,jj,jk,jpfer,Kbb) + 6.E-11 )
......@@ -250,7 +251,7 @@ CONTAINS
xfracal(ji,jj,jk) = MAX( 0.02, xfracal(ji,jj,jk) )
END_3D
!
DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1)
DO_3D( 0, 0, 0, 0, 1, jpkm1)
! denitrification factor computed from O2 levels
nitrfac(ji,jj,jk) = MAX( 0.e0, 0.4 * ( 6.e-6 - tr(ji,jj,jk,jpoxy,Kbb) ) &
& / ( oxymin + tr(ji,jj,jk,jpoxy,Kbb) ) )
......@@ -263,13 +264,45 @@ CONTAINS
END_3D
!
IF( lk_iomput .AND. knt == nrdttrc ) THEN ! save output diagnostics
CALL iom_put( "xfracal", xfracal(:,:,:) * tmask(:,:,:) ) ! euphotic layer deptht
CALL iom_put( "LNnut" , xlimphy(:,:,:) * tmask(:,:,:) ) ! Nutrient limitation term
CALL iom_put( "LDnut" , xlimdia(:,:,:) * tmask(:,:,:) ) ! Nutrient limitation term
CALL iom_put( "LNFe" , xlimnfe(:,:,:) * tmask(:,:,:) ) ! Iron limitation term
CALL iom_put( "LDFe" , xlimdfe(:,:,:) * tmask(:,:,:) ) ! Iron limitation term
CALL iom_put( "SIZEN" , sizen (:,:,:) * tmask(:,:,:) ) ! Iron limitation term
CALL iom_put( "SIZED" , sized (:,:,:) * tmask(:,:,:) ) ! Iron limitation term
!
ALLOCATE( zw3d(A2D(0),jpk) ) ; zw3d(A2D(0),jpk) = 0._wp
!
IF( iom_use ( "xfracal" ) ) THEN ! euphotic layer deptht
zw3d(A2D(0),1:jpkm1) = xfracal(A2D(0),1:jpkm1) * tmask(A2D(0),1:jpkm1)
CALL iom_put( "xfracal", zw3d)
ENDIF
!
IF( iom_use ( "LNnut" ) ) THEN ! Nutrient limitation term
zw3d(A2D(0),1:jpkm1) = xlimphy(A2D(0),1:jpkm1) * tmask(A2D(0),1:jpkm1)
CALL iom_put( "LNnut", zw3d)
ENDIF
!
IF( iom_use ( "LDnut" ) ) THEN ! Nutrient limitation term
zw3d(A2D(0),1:jpkm1) = xlimdia(A2D(0),1:jpkm1) * tmask(A2D(0),1:jpkm1)
CALL iom_put( "LDnut", zw3d)
ENDIF
!
IF( iom_use ( "LNFe" ) ) THEN ! Iron limitation term
zw3d(A2D(0),1:jpkm1) = xlimnfe(A2D(0),1:jpkm1) * tmask(A2D(0),1:jpkm1)
CALL iom_put( "LNFe", zw3d)
ENDIF
!
IF( iom_use ( "LDFe" ) ) THEN ! Iron limitation term
zw3d(A2D(0),1:jpkm1) = xlimdfe(A2D(0),1:jpkm1) * tmask(A2D(0),1:jpkm1)
CALL iom_put( "LDFe", zw3d)
ENDIF
!
IF( iom_use ( "SIZEN" ) ) THEN ! Iron limitation term
zw3d(A2D(0),1:jpkm1) = sizen(A2D(0),1:jpkm1) * tmask(A2D(0),1:jpkm1)
CALL iom_put( "SIZEN", zw3d)
ENDIF
!
IF( iom_use ( "SIZED" ) ) THEN ! Iron limitation term
zw3d(A2D(0),1:jpkm1) = sized(A2D(0),1:jpkm1) * tmask(A2D(0),1:jpkm1)
CALL iom_put( "SIZED", zw3d)
ENDIF
!
DEALLOCATE( zw3d )
ENDIF
!
IF( ln_timing ) CALL timing_stop('p4z_lim')
......@@ -355,16 +388,16 @@ CONTAINS
!!----------------------------------------------------------------------
!* Biological arrays for phytoplankton growth
ALLOCATE( xnanono3(jpi,jpj,jpk), xdiatno3(jpi,jpj,jpk), &
& xnanonh4(jpi,jpj,jpk), xdiatnh4(jpi,jpj,jpk), &
& xnanopo4(jpi,jpj,jpk), xdiatpo4(jpi,jpj,jpk), &
& xnanofer(jpi,jpj,jpk), xdiatfer(jpi,jpj,jpk), &
& xlimphy (jpi,jpj,jpk), xlimdia (jpi,jpj,jpk), &
& xlimnfe (jpi,jpj,jpk), xlimdfe (jpi,jpj,jpk), &
& xlimbac (jpi,jpj,jpk), xlimbacl(jpi,jpj,jpk), &
& concnfe (jpi,jpj,jpk), concdfe (jpi,jpj,jpk), &
& xqfuncfecn(jpi,jpj,jpk), xqfuncfecd(jpi,jpj,jpk), &
& xlimsi (jpi,jpj,jpk), STAT=p4z_lim_alloc )
ALLOCATE( xnanono3(A2D(0),jpk), xdiatno3(A2D(0),jpk), &
& xnanonh4(A2D(0),jpk), xdiatnh4(A2D(0),jpk), &
& xnanopo4(A2D(0),jpk), xdiatpo4(A2D(0),jpk), &
& xnanofer(A2D(0),jpk), xdiatfer(A2D(0),jpk), &
& xlimphy (A2D(0),jpk), xlimdia (A2D(0),jpk), &
& xlimnfe (A2D(0),jpk), xlimdfe (A2D(0),jpk), &
& xlimbac (A2D(0),jpk), xlimbacl(A2D(0),jpk), &
& concnfe (A2D(0),jpk), concdfe (A2D(0),jpk), &
& xqfuncfecn(A2D(0),jpk), xqfuncfecd(A2D(0),jpk), &
& xlimsi (A2D(0),jpk), STAT=p4z_lim_alloc )
!
IF( p4z_lim_alloc /= 0 ) CALL ctl_stop( 'STOP', 'p4z_lim_alloc : failed to allocate arrays.' )
!
......
......@@ -32,7 +32,7 @@ MODULE p4zlys
REAL(wp), PUBLIC :: kdca !: diss. rate constant calcite
REAL(wp), PUBLIC :: nca !: order of reaction for calcite dissolution
INTEGER :: rmtss ! number of seconds per month
INTEGER :: rmtss ! number of seconds per month
!! * Module variables
REAL(wp) :: calcon = 1.03E-2 !: mean calcite concentration [Ca2+] in sea water [mole/kg solution]
......@@ -66,20 +66,26 @@ CONTAINS
REAL(wp) :: zdispot, zrhd, zcalcon
REAL(wp) :: zomegaca, zexcess, zexcess0, zkd
CHARACTER (len=25) :: charout
REAL(wp), DIMENSION(jpi,jpj,jpk) :: zco3, zcaldiss, zhinit, zhi, zco3sat
!! CE later REAL(wp), DIMENSION(A2D(0),jpk) :: zco3, zcaldiss, zhinit, zhi, zco3sat
!! because of the routine solve_at_general in p4zche.F90
REAL(wp), DIMENSION(A2D(0),jpk) :: zco3, zcaldiss, zco3sat
REAL(wp), DIMENSION(A2D(0),jpk) :: zhinit, zhi
REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: zw3d
!!---------------------------------------------------------------------
!
IF( ln_timing ) CALL timing_start('p4z_lys')
!
zhinit (:,:,:) = hi(:,:,:) / ( rhd(:,:,:) + 1._wp )
DO_3D( 0, 0, 0, 0, 1, jpkm1)
zhinit(ji,jj,jk) = hi(ji,jj,jk) / ( rhd(ji,jj,jk) + 1._wp )
END_3D
!
! -------------------------------------------
! COMPUTE [CO3--] and [H+] CONCENTRATIONS
! -------------------------------------------
CALL solve_at_general( zhinit, zhi, Kbb )
DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1)
DO_3D( 0, 0, 0, 0, 1, jpkm1)
zco3(ji,jj,jk) = tr(ji,jj,jk,jpdic,Kbb) * ak13(ji,jj,jk) * ak23(ji,jj,jk) / (zhi(ji,jj,jk)**2 &
& + ak13(ji,jj,jk) * zhi(ji,jj,jk) + ak13(ji,jj,jk) * ak23(ji,jj,jk) + rtrn )
hi (ji,jj,jk) = zhi(ji,jj,jk) * ( rhd(ji,jj,jk) + 1._wp )
......@@ -91,7 +97,7 @@ CONTAINS
! MGCO3)
! ---------------------------------------------------------
DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1)
DO_3D( 0, 0, 0, 0, 1, jpkm1)
! DEVIATION OF [CO3--] FROM SATURATION VALUE
! Salinity dependance in zomegaca and divide by rhd to have good units
......@@ -125,16 +131,32 @@ CONTAINS
!
IF( lk_iomput .AND. knt == nrdttrc ) THEN
CALL iom_put( "PH" , -1. * LOG10( MAX( hi(:,:,:), rtrn ) ) * tmask(:,:,:) )
IF( iom_use( "CO3" ) ) THEN
zco3(:,:,jpk) = 0. ; CALL iom_put( "CO3" , zco3(:,:,:) * 1.e+3 * tmask(:,:,:) )
ALLOCATE( zw3d(A2D(0),jpk) ) ; zw3d(A2D(0),jpk) = 0._wp
IF( iom_use( "PH" ) ) THEN
DO_3D( 0, 0, 0, 0, 1, jpkm1)
zw3d(ji,jj,jk) = -1. * LOG10( MAX( hi(ji,jj,jk), rtrn ) ) * tmask(ji,jj,jk)
END_3D
CALL iom_put( "PH" , zw3d )
ENDIF
IF( iom_use( "CO3" ) ) THEN ! bicarbonate
DO_3D( 0, 0, 0, 0, 1, jpkm1)
zw3d(ji,jj,jk) = zco3(ji,jj,jk) * 1.e+3 * tmask(ji,jj,jk)
END_3D
CALL iom_put( "CO3", zw3d )
ENDIF
IF( iom_use( "CO3sat" ) ) THEN
zco3sat(:,:,jpk) = 0. ; CALL iom_put( "CO3sat", zco3sat(:,:,:) * 1.e+3 * tmask(:,:,:) )
IF( iom_use( "CO3sat" ) ) THEN ! calcite saturation
DO_3D( 0, 0, 0, 0, 1, jpkm1)
zw3d(ji,jj,jk) = zco3sat(ji,jj,jk) * 1.e+3 * tmask(ji,jj,jk)
END_3D
CALL iom_put( "CO3sat", zw3d )
ENDIF
IF( iom_use( "DCAL" ) ) THEN
zcaldiss(:,:,jpk) = 0. ; CALL iom_put( "DCAL" , zcaldiss(:,:,:) * 1.e+3 * rfact2r * tmask(:,:,:) )
IF( iom_use( "DCAL" ) ) THEN ! calcite dissolution
DO_3D( 0, 0, 0, 0, 1, jpkm1)
zw3d(ji,jj,jk) = zcaldiss(ji,jj,jk) * 1.e+3 * rfact2r * tmask(ji,jj,jk)
END_3D
CALL iom_put( "DCAL", zw3d )
ENDIF
DEALLOCATE( zw3d )
ENDIF
!
IF(sn_cfctl%l_prttrc) THEN ! print mean trends (used for debugging)
......@@ -185,6 +207,9 @@ CONTAINS
! Number of seconds per month
rmtss = nyear_len(1) * rday / raamo
!
! CE not really needed ; tempory, shoub be removed when quotan( A2D(0),jpk )
excess(:,:,:) = 0._wp
!
END SUBROUTINE p4z_lys_init
!!======================================================================
......
......@@ -89,10 +89,11 @@ CONTAINS
REAL(wp) :: zgrazfffp, zgrazfffg, zgrazffep, zgrazffeg, zrum, zcodel, zargu, zval, zdep
REAL(wp) :: zsigma, zdiffdn, ztmp1, ztmp2, ztmp3, ztmp4, ztmptot, zmigthick
CHARACTER (len=25) :: charout
REAL(wp), DIMENSION(jpi,jpj,jpk) :: zgrazing, zfezoo2
REAL(wp), DIMENSION(jpi,jpj,jpk) :: zgrarem, zgraref, zgrapoc, zgrapof, zgrabsi
REAL(wp), DIMENSION(A2D(0),jpk) :: zgrazing, zfezoo2
REAL(wp), DIMENSION(A2D(0),jpk) :: zgrarem, zgraref, zgrapoc, zgrapof, zgrabsi
REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zgramigrem, zgramigref, zgramigpoc, zgramigpof
REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zgramigbsi
REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: zw3d
!!---------------------------------------------------------------------
!
IF( ln_timing ) CALL timing_start('p4z_meso')
......@@ -108,7 +109,7 @@ CONTAINS
! ---------------------------------------------
IF (ln_dvm_meso) CALL p4z_meso_depmig( Kbb, Kmm )
!
DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1)
DO_3D( 0, 0, 0, 0, 1, jpkm1)
zcompam = MAX( ( tr(ji,jj,jk,jpmes,Kbb) - 1.e-9 ), 0.e0 )
zfact = xstep * tgfunc2(ji,jj,jk) * zcompam
......@@ -345,7 +346,7 @@ CONTAINS
! This fraction is sumed over the euphotic zone and is removed from
! the fluxes driven by mesozooplankton in the euphotic zone.
! --------------------------------------------------------------------
DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpk)
DO_3D( 0, 0, 0, 0, 1, jpk)
zmigreltime = (1. - strn(ji,jj))
zmigthick = (1. - zmigreltime ) * e3t(ji,jj,jk,Kmm) * tmask(ji,jj,jk)
IF ( gdept(ji,jj,jk,Kmm) <= heup(ji,jj) ) THEN
......@@ -366,7 +367,7 @@ CONTAINS
! The inorganic and organic fluxes induced by migrating organisms are added at the
! the migration depth (corresponding indice is set by kmig)
! --------------------------------------------------------------------------------
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
DO_2D( 0, 0, 0, 0 )
IF( tmask(ji,jj,1) == 1.) THEN
jkt = kmig(ji,jj)
zdep = 1. / e3t(ji,jj,jkt,Kmm)
......@@ -387,7 +388,7 @@ CONTAINS
! Update the arrays TRA which contain the biological sources and sinks
! This only concerns the variables which are affected by DVM (inorganic
! nutrients, DOC agands, and particulate organic carbon).
DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpk)
DO_3D( 0, 0, 0, 0, 1, jpk)
tr(ji,jj,jk,jppo4,Krhs) = tr(ji,jj,jk,jppo4,Krhs) + zgrarem(ji,jj,jk) * sigma2
tr(ji,jj,jk,jpnh4,Krhs) = tr(ji,jj,jk,jpnh4,Krhs) + zgrarem(ji,jj,jk) * sigma2
tr(ji,jj,jk,jpdoc,Krhs) = tr(ji,jj,jk,jpdoc,Krhs) + zgrarem(ji,jj,jk) * ( 1. - sigma2 )
......@@ -408,11 +409,30 @@ CONTAINS
!
! Write the output
IF( lk_iomput .AND. knt == nrdttrc ) THEN
CALL iom_put( "PCAL" , prodcal(:,:,:) * 1.e+3 * rfact2r * tmask(:,:,:) ) ! Calcite production
CALL iom_put( "GRAZ2" , zgrazing(:,:,:) * 1.e+3 * rfact2r * tmask(:,:,:) ) ! Total grazing of phyto by zoo
CALL iom_put( "FEZOO2", zfezoo2(:,:,:) * 1e9 * 1.e+3 * rfact2r * tmask(:,:,:) )
IF( ln_ligand ) &
& CALL iom_put( "LPRODZ2", zgrarem(ji,jj,jk) * ( 1. - sigma2 ) * ldocz * 1e9 * 1.e+3 * rfact2r * tmask(:,:,:) )
!
ALLOCATE( zw3d(A2D(0),jpk) ) ; zw3d(A2D(0),jpk) = 0._wp
!
IF( iom_use ( "PCAL" ) ) THEN ! Calcite production
zw3d(A2D(0),1:jpkm1) = prodcal(A2D(0),1:jpkm1) * 1.e+3 * rfact2r * tmask(A2D(0),1:jpkm1)
CALL iom_put( "PCAL", zw3d )
ENDIF
!
IF( iom_use ( "GRAZ2" ) ) THEN ! Total grazing of phyto by zoo
zw3d(A2D(0),1:jpkm1) = zgrazing(A2D(0),1:jpkm1) * 1.e+3 * rfact2r * tmask(A2D(0),1:jpkm1)
CALL iom_put( "GRAZ2", zw3d )
ENDIF
!
IF( iom_use ( "FEZOO2" ) ) THEN
zw3d(A2D(0),1:jpkm1) = zfezoo2(A2D(0),1:jpkm1) * 1e9 * 1.e+3 * rfact2r * tmask(A2D(0),1:jpkm1)
CALL iom_put( "FEZOO2", zw3d )
ENDIF
!
IF( ln_ligand .AND. iom_use( "LPRODZ2" ) ) THEN
zw3d(A2D(0),1:jpkm1) = zgrarem(A2D(0),1:jpkm1) * ( 1. - sigma2 ) * 1.e+3 * rfact2r * tmask(A2D(0),1:jpkm1)
CALL iom_put( "LPRODZ2" , zw3d * ldocz * 1e9 )
ENDIF
!
DEALLOCATE( zw3d )
ENDIF
!
IF(sn_cfctl%l_prttrc) THEN ! print mean trends (used for debugging)
......@@ -502,7 +522,7 @@ CONTAINS
INTEGER :: ji, jj, jk
!
REAL(wp) :: ztotchl, z1dep
REAL(wp), DIMENSION(jpi,jpj) :: oxymoy, tempmoy, zdepmoy
REAL(wp), DIMENSION(A2D(0)) :: oxymoy, tempmoy, zdepmoy
!!---------------------------------------------------------------------
!
......@@ -517,7 +537,7 @@ CONTAINS
! Compute the averaged values of oxygen, temperature over the domain
! 150m to 500 m depth.
! ------------------------------------------------------------------
DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpk)
DO_3D( 0, 0, 0, 0, 1, jpk)
IF( tmask(ji,jj,jk) == 1.) THEN
IF( gdept(ji,jj,jk,Kmm) >= 150. .AND. gdept(ji,jj,jk,kmm) <= 500.) THEN
oxymoy(ji,jj) = oxymoy(ji,jj) + tr(ji,jj,jk,jpoxy,Kbb) * 1E6 * e3t(ji,jj,jk,Kmm)
......@@ -530,7 +550,7 @@ CONTAINS
! Compute the difference between surface values and the mean values in the mesopelagic
! domain
! ------------------------------------------------------------------------------------
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
DO_2D( 0, 0, 0, 0 )
z1dep = 1. / ( zdepmoy(ji,jj) + rtrn )
oxymoy(ji,jj) = tr(ji,jj,1,jpoxy,Kbb) * 1E6 - oxymoy(ji,jj) * z1dep
tempmoy(ji,jj) = ts(ji,jj,1,jp_tem,Kmm) - tempmoy(ji,jj) * z1dep
......@@ -539,7 +559,7 @@ CONTAINS
! Computation of the migration depth based on the parameterization of
! Bianchi et al. (2013)
! -------------------------------------------------------------------
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
DO_2D( 0, 0, 0, 0 )
IF( tmask(ji,jj,1) == 1. ) THEN
ztotchl = ( tr(ji,jj,1,jpnch,Kbb) + tr(ji,jj,1,jpdch,Kbb) ) * 1E6
depmig(ji,jj) = 398. - 0.56 * oxymoy(ji,jj) -115. * log10(ztotchl) + 0.36 * hmld(ji,jj) -2.4 * tempmoy(ji,jj)
......@@ -548,7 +568,7 @@ CONTAINS
!
! Computation of the corresponding jk indice
! ------------------------------------------
DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1)
DO_3D( 0, 0, 0, 0, 1, jpkm1)
IF( depmig(ji,jj) >= gdepw(ji,jj,jk,Kmm) .AND. depmig(ji,jj) < gdepw(ji,jj,jk+1,Kmm) ) THEN
kmig(ji,jj) = jk
ENDIF
......@@ -560,7 +580,7 @@ CONTAINS
! to 0. Thus, to avoid that problem, the migration depth is adjusted so
! that it falls above the OMZ
! -----------------------------------------------------------------------
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
DO_2D( 0, 0, 0, 0 )
IF( tr(ji,jj,kmig(ji,jj),jpoxy,Kbb) < 5E-6 ) THEN
DO jk = kmig(ji,jj),1,-1
IF( tr(ji,jj,jk,jpoxy,Kbb) >= 5E-6 .AND. tr(ji,jj,jk+1,jpoxy,Kbb) < 5E-6) THEN
......@@ -580,7 +600,7 @@ CONTAINS
!! *** ROUTINE p4z_meso_alloc ***
!!----------------------------------------------------------------------
!
ALLOCATE( depmig(jpi,jpj), kmig(jpi,jpj), STAT= p4z_meso_alloc )
ALLOCATE( depmig(A2D(0)), kmig(A2D(0)), STAT= p4z_meso_alloc )
!
IF( p4z_meso_alloc /= 0 ) CALL ctl_stop( 'STOP', 'p4z_meso_alloc : failed to allocate arrays.' )
!
......
......@@ -78,20 +78,15 @@ CONTAINS
REAL(wp) :: zrespz, ztortz, zgrasratf, zgrasratn
REAL(wp) :: zgraznc, zgrazpoc, zgrazdc, zgrazpof, zgrazdf, zgraznf
REAL(wp) :: zsigma, zdiffdn, ztmp1, ztmp2, ztmp3, ztmptot, zproport
REAL(wp), DIMENSION(jpi,jpj,jpk) :: zgrazing, zfezoo
REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: zzligprod
REAL(wp), DIMENSION(A2D(0),jpk) :: zgrazing, zfezoo
REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: zw3d
CHARACTER (len=25) :: charout
!!---------------------------------------------------------------------
!
IF( ln_timing ) CALL timing_start('p4z_micro')
!
IF (ln_ligand) THEN
ALLOCATE( zzligprod(jpi,jpj,jpk) )
zzligprod(:,:,:) = 0._wp
ENDIF
!
DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1)
DO_3D( 0, 0, 0, 0, 1, jpkm1)
zcompaz = MAX( ( tr(ji,jj,jk,jpzoo,Kbb) - 1.e-9 ), 0.e0 )
zfact = xstep * tgfunc2(ji,jj,jk) * zcompaz
......@@ -215,7 +210,6 @@ CONTAINS
!
IF( ln_ligand ) THEN
tr(ji,jj,jk,jplgw,Krhs) = tr(ji,jj,jk,jplgw,Krhs) + (zgrarem - zgrarsig) * ldocz
zzligprod(ji,jj,jk) = (zgrarem - zgrarsig) * ldocz
ENDIF
!
tr(ji,jj,jk,jpoxy,Krhs) = tr(ji,jj,jk,jpoxy,Krhs) - o2ut * zgrarsig
......@@ -257,15 +251,25 @@ CONTAINS
END_3D
!
IF( lk_iomput .AND. knt == nrdttrc ) THEN
IF( iom_use("GRAZ1") ) THEN ! Total grazing of phyto by zooplankton
zgrazing(:,:,jpk) = 0._wp ; CALL iom_put( "GRAZ1" , zgrazing(:,:,:) * 1.e+3 * rfact2r * tmask(:,:,:) )
ENDIF
IF( iom_use("FEZOO") ) THEN
zfezoo (:,:,jpk) = 0._wp ; CALL iom_put( "FEZOO", zfezoo(:,:,:) * 1e9 * 1.e+3 * rfact2r * tmask(:,:,:) )
ENDIF
IF( ln_ligand ) THEN
zzligprod(:,:,jpk) = 0._wp ; CALL iom_put( "LPRODZ", zzligprod(:,:,:) * 1e9 * 1.e+3 * rfact2r * tmask(:,:,:))
ENDIF
!
ALLOCATE( zw3d(A2D(0),jpk) ) ; zw3d(A2D(0),jpk) = 0._wp
!
IF( iom_use ( "GRAZ1" ) ) THEN ! Total grazing of phyto by zooplankton
zw3d(A2D(0),1:jpkm1) = zgrazing(A2D(0),1:jpkm1) * 1.e+3 * rfact2r * tmask(A2D(0),1:jpkm1)
CALL iom_put( "GRAZ1", zw3d )
ENDIF
!
IF( iom_use ( "FEZOO" ) ) THEN
zw3d(A2D(0),1:jpkm1) = zfezoo(A2D(0),1:jpkm1) * 1e9 * 1.e+3 * rfact2r * tmask(A2D(0),1:jpkm1)
CALL iom_put( "FEZOO", zw3d )
ENDIF
!
IF( ln_ligand .AND. iom_use( "LPRODZ" ) ) THEN
zw3d(A2D(0),1:jpkm1) = (zgrarem - zgrarsig) * 1.e+3 * rfact2r * tmask(A2D(0),1:jpkm1)
CALL iom_put( "LPRODZ", zw3d * ldocz * 1e9 )
ENDIF
!
DEALLOCATE( zw3d )
ENDIF
!
IF(sn_cfctl%l_prttrc) THEN ! print mean trends (used for debugging)
......
......@@ -74,7 +74,7 @@ CONTAINS
IF( ln_timing ) CALL timing_start('p4z_mort_nano')
!
prodcal(:,:,:) = 0._wp ! calcite production variable set to zero
DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1)
DO_3D( 0, 0, 0, 0, 1, jpkm1)
zcompaph = MAX( ( tr(ji,jj,jk,jpphy,Kbb) - 1e-9 ), 0.e0 )
! Quadratic mortality of nano due to aggregation during
......@@ -152,7 +152,7 @@ CONTAINS
! This is due to the production of EPS by stressed cells
! -------------------------------------------------------------
DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1)
DO_3D( 0, 0, 0, 0, 1, jpkm1)
zcompadi = MAX( ( tr(ji,jj,jk,jpdia,Kbb) - 1e-9), 0. )
......
......@@ -63,10 +63,12 @@ CONTAINS
INTEGER :: irgb
REAL(wp) :: zchl
REAL(wp) :: zc0 , zc1 , zc2, zc3, z1_dep
REAL(wp), ALLOCATABLE, DIMENSION(:,: ) :: zetmp5
REAL(wp), DIMENSION(jpi,jpj ) :: zdepmoy, zetmp1, zetmp2, zetmp3, zetmp4
REAL(wp), DIMENSION(jpi,jpj ) :: zqsr100, zqsr_corr
REAL(wp), DIMENSION(jpi,jpj,jpk) :: zpar, ze0, ze1, ze2, ze3
REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zetmp5
REAL(wp), DIMENSION(A2D(0) ) :: zdepmoy, zetmp1, zetmp2, zetmp3, zetmp4
REAL(wp), DIMENSION(A2D(0) ) :: zqsr100, zqsr_corr
REAL(wp), DIMENSION(A2D(0),jpk) :: zpar, ze0, ze1, ze2, ze3
REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: zw3d
REAL(wp), ALLOCATABLE, DIMENSION(:,: ) :: zw2d
!!---------------------------------------------------------------------
!
IF( ln_timing ) CALL timing_start('p4z_opt')
......@@ -88,7 +90,7 @@ CONTAINS
!
! Computation of the light attenuation parameters based on a
! look-up table
DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1)
DO_3D( 0, 0, 0, 0, 1, jpkm1)
zchl = ( tr(ji,jj,jk,jpnch,Kbb) + tr(ji,jj,jk,jpdch,Kbb) + rtrn ) * 1.e6
IF( ln_p5z ) zchl = zchl + tr(ji,jj,jk,jppch,Kbb) * 1.e6
zchl = MIN( 10. , MAX( 0.05, zchl ) )
......@@ -116,36 +118,40 @@ CONTAINS
! not fully correct with LIM3 and SI3 but no information is
! currently available to do a better job. SHould be improved in the
! (near) future.
zqsr_corr(:,:) = qsr_mean(:,:) / ( 1.-fr_i(:,:) + rtrn )
DO_2D( 0, 0, 0, 0 )
zqsr_corr(ji,jj) = qsr_mean(ji,jj) / ( 1.-fr_i(ji,jj) + rtrn )
END_2D
!
CALL p4z_opt_par( kt, Kmm, zqsr_corr, ze1, ze2, ze3, pqsr100 = zqsr100 )
!
! Used PAR is computed for each phytoplankton species
! etot_ndcy is PAR at level jk averaged over 24h.
! Due to their size, they have different light absorption characteristics
DO jk = 1, nksr
etot_ndcy(:,:,jk) = ze1(:,:,jk) + ze2(:,:,jk) + ze3(:,:,jk)
END DO
DO_3D( 0, 0, 0, 0, 1, nksr )
etot_ndcy(ji,jj,jk) = ze1(ji,jj,jk) + ze2(ji,jj,jk) + ze3(ji,jj,jk)
END_3D
!
! SW over the ice free zone of the grid cell. This assumes that
! SW is zero below sea ice which is a very crude assumption that is
! not fully correct with LIM3 and SI3 but no information is
! currently available to do a better job. SHould be improved in the
! (near) future.
zqsr_corr(:,:) = qsr(:,:) / ( 1.-fr_i(:,:) + rtrn )
DO_2D( 0, 0, 0, 0 )
zqsr_corr(ji,jj) = qsr(ji,jj) / ( 1.-fr_i(ji,jj) + rtrn )
END_2D
!
CALL p4z_opt_par( kt, Kmm, zqsr_corr, ze1, ze2, ze3 )
!
! Total PAR computation at level jk that includes the diurnal cycle
DO jk = 1, nksr
etot (:,:,jk) = ze1(:,:,jk) + ze2(:,:,jk) + ze3(:,:,jk)
enano(:,:,jk) = 1.85 * ze1(:,:,jk) + 0.69 * ze2(:,:,jk) + 0.46 * ze3(:,:,jk)
ediat(:,:,jk) = 1.62 * ze1(:,:,jk) + 0.74 * ze2(:,:,jk) + 0.63 * ze3(:,:,jk)
END DO
DO_3D( 0, 0, 0, 0, 1, nksr )
etot (ji,jj,jk) = ze1(ji,jj,jk) + ze2(ji,jj,jk) + ze3(ji,jj,jk)
enano(ji,jj,jk) = 1.85 * ze1(ji,jj,jk) + 0.69 * ze2(ji,jj,jk) + 0.46 * ze3(ji,jj,jk)
ediat(ji,jj,jk) = 1.62 * ze1(ji,jj,jk) + 0.74 * ze2(ji,jj,jk) + 0.63 * ze3(ji,jj,jk)
END_3D
IF( ln_p5z ) THEN
DO jk = 1, nksr
epico (:,:,jk) = 1.94 * ze1(:,:,jk) + 0.66 * ze2(:,:,jk) + 0.4 * ze3(:,:,jk)
END DO
DO_3D( 0, 0, 0, 0, 1, nksr )
epico(ji,jj,jk) = 1.94 * ze1(ji,jj,jk) + 0.66 * ze2(ji,jj,jk) + 0.4 * ze3(ji,jj,jk)
END_3D
ENDIF
ELSE ! No diurnal cycle in PISCES
......@@ -157,22 +163,24 @@ CONTAINS
! not fully correct with LIM3 and SI3 but no information is
! currently available to do a better job. SHould be improved in the
! (near) future.
zqsr_corr(:,:) = qsr_mean(:,:) / ( 1.-fr_i(:,:) + rtrn )
DO_2D( 0, 0, 0, 0 )
zqsr_corr(ji,jj) = qsr_mean(ji,jj) / ( 1.-fr_i(ji,jj) + rtrn )
END_2D
!
CALL p4z_opt_par( kt, Kmm, zqsr_corr, ze1, ze2, ze3, pqsr100 = zqsr100 )
!
! Used PAR is computed for each phytoplankton species
! etot_ndcy is PAR at level jk averaged over 24h.
! Due to their size, they have different light absorption characteristics
DO jk = 1, nksr
etot_ndcy(:,:,jk) = ze1(:,:,jk) + ze2(:,:,jk) + ze3(:,:,jk)
enano (:,:,jk) = 1.85 * ze1(:,:,jk) + 0.69 * ze2(:,:,jk) + 0.46 * ze3(:,:,jk)
ediat (:,:,jk) = 1.62 * ze1(:,:,jk) + 0.74 * ze2(:,:,jk) + 0.63 * ze3(:,:,jk)
END DO
DO_3D( 0, 0, 0, 0, 1, nksr )
etot_ndcy(ji,jj,jk) = ze1(ji,jj,jk) + ze2(ji,jj,jk) + ze3(ji,jj,jk)
enano (ji,jj,jk) = 1.85 * ze1(ji,jj,jk) + 0.69 * ze2(ji,jj,jk) + 0.46 * ze3(ji,jj,jk)
ediat (ji,jj,jk) = 1.62 * ze1(ji,jj,jk) + 0.74 * ze2(ji,jj,jk) + 0.63 * ze3(ji,jj,jk)
END_3D
IF( ln_p5z ) THEN
DO jk = 1, nksr
epico (:,:,jk) = 1.94 * ze1(:,:,jk) + 0.66 * ze2(:,:,jk) + 0.4 * ze3(:,:,jk)
END DO
DO_3D( 0, 0, 0, 0, 1, nksr )
epico(ji,jj,jk) = 1.94 * ze1(ji,jj,jk) + 0.66 * ze2(ji,jj,jk) + 0.4 * ze3(ji,jj,jk)
END_3D
ENDIF
!
! SW over the ice free zone of the grid cell. This assumes that
......@@ -180,14 +188,16 @@ CONTAINS
! not fully correct with LIM3 and SI3 but no information is
! currently available to do a better job. SHould be improved in the
! (near) future.
zqsr_corr(:,:) = qsr(:,:) / ( 1.-fr_i(:,:) + rtrn )
DO_2D( 0, 0, 0, 0 )
zqsr_corr(ji,jj) = qsr(ji,jj) / ( 1.-fr_i(ji,jj) + rtrn )
END_2D
!
CALL p4z_opt_par( kt, Kmm, zqsr_corr, ze1, ze2, ze3 )
!
! Total PAR computation at level jk that includes the diurnal cycle
DO jk = 1, nksr
etot(:,:,jk) = ze1(:,:,jk) + ze2(:,:,jk) + ze3(:,:,jk)
END DO
DO_3D( 0, 0, 0, 0, 1, nksr )
etot(ji,jj,jk) = ze1(ji,jj,jk) + ze2(ji,jj,jk) + ze3(ji,jj,jk)
END_3D
ENDIF
!
ELSE ! no diurnal cycle
......@@ -198,22 +208,24 @@ CONTAINS
! not fully correct with LIM3 and SI3 but no information is
! currently available to do a better job. SHould be improved in the
! (near) future.
zqsr_corr(:,:) = qsr(:,:) / ( 1.-fr_i(:,:) + rtrn )
DO_2D( 0, 0, 0, 0 )
zqsr_corr(ji,jj) = qsr(ji,jj) / ( 1.-fr_i(ji,jj) + rtrn )
END_2D
!
CALL p4z_opt_par( kt, Kmm, zqsr_corr, ze1, ze2, ze3, pqsr100 = zqsr100 )
!
! Used PAR is computed for each phytoplankton species
! Due to their size, they have different light absorption characteristics
DO jk = 1, nksr
etot (:,:,jk) = ze1(:,:,jk) + ze2(:,:,jk) + ze3(:,:,jk) ! Total PAR
enano(:,:,jk) = 1.85 * ze1(:,:,jk) + 0.69 * ze2(:,:,jk) + 0.46 * ze3(:,:,jk) ! Nanophytoplankton
ediat(:,:,jk) = 1.62 * ze1(:,:,jk) + 0.74 * ze2(:,:,jk) + 0.63 * ze3(:,:,jk) ! Diatoms
END DO
DO_3D( 0, 0, 0, 0, 1, nksr )
etot (ji,jj,jk) = ze1(ji,jj,jk) + ze2(ji,jj,jk) + ze3(ji,jj,jk)
enano(ji,jj,jk) = 1.85 * ze1(ji,jj,jk) + 0.69 * ze2(ji,jj,jk) + 0.46 * ze3(ji,jj,jk)
ediat(ji,jj,jk) = 1.62 * ze1(ji,jj,jk) + 0.74 * ze2(ji,jj,jk) + 0.63 * ze3(ji,jj,jk)
END_3D
IF( ln_p5z ) THEN
DO jk = 1, nksr
epico(:,:,jk) = 1.94 * ze1(:,:,jk) + 0.66 * ze2(:,:,jk) + 0.4 * ze3(:,:,jk) ! Picophytoplankton (PISCES-QUOTA)
END DO
DO_3D( 0, 0, 0, 0, 1, nksr )
epico(ji,jj,jk) = 1.94 * ze1(ji,jj,jk) + 0.66 * ze2(ji,jj,jk) + 0.4 * ze3(ji,jj,jk) ! Picophytoplankton (PISCES-QUOTA)
END_3D
ENDIF
etot_ndcy(:,:,:) = etot(:,:,:)
ENDIF
......@@ -224,10 +236,12 @@ CONTAINS
! ! ------------------------
CALL p4z_opt_par( kt, Kmm, qsr, ze1, ze2, ze3, pe0=ze0 )
!
etot3(:,:,1) = qsr(:,:) * tmask(:,:,1)
DO jk = 2, nksr + 1
etot3(:,:,jk) = ( ze0(:,:,jk) + ze1(:,:,jk) + ze2(:,:,jk) + ze3(:,:,jk) ) * tmask(:,:,jk)
END DO
DO_2D( 0, 0, 0, 0 )
etot3(ji,jj,1) = qsr(ji,jj) * tmask(ji,jj,1)
END_2D
DO_3D( 0, 0, 0, 0, 2, nksr+1 )
etot3(ji,jj,jk) = ( ze0(ji,jj,jk) + ze1(ji,jj,jk) + ze2(ji,jj,jk) + ze3(ji,jj,jk) ) * tmask(ji,jj,jk)
END_3D
! ! ------------------------
ENDIF
......@@ -236,11 +250,13 @@ CONTAINS
! (1) The classical definition based on the relative threshold value
! (2) An alternative definition based on a absolute threshold value.
! -------------------------------------------------------------------
neln(:,:) = 1
heup (:,:) = gdepw(:,:,2,Kmm)
heup_01(:,:) = gdepw(:,:,2,Kmm)
DO_2D( 0, 0, 0, 0 )
neln (ji,jj) = 1
heup (ji,jj) = gdepw(ji,jj,2,Kmm)
heup_01(ji,jj) = gdepw(ji,jj,2,Kmm)
END_2D
DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 2, nksr)
DO_3D( 0, 0, 0, 0, 2, nksr)
IF( etot_ndcy(ji,jj,jk) * tmask(ji,jj,jk) >= zqsr100(ji,jj) ) THEN
neln(ji,jj) = jk+1 ! Euphotic level : 1rst T-level strictly below Euphotic layer
! ! nb: ensure the compatibility with nmld_trc definition in trd_mld_trc_zint
......@@ -261,7 +277,7 @@ CONTAINS
zetmp1 (:,:) = 0.e0
zetmp2 (:,:) = 0.e0
DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, nksr)
DO_3D( 0, 0, 0, 0, 1, nksr)
IF( gdepw(ji,jj,jk+1,Kmm) <= hmld(ji,jj) ) THEN
zetmp1 (ji,jj) = zetmp1 (ji,jj) + etot (ji,jj,jk) * e3t(ji,jj,jk,Kmm) ! Actual PAR for remineralisation
zetmp2 (ji,jj) = zetmp2 (ji,jj) + etot_ndcy(ji,jj,jk) * e3t(ji,jj,jk,Kmm) ! Par averaged over 24h for production
......@@ -272,7 +288,7 @@ CONTAINS
emoy(:,:,:) = etot(:,:,:) ! remineralisation
zpar(:,:,:) = etot_ndcy(:,:,:) ! diagnostic : PAR with no diurnal cycle
!
DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, nksr)
DO_3D( 0, 0, 0, 0, 1, nksr)
IF( gdepw(ji,jj,jk+1,Kmm) <= hmld(ji,jj) ) THEN
z1_dep = 1. / ( zdepmoy(ji,jj) + rtrn )
emoy (ji,jj,jk) = zetmp1(ji,jj) * z1_dep
......@@ -286,7 +302,7 @@ CONTAINS
zetmp3 (:,:) = 0.e0
zetmp4 (:,:) = 0.e0
!
DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, nksr)
DO_3D( 0, 0, 0, 0, 1, nksr)
IF( gdepw(ji,jj,jk+1,Kmm) <= MIN(hmld(ji,jj), heup_01(ji,jj)) ) THEN
zetmp3 (ji,jj) = zetmp3 (ji,jj) + enano (ji,jj,jk) * e3t(ji,jj,jk,Kmm) ! Nanophytoplankton
zetmp4 (ji,jj) = zetmp4 (ji,jj) + ediat (ji,jj,jk) * e3t(ji,jj,jk,Kmm) ! Diatoms
......@@ -296,7 +312,7 @@ CONTAINS
enanom(:,:,:) = enano(:,:,:)
ediatm(:,:,:) = ediat(:,:,:)
!
DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, nksr)
DO_3D( 0, 0, 0, 0, 1, nksr)
IF( gdepw(ji,jj,jk+1,Kmm) <= hmld(ji,jj) ) THEN
z1_dep = 1. / ( zdepmoy(ji,jj) + rtrn )
enanom(ji,jj,jk) = zetmp3(ji,jj) * z1_dep
......@@ -306,8 +322,8 @@ CONTAINS
!
IF( ln_p5z ) THEN
! Picophytoplankton when using PISCES-QUOTA
ALLOCATE( zetmp5(jpi,jpj) ) ; zetmp5 (:,:) = 0.e0
DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, nksr)
ALLOCATE( zetmp5(A2D(0)) ) ; zetmp5 (:,:) = 0.e0
DO_3D( 0, 0, 0, 0, 1, nksr)
IF( gdepw(ji,jj,jk+1,Kmm) <= MIN(hmld(ji,jj), heup_01(ji,jj)) ) THEN
zetmp5(ji,jj) = zetmp5 (ji,jj) + epico(ji,jj,jk) * e3t(ji,jj,jk,Kmm)
ENDIF
......@@ -315,7 +331,7 @@ CONTAINS
!
epicom(:,:,:) = epico(:,:,:)
!
DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, nksr)
DO_3D( 0, 0, 0, 0, 1, nksr)
IF( gdepw(ji,jj,jk+1,Kmm) <= hmld(ji,jj) ) THEN
z1_dep = 1. / ( zdepmoy(ji,jj) + rtrn )
epicom(ji,jj,jk) = zetmp5(ji,jj) * z1_dep
......@@ -325,10 +341,18 @@ CONTAINS
ENDIF
!
IF( lk_iomput .AND. knt == nrdttrc ) THEN
CALL iom_put( "Heup" , heup(:,: ) * tmask(:,:,1) ) ! euphotic layer deptht
IF( iom_use( "Heup" ) ) THEN
ALLOCATE( zw2d(A2D(0)) )
zw2d(A2D(0)) = heup(A2D(0)) * tmask(A2D(0),1)
CALL iom_put( "Heup", zw2d ) ! Euphotic layer depth
DEALLOCATE( zw2d )
ENDIF
IF( iom_use( "PAR" ) ) THEN
zpar(:,:,1) = zpar(:,:,1) * ( 1._wp - fr_i(:,:) )
CALL iom_put( "PAR", zpar(:,:,:) * tmask(:,:,:) ) ! Photosynthetically Available Radiation
ALLOCATE( zw3d(A2D(0),jpk) ) ; zw3d(A2D(0),jpk) = 0._wp
zw3d(A2D(0),1) = zpar(A2D(0),1) * ( 1._wp - fr_i(A2D(0)) ) * tmask(A2D(0),1)
zw3d(A2D(0),2:jpkm1) = zpar(A2D(0),2:jpkm1) * tmask(A2D(0),2:jpkm1)
CALL iom_put( "PAR", zw3d ) ! Photosynthetically Available Radiation
DEALLOCATE( zw3d )
ENDIF
ENDIF
!
......@@ -345,15 +369,15 @@ CONTAINS
!! for a given shortwave radiation
!!
!!----------------------------------------------------------------------
INTEGER , INTENT(in) :: kt ! ocean time-step
INTEGER , INTENT(in) :: Kmm ! ocean time-index
REAL(wp), DIMENSION(jpi,jpj) , INTENT(in ) :: pqsr ! shortwave
REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pe1 , pe2 , pe3 ! PAR ( R-G-B)
REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout), OPTIONAL :: pe0 !
REAL(wp), DIMENSION(jpi,jpj) , INTENT( out), OPTIONAL :: pqsr100 !
INTEGER , INTENT(in) :: kt ! ocean time-step
INTEGER , INTENT(in) :: Kmm ! ocean time-index
REAL(wp), DIMENSION(A2D(0)) , INTENT(in ) :: pqsr ! shortwave
REAL(wp), DIMENSION(A2D(0),jpk), INTENT(inout) :: pe1 , pe2 , pe3 ! PAR ( R-G-B)
REAL(wp), DIMENSION(A2D(0),jpk), INTENT(inout), OPTIONAL :: pe0 !
REAL(wp), DIMENSION(A2D(0)) , INTENT( out), OPTIONAL :: pqsr100 !
!
INTEGER :: ji, jj, jk ! dummy loop indices
REAL(wp), DIMENSION(jpi,jpj) :: zqsr ! shortwave
REAL(wp), DIMENSION(A2D(0)) :: zqsr ! shortwave
!!----------------------------------------------------------------------
! Real shortwave
......@@ -371,7 +395,7 @@ CONTAINS
pe2(:,:,1) = zqsr(:,:)
pe3(:,:,1) = zqsr(:,:)
!
DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 2, nksr + 1)
DO_3D( 0, 0, 0, 0, 2, nksr + 1)
pe0(ji,jj,jk) = pe0(ji,jj,jk-1) * EXP( -e3t(ji,jj,jk-1,Kmm) * xsi0r )
pe1(ji,jj,jk) = pe1(ji,jj,jk-1) * EXP( -ekb (ji,jj,jk-1 ) )
pe2(ji,jj,jk) = pe2(ji,jj,jk-1) * EXP( -ekg (ji,jj,jk-1 ) )
......@@ -384,7 +408,7 @@ CONTAINS
pe2(:,:,1) = zqsr(:,:) * EXP( -0.5 * ekg(:,:,1) )
pe3(:,:,1) = zqsr(:,:) * EXP( -0.5 * ekr(:,:,1) )
!
DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 2, nksr)
DO_3D( 0, 0, 0, 0, 2, nksr)
pe1(ji,jj,jk) = pe1(ji,jj,jk-1) * EXP( -0.5 * ( ekb(ji,jj,jk-1) + ekb(ji,jj,jk) ) )
pe2(ji,jj,jk) = pe2(ji,jj,jk-1) * EXP( -0.5 * ( ekg(ji,jj,jk-1) + ekg(ji,jj,jk) ) )
pe3(ji,jj,jk) = pe3(ji,jj,jk-1) * EXP( -0.5 * ( ekr(ji,jj,jk-1) + ekr(ji,jj,jk) ) )
......@@ -419,7 +443,9 @@ CONTAINS
IF( ln_varpar ) THEN
IF( kt == nit000 .OR. ( kt /= nit000 .AND. ntimes_par > 1 ) ) THEN
CALL fld_read( kt, 1, sf_par )
par_varsw(:,:) = ( sf_par(1)%fnow(:,:,1) ) / 3.0
DO_2D( 0, 0, 0, 0 )
par_varsw(ji,jj) = ( sf_par(1)%fnow(ji,jj,1) ) / 3.0
END_2D
ENDIF
ENDIF
!
......@@ -479,7 +505,7 @@ CONTAINS
IF(lwp) WRITE(numout,*)
IF(lwp) WRITE(numout,*) ' ==>>> initialize variable par fraction (ln_varpar=T)'
!
ALLOCATE( par_varsw(jpi,jpj) )
ALLOCATE( par_varsw(A2D(0)) )
!
ALLOCATE( sf_par(1), STAT=ierr ) !* allocate and fill sf_sst (forcing structure) with sn_sst
IF( ierr > 0 ) CALL ctl_stop( 'STOP', 'p4z_opt_init: unable to allocate sf_par structure' )
......@@ -510,8 +536,7 @@ CONTAINS
!! *** ROUTINE p4z_opt_alloc ***
!!----------------------------------------------------------------------
!
ALLOCATE( ekb(jpi,jpj,jpk), ekr(jpi,jpj,jpk), &
ekg(jpi,jpj,jpk), STAT= p4z_opt_alloc )
ALLOCATE( ekb(A2D(0),jpk), ekr(A2D(0),jpk), ekg(A2D(0),jpk), STAT= p4z_opt_alloc )
!
IF( p4z_opt_alloc /= 0 ) CALL ctl_stop( 'STOP', 'p4z_opt_alloc : failed to allocate arrays.' )
!
......
......@@ -74,9 +74,10 @@ CONTAINS
REAL(wp) :: zofer2, zofer3
REAL(wp) :: zrfact2
CHARACTER (len=25) :: charout
REAL(wp), DIMENSION(jpi,jpj ) :: totprod, totthick, totcons
REAL(wp), DIMENSION(jpi,jpj,jpk) :: zremipoc, zremigoc, zorem3, ztremint, zfolimi
REAL(wp), DIMENSION(jpi,jpj,jpk,jcpoc) :: alphag
REAL(wp), DIMENSION(A2D(0) ) :: totprod, totthick, totcons
REAL(wp), DIMENSION(A2D(0),jpk) :: zremipoc, zremigoc, zorem3, ztremint, zfolimi
REAL(wp), DIMENSION(A2D(0),jpk,jcpoc) :: alphag
REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: zw3d
!!---------------------------------------------------------------------
!
IF( ln_timing ) CALL timing_start('p4z_poc')
......@@ -118,7 +119,7 @@ CONTAINS
! a standard parameterisation with a constant lability
! -----------------------------------------------------------------------
ztremint(:,:,:) = zremigoc(:,:,:)
DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1)
DO_3D( 0, 0, 0, 0, 1, jpkm1)
IF (tmask(ji,jj,jk) == 1.) THEN
zdep = hmld(ji,jj)
!
......@@ -204,7 +205,7 @@ CONTAINS
IF( ln_p4z ) THEN
! The standard PISCES part
DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1)
DO_3D( 0, 0, 0, 0, 1, jpkm1)
! POC degradation by bacterial activity. It is a function
! of the mean lability and of temperature. This also includes
! shrinking of particles due to the bacterial activity
......@@ -226,7 +227,7 @@ CONTAINS
zfolimi(ji,jj,jk) = zofer2
END_3D
ELSE
DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1)
DO_3D( 0, 0, 0, 0, 1, jpkm1)
! POC degradation by bacterial activity. It is a function
! of the mean lability and of temperature. This also includes
! shrinking of particles due to the bacterial activity
......@@ -274,7 +275,7 @@ CONTAINS
! intregrated production and consumption of POC in the mixed layer
! ----------------------------------------------------------------
DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1)
DO_3D( 0, 0, 0, 0, 1, jpkm1)
zdep = hmld(ji,jj)
IF (tmask(ji,jj,jk) == 1. .AND. gdept(ji,jj,jk,Kmm) <= zdep ) THEN
totprod(ji,jj) = totprod(ji,jj) + prodpoc(ji,jj,jk) * e3t(ji,jj,jk,Kmm) * rday/ rfact2
......@@ -290,7 +291,7 @@ CONTAINS
! mixing.
! ---------------------------------------------------------------------
ztremint(:,:,:) = zremipoc(:,:,:)
DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1)
DO_3D( 0, 0, 0, 0, 1, jpkm1)
IF (tmask(ji,jj,jk) == 1.) THEN
zdep = hmld(ji,jj)
alphat = 0.0
......@@ -323,7 +324,7 @@ CONTAINS
! that since we need the lability spectrum of GOC, GOC spectrum
! should be determined before.
! -----------------------------------------------------------------------
DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 2, jpkm1)
DO_3D( 0, 0, 0, 0, 2, jpkm1)
IF (tmask(ji,jj,jk) == 1.) THEN
zdep = hmld(ji,jj)
IF( gdept(ji,jj,jk,Kmm) > zdep ) THEN
......@@ -397,7 +398,7 @@ CONTAINS
ENDIF
IF( ln_p4z ) THEN
DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1)
DO_3D( 0, 0, 0, 0, 1, jpkm1)
IF (tmask(ji,jj,jk) == 1.) THEN
! POC disaggregation by turbulence and bacterial activity.It is a function
! of the mean lability and of temperature
......@@ -416,7 +417,7 @@ CONTAINS
ENDIF
END_3D
ELSE
DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1)
DO_3D( 0, 0, 0, 0, 1, jpkm1)
! POC disaggregation by turbulence and bacterial activity.It is a function
! of the mean lability and of temperature
!--------------------------------------------------------
......@@ -443,9 +444,25 @@ CONTAINS
IF( lk_iomput ) THEN
IF( knt == nrdttrc ) THEN
zrfact2 = 1.e3 * rfact2r
CALL iom_put( "REMINP" , zremipoc(:,:,:) * tmask(:,:,:) ) ! Remineralisation rate of small particles
CALL iom_put( "REMING" , zremigoc(:,:,:) * tmask(:,:,:) ) ! Remineralisation rate of large particles
CALL iom_put( "REMINF" , zfolimi(:,:,:) * tmask(:,:,:) * 1.e+9 * zrfact2 ) ! Remineralisation of biogenic particulate iron
!
ALLOCATE( zw3d(A2D(0),jpk) ) ; zw3d(A2D(0),jpk) = 0._wp
!
IF( iom_use ( "REMINP" ) ) THEN ! Remineralisation rate of small particles
zw3d(A2D(0),1:jpkm1) = zremipoc(A2D(0),1:jpkm1) * tmask(A2D(0),1:jpkm1)
CALL iom_put( "REMINP", zw3d )
ENDIF
!
IF( iom_use ( "REMING" ) ) THEN ! Remineralisation rate of large particles
zw3d(A2D(0),1:jpkm1) = zremigoc(A2D(0),1:jpkm1) * tmask(A2D(0),1:jpkm1)
CALL iom_put( "REMING", zw3d )
ENDIF
!
IF( iom_use ( "REMINF" ) ) THEN ! Remineralisation of biogenic particulate iron
zw3d(A2D(0),1:jpkm1) = zfolimi(A2D(0),1:jpkm1) * tmask(A2D(0),1:jpkm1) * 1.e+9 * zrfact2
CALL iom_put( "REMINF", zw3d )
ENDIF
!
DEALLOCATE ( zw3d )
ENDIF
ENDIF
......@@ -508,7 +525,7 @@ CONTAINS
! Discretization along the lability space
! ---------------------------------------
!
ALLOCATE( alphan(jcpoc) , reminp(jcpoc) , alphap(jpi,jpj,jpk,jcpoc) )
ALLOCATE( alphan(jcpoc) , reminp(jcpoc) , alphap(A2D(0),jpk,jcpoc) )
!
IF (jcpoc > 1) THEN ! Case when more than one lability class is used
!
......
......@@ -77,13 +77,13 @@ CONTAINS
REAL(wp) :: zratiosi, zmaxsi, zlimfac, zsizetmp, zfecnm, zfecdm
REAL(wp) :: zprod, zval
CHARACTER (len=25) :: charout
REAL(wp), DIMENSION(jpi,jpj,jpk) :: zprmaxn,zprmaxd
REAL(wp), DIMENSION(jpi,jpj,jpk) :: zpislopeadn, zpislopeadd, zysopt
REAL(wp), DIMENSION(jpi,jpj,jpk) :: zprdia, zprbio, zprchld, zprchln
REAL(wp), DIMENSION(jpi,jpj,jpk) :: zprorcan, zprorcad, zprofed, zprofen
REAL(wp), DIMENSION(jpi,jpj,jpk) :: zpronewn, zpronewd
REAL(wp), DIMENSION(jpi,jpj,jpk) :: zmxl_fac, zmxl_chl
REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: zpligprod
REAL(wp), DIMENSION(A2D(0),jpk) :: zprmaxn,zprmaxd
REAL(wp), DIMENSION(A2D(0),jpk) :: zpislopeadn, zpislopeadd, zysopt
REAL(wp), DIMENSION(A2D(0),jpk) :: zprdia, zprbio, zprchld, zprchln
REAL(wp), DIMENSION(A2D(0),jpk) :: zprorcan, zprorcad, zprofed, zprofen
REAL(wp), DIMENSION(A2D(0),jpk) :: zpronewn, zpronewd
REAL(wp), DIMENSION(A2D(0),jpk) :: zmxl_fac, zmxl_chl
REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: zw3d
!!---------------------------------------------------------------------
!
IF( ln_timing ) CALL timing_start('p4z_prod')
......@@ -110,7 +110,7 @@ CONTAINS
! -------------------------------------------------------------------------
IF ( ln_p4z_dcyc ) THEN ! Diurnal cycle in PISCES
DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1)
DO_3D( 0, 0, 0, 0, 1, jpkm1)
IF( etot_ndcy(ji,jj,jk) > 1.E-3 ) THEN
IF( gdepw(ji,jj,jk+1,Kmm) <= hmld(ji,jj) ) THEN
zval = MIN(1., heup_01(ji,jj) / ( hmld(ji,jj) + rtrn ))
......@@ -121,7 +121,7 @@ CONTAINS
END_3D
ELSE ! No diurnal cycle in PISCES
DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1)
DO_3D( 0, 0, 0, 0, 1, jpkm1)
IF( etot_ndcy(ji,jj,jk) > 1.E-3 ) THEN
zval = MAX( 1., strn(ji,jj) )
IF( gdepw(ji,jj,jk+1,Kmm) <= hmld(ji,jj) ) THEN
......@@ -141,7 +141,7 @@ CONTAINS
! to exclude the effect of nutrient limitation and temperature in the PI
! curve following Vichi et al. (2007)
! -----------------------------------------------------------------------
DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1)
DO_3D( 0, 0, 0, 0, 1, jpkm1)
IF( etot_ndcy(ji,jj,jk) > 1.E-3 ) THEN
ztn = MAX( 0., ts(ji,jj,jk,jp_tem,Kmm) - 15. )
zadap = xadap * ztn / ( 2.+ ztn )
......@@ -163,7 +163,7 @@ CONTAINS
ENDIF
END_3D
DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1)
DO_3D( 0, 0, 0, 0, 1, jpkm1)
IF( etot_ndcy(ji,jj,jk) > 1.E-3 ) THEN
! Computation of production function for Carbon
! Actual light levels are used here
......@@ -190,7 +190,7 @@ CONTAINS
! Computation of a proxy of the N/C quota from nutrient limitation
! and light limitation. Steady state is assumed to allow the computation
! ----------------------------------------------------------------------
DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1)
DO_3D( 0, 0, 0, 0, 1, jpkm1)
zval = MIN( xnanopo4(ji,jj,jk), ( xnanonh4(ji,jj,jk) + xnanono3(ji,jj,jk) ) ) &
& * zprmaxn(ji,jj,jk) / ( zprbio(ji,jj,jk) + rtrn )
quotan(ji,jj,jk) = MIN( 1., 0.3 + 0.7 * zval )
......@@ -200,7 +200,7 @@ CONTAINS
END_3D
DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1)
DO_3D( 0, 0, 0, 0, 1, jpkm1)
IF( etot_ndcy(ji,jj,jk) > 1.E-3 ) THEN
! Si/C of diatoms
......@@ -234,14 +234,14 @@ CONTAINS
! Sea-ice effect on production
! No production is assumed below sea ice
! --------------------------------------
DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1)
DO_3D( 0, 0, 0, 0, 1, jpkm1)
zprbio(ji,jj,jk) = zprbio(ji,jj,jk) * ( 1. - fr_i(ji,jj) )
zprdia(ji,jj,jk) = zprdia(ji,jj,jk) * ( 1. - fr_i(ji,jj) )
END_3D
! Computation of the various production and nutrient uptake terms
! ---------------------------------------------------------------
DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1)
DO_3D( 0, 0, 0, 0, 1, jpkm1)
IF( etot_ndcy(ji,jj,jk) > 1.E-3 ) THEN
! production terms for nanophyto. (C)
zprorcan(ji,jj,jk) = zprbio(ji,jj,jk) * xlimphy(ji,jj,jk) * tr(ji,jj,jk,jpphy,Kbb) * rfact2
......@@ -303,7 +303,7 @@ CONTAINS
! Computation of the chlorophyll production terms
! The parameterization is taken from Geider et al. (1997)
! -------------------------------------------------------
DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1)
DO_3D( 0, 0, 0, 0, 1, jpkm1)
IF( etot_ndcy(ji,jj,jk) > 1.E-3 ) THEN
! production terms for nanophyto. ( chlorophyll )
znanotot = enanom(ji,jj,jk) / ( zmxl_chl(ji,jj,jk) + rtrn )
......@@ -326,7 +326,7 @@ CONTAINS
END_3D
! Update the arrays TRA which contain the biological sources and sinks
DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1)
DO_3D( 0, 0, 0, 0, 1, jpkm1)
IF( etot_ndcy(ji,jj,jk) > 1.E-3 ) THEN
zpptot = zprorcan(ji,jj,jk) + zprorcad(ji,jj,jk)
zpnewtot = zpronewn(ji,jj,jk) + zpronewd(ji,jj,jk)
......@@ -362,7 +362,7 @@ CONTAINS
! Shaked et al. (2020)
! -------------------------------------------------------------------------
IF( ln_ligand ) THEN
DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1)
DO_3D( 0, 0, 0, 0, 1, jpkm1)
IF( etot_ndcy(ji,jj,jk) > 1.E-3 ) THEN
zproddoc = excretd * zprorcad(ji,jj,jk) + excretn * zprorcan(ji,jj,jk)
zprodfer = texcretn * zprofen(ji,jj,jk) + texcretd * zprofed(ji,jj,jk)
......@@ -376,38 +376,108 @@ CONTAINS
! Output of the diagnostics
! Total primary production per year
IF( iom_use( "tintpp" ) .OR. ( ln_check_mass .AND. kt == nitend .AND. knt == nrdttrc ) ) &
& tpp = glob_sum( 'p4zprod', ( zprorcan(:,:,:) + zprorcad(:,:,:) ) * cvol(:,:,:) )
IF( iom_use( "tintpp" ) .OR. ( ln_check_mass .AND. kt == nitend .AND. knt == nrdttrc ) ) THEN
ALLOCATE( zw3d(A2D(0),jpk) ) ; zw3d(A2D(0),jpk) = 0._wp
zw3d(A2D(0),1:jpkm1) = ( zprorcan(A2D(0),1:jpkm1) + zprorcad(A2D(0),1:jpkm1) ) * cvol(A2D(0),1:jpkm1)
tpp = glob_sum( 'p4zprod', zw3d )
DEALLOCATE ( zw3d )
ENDIF
IF( lk_iomput .AND. knt == nrdttrc ) THEN
zfact = 1.e+3 * rfact2r ! conversion from mol/l/kt to mol/m3/s
!
CALL iom_put( "PPPHYN" , zprorcan(:,:,:) * zfact * tmask(:,:,:) ) ! primary production by nanophyto
CALL iom_put( "PPPHYD" , zprorcad(:,:,:) * zfact * tmask(:,:,:) ) ! primary production by diatomes
CALL iom_put( "PPNEWN" , zpronewn(:,:,:) * zfact * tmask(:,:,:) ) ! new primary production by nanophyto
CALL iom_put( "PPNEWD" , zpronewd(:,:,:) * zfact * tmask(:,:,:) ) ! new primary production by diatomes
CALL iom_put( "PBSi" , zprorcad(:,:,:) * zfact * tmask(:,:,:) * zysopt(:,:,:) ) ! biogenic silica production
CALL iom_put( "PFeN" , zprofen(:,:,:) * zfact * tmask(:,:,:) ) ! biogenic iron production by nanophyto
CALL iom_put( "PFeD" , zprofed(:,:,:) * zfact * tmask(:,:,:) ) ! biogenic iron production by diatomes
ALLOCATE( zw3d(A2D(0),jpk) ) ; zw3d(A2D(0),jpk) = 0._wp
!
IF( iom_use ( "PPPHYN" ) ) THEN ! primary production by nanophyto
zw3d(A2D(0),1:jpkm1) = zprorcan(A2D(0),1:jpkm1) * zfact * tmask(A2D(0),1:jpkm1)
CALL iom_put( "PPPHYN", zw3d )
ENDIF
!
IF( iom_use ( "PPPHYD" ) ) THEN ! primary production by diatomes
zw3d(A2D(0),1:jpkm1) = zprorcad(A2D(0),1:jpkm1) * zfact * tmask(A2D(0),1:jpkm1)
CALL iom_put( "PPPHYD", zw3d )
ENDIF
!
IF( iom_use ( "PPNEWN" ) ) THEN ! new primary production by nano
zw3d(A2D(0),1:jpkm1) = zpronewn(A2D(0),1:jpkm1) * zfact * tmask(A2D(0),1:jpkm1)
CALL iom_put( "PPNEWN", zw3d )
ENDIF
!
IF( iom_use ( "PPNEWD" ) ) THEN ! new primary production by diatomes
zw3d(A2D(0),1:jpkm1) = zpronewd(A2D(0),1:jpkm1) * zfact * tmask(A2D(0),1:jpkm1)
CALL iom_put( "PPNEWD", zw3d )
ENDIF
!
IF( iom_use ( "PBSi" ) ) THEN ! biogenic silica production
zw3d(A2D(0),1:jpkm1) = zprorcad(A2D(0),1:jpkm1) * zysopt(A2D(0),1:jpkm1) * zfact * tmask(A2D(0),1:jpkm1)
CALL iom_put( "PBSi", zw3d )
ENDIF
!
IF( iom_use ( "PFeN" ) ) THEN ! biogenic iron production by nanophyto
zw3d(A2D(0),1:jpkm1) = zprofen(A2D(0),1:jpkm1) * zfact * tmask(A2D(0),1:jpkm1)
CALL iom_put( "PFeN", zw3d )
ENDIF
!
IF( iom_use ( "PFeD" ) ) THEN ! biogenic iron production by nanophyto
zw3d(A2D(0),1:jpkm1) = zprofed(A2D(0),1:jpkm1) * zfact * tmask(A2D(0),1:jpkm1)
CALL iom_put( "PFeD", zw3d )
ENDIF
!
IF( ln_ligand .AND. ( iom_use( "LPRODP" ) .OR. iom_use( "LDETP" ) ) ) THEN
ALLOCATE( zpligprod(jpi,jpj,jpk) )
zpligprod(:,:,:) = excretd * zprorcad(:,:,:) + excretn * zprorcan(:,:,:)
CALL iom_put( "LPRODP" , zpligprod(:,:,:) * ldocp * 1e9 * zfact * tmask(:,:,:) )
zw3d(A2D(0),1:jpkm1) = ( excretd * zprorcad(A2D(0),1:jpkm1) + excretn * zprorcan(A2D(0),1:jpkm1) ) &
& * zfact * tmask(A2D(0),1:jpkm1)
CALL iom_put( "LPRODP" , zw3d * ldocp * 1e9 )
!
zpligprod(:,:,:) = ( texcretn * zprofen(:,:,:) + texcretd * zprofed(:,:,:) ) &
& * plig(:,:,:) / ( rtrn + plig(:,:,:) + 75.0 * (1.0 - plig(:,:,:) ) )
CALL iom_put( "LDETP" , zpligprod(:,:,:) * lthet * 1e9 * zfact * tmask(:,:,:) )
DEALLOCATE( zpligprod )
zw3d(A2D(0),1:jpkm1) = ( texcretn * zprofen(A2D(0),1:jpkm1) + texcretd * zprofed(A2D(0),1:jpkm1) ) &
& * plig(A2D(0),1:jpkm1) / ( rtrn + plig(A2D(0),1:jpkm1) + 75.0 * (1.0 - plig(A2D(0),1:jpkm1) ) ) &
& * zfact * tmask(A2D(0),1:jpkm1)
CALL iom_put( "LDETP" , zw3d * lthet * 1e9 )
ENDIF
!
IF( iom_use ( "Mumax" ) ) THEN ! Maximum growth rate
zw3d(A2D(0),1:jpkm1) = zprmaxn(A2D(0),1:jpkm1) * tmask(A2D(0),1:jpkm1)
CALL iom_put( "Mumax", zw3d )
ENDIF
CALL iom_put( "Mumax" , zprmaxn(:,:,:) * tmask(:,:,:) ) ! Maximum growth rate
CALL iom_put( "MuN" , zprbio(:,:,:) * xlimphy(:,:,:) * tmask(:,:,:) ) ! Realized growth rate for nanophyto
CALL iom_put( "MuD" , zprdia(:,:,:) * xlimdia(:,:,:) * tmask(:,:,:) ) ! Realized growth rate for diatoms
CALL iom_put( "LNlight" , zprbio (:,:,:) / (zprmaxn(:,:,:) + rtrn) * tmask(:,:,:) ) ! light limitation term
CALL iom_put( "LDlight" , zprdia (:,:,:) / (zprmaxd(:,:,:) + rtrn) * tmask(:,:,:) )
CALL iom_put( "TPP" , ( zprorcan(:,:,:) + zprorcad(:,:,:) ) * zfact * tmask(:,:,:) ) ! total primary production
CALL iom_put( "TPNEW" , ( zpronewn(:,:,:) + zpronewd(:,:,:) ) * zfact * tmask(:,:,:) ) ! total new production
CALL iom_put( "TPBFE" , ( zprofen(:,:,:) + zprofed(:,:,:) ) * zfact * tmask(:,:,:) ) ! total biogenic iron production
CALL iom_put( "tintpp" , tpp * zfact ) ! global total integrated primary production molC/s
!
IF( iom_use ( "MuN" ) ) THEN ! Realized growth rate for nanophyto
zw3d(A2D(0),1:jpkm1) = zprbio(A2D(0),1:jpkm1) * xlimphy(A2D(0),1:jpkm1) * tmask(A2D(0),1:jpkm1)
CALL iom_put( "MuN", zw3d )
ENDIF
!
IF( iom_use ( "MuD" ) ) THEN ! Realized growth rate for diatoms
zw3d(A2D(0),1:jpkm1) = zprdia(A2D(0),1:jpkm1) * xlimdia(A2D(0),1:jpkm1) * tmask(A2D(0),1:jpkm1)
CALL iom_put( "MuD", zw3d )
ENDIF
!
IF( iom_use ( "LNlight" ) ) THEN ! light limitation term for nano
zw3d(A2D(0),1:jpkm1) = zprbio(A2D(0),1:jpkm1) / ( zprmaxn(A2D(0),1:jpkm1) + rtrn ) * tmask(A2D(0),1:jpkm1)
CALL iom_put( "LNlight", zw3d )
ENDIF
!
IF( iom_use ( "LDlight" ) ) THEN ! light limitation term for diatomes
zw3d(A2D(0),1:jpkm1) = zprdia(A2D(0),1:jpkm1) / ( zprmaxd(A2D(0),1:jpkm1) + rtrn ) * tmask(A2D(0),1:jpkm1)
CALL iom_put( "LDlight", zw3d )
ENDIF
!
IF( iom_use ( "TPP" ) ) THEN ! total primary production
zw3d(A2D(0),1:jpkm1) = ( zprorcan(A2D(0),1:jpkm1) + zprorcad(A2D(0),1:jpkm1) ) * zfact * tmask(A2D(0),1:jpkm1)
CALL iom_put( "TPP", zw3d )
ENDIF
!
IF( iom_use ( "TPNEW" ) ) THEN ! total new production
zw3d(A2D(0),1:jpkm1) = ( zpronewn(A2D(0),1:jpkm1) + zpronewd(A2D(0),1:jpkm1) ) * zfact * tmask(A2D(0),1:jpkm1)
CALL iom_put( "TPNEW", zw3d )
ENDIF
!
IF( iom_use ( "TPBFE" ) ) THEN ! total biogenic iron production
zw3d(A2D(0),1:jpkm1) = ( zprofen(A2D(0),1:jpkm1) + zprofed(A2D(0),1:jpkm1) ) * zfact * tmask(A2D(0),1:jpkm1)
CALL iom_put( "TPBFE", zw3d )
ENDIF
!
IF( iom_use( "tintpp") ) CALL iom_put( "tintpp" , tpp * zfact ) ! global total integrated primary production molC/s
!
DEALLOCATE( zw3d )
ENDIF
IF(sn_cfctl%l_prttrc) THEN ! print mean trends (used for debugging)
......@@ -472,6 +542,9 @@ CONTAINS
texcretn = 1._wp - excretn
texcretd = 1._wp - excretd
tpp = 0._wp
! CE not really needed ; tempory, shoub be removed when quotan( A2D(0),jpk )
! quotan(:,:,:) = 0._wp
! quotad(:,:,:) = 0._wp
!
END SUBROUTINE p4z_prod_init
......@@ -480,7 +553,7 @@ CONTAINS
!!----------------------------------------------------------------------
!! *** ROUTINE p4z_prod_alloc ***
!!----------------------------------------------------------------------
ALLOCATE( quotan(jpi,jpj,jpk), quotad(jpi,jpj,jpk), STAT = p4z_prod_alloc )
ALLOCATE( quotan(A2D(0),jpk), quotad(A2D(0),jpk), STAT = p4z_prod_alloc )
!
IF( p4z_prod_alloc /= 0 ) CALL ctl_stop( 'STOP', 'p4z_prod_alloc : failed to allocate arrays.' )
!
......
......@@ -73,11 +73,12 @@ CONTAINS
INTEGER :: ji, jj, jk
REAL(wp) :: zremik, zremikc, zremikn, zremikp, zsiremin, zfact
REAL(wp) :: zsatur, zsatur2, znusil, znusil2, zdep, zdepmin, zfactdep
REAL(wp) :: zbactfer, zonitr, zrfact2
REAL(wp) :: zbactfer, zonitr
REAL(wp) :: zammonic, zoxyremc, zosil, ztem, zdenitnh4, zolimic
CHARACTER (len=25) :: charout
REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdepbac, zolimi, zfacsi, zfacsib, zdepeff, zfebact
REAL(wp), DIMENSION(jpi,jpj ) :: ztempbac
REAL(wp), DIMENSION(A2D(0),jpk) :: zdepbac, zolimi, zfacsi, zfacsib, zdepeff, zfebact
REAL(wp), DIMENSION(A2D(0) ) :: ztempbac
REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: zw3d
!!---------------------------------------------------------------------
!
IF( ln_timing ) CALL timing_start('p4z_rem')
......@@ -90,11 +91,11 @@ CONTAINS
! Computation of the mean bacterial concentration
! this parameterization has been deduced from a model version
! that was modeling explicitely bacteria. This is a very old param
! that was modeling explicitely bacteria. This is a very old parame
! that will be very soon updated based on results from a much more
! recent version of PISCES with bacteria.
! ----------------------------------------------------------------
DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1)
DO_3D( 0, 0, 0, 0, 1, jpkm1)
zdep = MAX( hmld(ji,jj), heup_01(ji,jj) )
IF ( gdept(ji,jj,jk,Kmm) < zdep ) THEN
zdepbac(ji,jj,jk) = 0.6 * ( MAX(0.0, tr(ji,jj,jk,jpzoo,Kbb) + tr(ji,jj,jk,jpmes,Kbb) ) * 1.0E6 )**0.6 * 1.E-6
......@@ -108,7 +109,7 @@ CONTAINS
ENDIF
END_3D
DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1)
DO_3D( 0, 0, 0, 0, 1, jpkm1)
! DOC ammonification. Depends on depth, phytoplankton biomass
! and a limitation term which is supposed to be a parameterization of the bacterial activity.
! --------------------------------------------------------------------------
......@@ -152,7 +153,7 @@ CONTAINS
ENDIF
END_3D
DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1)
DO_3D( 0, 0, 0, 0, 1, jpkm1)
! NH4 nitrification to NO3. Ceased for oxygen concentrations
! below 2 umol/L. Inhibited at strong light
! ----------------------------------------------------------
......@@ -174,7 +175,7 @@ CONTAINS
CALL prt_ctl(tab4d_1=tr(:,:,:,:,Krhs), mask1=tmask, clinfo=ctrcnm)
ENDIF
DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1)
DO_3D( 0, 0, 0, 0, 1, jpkm1)
! Bacterial uptake of iron. No iron is available in DOC. So
! Bacteries are obliged to take up iron from the water. Some
......@@ -202,7 +203,7 @@ CONTAINS
! Initialization of the array which contains the labile fraction
! of bSi. Set to a constant in the upper ocean
! ---------------------------------------------------------------
DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1)
DO_3D( 0, 0, 0, 0, 1, jpkm1)
! Remineralization rate of BSi dependent on T and saturation
! The parameterization is taken from Ridgwell et al. (2002)
! ---------------------------------------------------------
......@@ -236,20 +237,34 @@ CONTAINS
WRITE(charout, FMT="('rem3')")
CALL prt_ctl_info( charout, cdcomp = 'top' )
CALL prt_ctl(tab4d_1=tr(:,:,:,:,Krhs), mask1=tmask, clinfo=ctrcnm)
ENDIF
ENDIF
IF( knt == nrdttrc ) THEN
zrfact2 = 1.e+3 * rfact2r ! conversion from mol/l/kt to mol/m3/s
!
IF( iom_use( "REMIN" ) ) THEN ! Remineralisation rate
zolimi(:,:,jpk) = 0. ; CALL iom_put( "REMIN" , zolimi(:,:,:) * tmask(:,:,:) * zrfact2 )
ALLOCATE( zw3d(A2D(0),jpk) ) ; zw3d(A2D(0),jpk) = 0._wp
!
IF( iom_use ( "REMIN" ) ) THEN ! Remineralisation rate
zw3d(A2D(0),1:jpkm1) = zolimi(A2D(0),1:jpkm1) * rfact2r * tmask(A2D(0),1:jpkm1)
CALL iom_put( "REMIN", zw3d )
ENDIF
CALL iom_put( "DENIT" , denitr(:,:,:) * rdenit * rno3 * tmask(:,:,:) * zrfact2 ) ! Denitrification
IF( iom_use( "BACT" ) ) THEN ! Bacterial biomass
zdepbac(:,:,jpk) = 0. ; CALL iom_put( "BACT", zdepbac(:,:,:) * 1.E6 * tmask(:,:,:) )
!
IF( iom_use ( "BACT" ) ) THEN ! Bacterial biomass
zw3d(A2D(0),1:jpkm1) = zdepbac(A2D(0),1:jpkm1) * 1.E6 * tmask(A2D(0),1:jpkm1)
CALL iom_put( "BACT", zw3d )
ENDIF
CALL iom_put( "FEBACT" , zfebact(:,:,:) * 1E9 * tmask(:,:,:) * zrfact2 )
ENDIF
!
IF( iom_use ( "FEBACT" ) ) THEN
zw3d(A2D(0),1:jpkm1) = zfebact(A2D(0),1:jpkm1) * 1E9 * rfact2r * tmask(A2D(0),1:jpkm1)
CALL iom_put( "FEBACT", zw3d )
ENDIF
!
IF( iom_use( "DENIT" ) ) THEN ! Denitrification
zw3d(A2D(0),1:jpkm1) = denitr(A2D(0),1:jpkm1) * 1E+3 * rfact2r * rno3 * tmask(A2D(0),1:jpkm1)
CALL iom_put( "DENIT", zw3d )
ENDIF
!
DEALLOCATE( zw3d )
ENDIF
!
IF( ln_timing ) CALL timing_stop('p4z_rem')
!
......