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 231 additions and 171 deletions
......@@ -217,7 +217,6 @@ CONTAINS
IF( lk_oasis ) THEN ; cxios_context = 'sas' ! when coupling SAS to OCE
ELSE ; cxios_context = 'nemo' !
ENDIF
nn_hls = 1
!
l_SAS = .TRUE. ! used in domain:dom_nam
!
......@@ -390,6 +389,7 @@ CONTAINS
#if defined key_agrif
uu(:,:,:,:) = 0.0_wp ; vv(:,:,:,:) = 0.0_wp ; ts(:,:,:,:,:) = 0.0_wp ! needed for interp done at initialization phase
uu_b(:,:,:) = 0.0_wp ; vv_b(:,:,:) = 0.0_wp
#endif
! ! external forcing
CALL sbc_init( Nbb, Nnn, Naa ) ! Forcings : surface module
......
......@@ -36,6 +36,9 @@ MODULE stpctl
INTEGER, PARAMETER :: jpvar = 3
INTEGER :: nrunid ! netcdf file id
INTEGER, DIMENSION(jpvar) :: nvarid ! netcdf variable id
!! * Substitutions
# include "do_loop_substitute.h90"
!!----------------------------------------------------------------------
!! NEMO/SAS 4.0 , NEMO Consortium (2018)
!! $Id: stpctl.F90 14433 2021-02-11 08:06:49Z smasson $
......@@ -77,8 +80,8 @@ CONTAINS
IF( nstop > 0 .AND. ngrdstop > -1 ) RETURN ! stpctl was already called by a child grid
!
ll_wrtstp = ( MOD( kt-nit000, sn_cfctl%ptimincr ) == 0 ) .OR. ( kt == nitend )
ll_colruns = ll_wrtstp .AND. sn_cfctl%l_runstat .AND. jpnij > 1
ll_wrtruns = ( ll_colruns .OR. jpnij == 1 ) .AND. lwm
ll_colruns = sn_cfctl%l_runstat .AND. ll_wrtstp .AND. jpnij > 1
ll_wrtruns = sn_cfctl%l_runstat .AND. ll_wrtstp .AND. lwm
!
IF( kt == nit000 ) THEN
!
......@@ -132,7 +135,7 @@ CONTAINS
!
zmax(1) = MAXVAL( vt_i (:,:) , mask = llmsk ) ! max ice thickness
zmax(2) = MAXVAL( ABS( u_ice(:,:) ) , mask = llmsk ) ! max ice velocity (zonal only)
zmax(3) = MAXVAL( -tm_i (:,:) + rt0, mask = llmsk ) ! min ice temperature (in degC)
zmax(3) = MAXVAL( -tm_i (:,:) + rt0, mask = llmsk(A2D(0)) ) ! min ice temperature (in degC)
zmax(jpvar+1) = REAL( nstop, wp ) ! stop indicator
!
! !== get global extrema ==!
......@@ -172,9 +175,9 @@ CONTAINS
! first: close the netcdf file, so we can read it
IF( lwm .AND. kt /= nitend ) istatus = NF90_CLOSE(nrunid)
! get global loc on the min/max
CALL mpp_maxloc( 'stpctl', vt_i(:,:) , llmsk, zzz, iloc(1:2,1) ) ! mpp_maxloc ok if mask = F
CALL mpp_maxloc( 'stpctl',ABS( u_ice(:,:) ) , llmsk, zzz, iloc(1:2,2) )
CALL mpp_minloc( 'stpctl', tm_i(:,:) - rt0, llmsk, zzz, iloc(1:2,3) )
CALL mpp_maxloc( 'stpctl', vt_i(:,:) , llmsk , zzz, iloc(1:2,1) ) ! mpp_maxloc ok if mask = F
CALL mpp_maxloc( 'stpctl',ABS( u_ice(:,:) ) , llmsk , zzz, iloc(1:2,2) )
CALL mpp_minloc( 'stpctl', tm_i(:,:) - rt0, llmsk(A2D(0)), zzz, iloc(1:2,3) )
! find which subdomain has the max.
iareamin(:) = jpnij+1 ; iareamax(:) = 0 ; iareasum(:) = 0
DO ji = 1, jptst
......@@ -187,11 +190,11 @@ CONTAINS
CALL mpp_sum( "stpctl", iareasum ) ! sum over the global domain
ELSE ! find local min and max locations:
! if we are here, this means that the subdomain contains some oce points -> no need to test the mask used in maxloc
iloc(1:2,1) = MAXLOC( vt_i(:,:) , mask = llmsk )
iloc(1:2,2) = MAXLOC( ABS( u_ice(:,:) ) , mask = llmsk )
iloc(1:2,3) = MINLOC( tm_i(:,:) - rt0, mask = llmsk )
iloc(1:2,1) = MAXLOC( vt_i(:,:) , mask = llmsk )
iloc(1:2,2) = MAXLOC( ABS( u_ice(:,:) ) , mask = llmsk )
iloc(1:2,3) = MINLOC( tm_i(:,:) - rt0, mask = llmsk(A2D(0)) )
DO ji = 1, jptst ! local domain indices ==> global domain indices, excluding halos
iloc(1:2,ji) = (/ mig0(iloc(1,ji)), mjg0(iloc(2,ji)) /)
iloc(1:2,ji) = (/ mig(iloc(1,ji),0), mjg(iloc(2,ji),0) /)
END DO
iareamin(:) = narea ; iareamax(:) = narea ; iareasum(:) = 1 ! this is local information
ENDIF
......
......@@ -75,8 +75,8 @@ CONTAINS
IF( nstop > 0 .AND. ngrdstop > -1 ) RETURN ! stpctl was already called by a child grid
!
ll_wrtstp = ( MOD( kt-nit000, sn_cfctl%ptimincr ) == 0 ) .OR. ( kt == nitend )
ll_colruns = ll_wrtstp .AND. sn_cfctl%l_runstat .AND. jpnij > 1
ll_wrtruns = ( ll_colruns .OR. jpnij == 1 ) .AND. lwm
ll_colruns = sn_cfctl%l_runstat .AND. ll_wrtstp .AND. jpnij > 1
ll_wrtruns = sn_cfctl%l_runstat .AND. ll_wrtstp .AND. lwm
!
IF( kt == nit000 ) THEN
!
......@@ -185,7 +185,7 @@ CONTAINS
llmsk(Nis0:Nie0,Njs0:Nje0,:) = umask(Nis0:Nie0,Njs0:Nje0,:) == 1._wp ! define only the inner domain
iloc(1:3,2) = MAXLOC( ABS( uu(:,:,:, Kmm)), mask = llmsk(:,:,:) )
DO ji = 1, jptst ! local domain indices ==> global domain indices, excluding halos
iloc(1:2,ji) = (/ mig0(iloc(1,ji)), mjg0(iloc(2,ji)) /)
iloc(1:2,ji) = (/ mig(iloc(1,ji),0), mjg(iloc(2,ji),0) /)
END DO
iareamin(:) = narea ; iareamax(:) = narea ; iareasum(:) = 1 ! this is local information
ENDIF
......
......@@ -114,9 +114,9 @@ CONTAINS
zrhs_u = - grav * ( ssh(ji+1,jj,Nnn) - ssh(ji,jj,Nnn) ) * r1_e1u(ji,jj)
zrhs_v = - grav * ( ssh(ji,jj+1,Nnn) - ssh(ji,jj,Nnn) ) * r1_e2v(ji,jj)
! ! wind stress and layer friction
zrhs_u = zrhs_u + z1_2rho0 * ( utau_b(ji,jj) + utau(ji,jj) ) / e3u(ji,jj,jk,Nnn) &
zrhs_u = zrhs_u + z1_2rho0 * ( utau_b(ji,jj) + utauU(ji,jj) ) / e3u(ji,jj,jk,Nnn) &
& - rn_rfr * uu(ji,jj,jk,Nbb)
zrhs_v = zrhs_v + z1_2rho0 * ( vtau_b(ji,jj) + vtau(ji,jj) ) / e3v(ji,jj,jk,Nnn) &
zrhs_v = zrhs_v + z1_2rho0 * ( vtau_b(ji,jj) + vtauV(ji,jj) ) / e3v(ji,jj,jk,Nnn) &
& - rn_rfr * vv(ji,jj,jk,Nbb)
! ! ==> RHS
uu(ji,jj,jk,Nrhs) = uu(ji,jj,jk,Nrhs) + zrhs_u
......
......@@ -135,9 +135,9 @@ CONTAINS
zrhs_v = - grav * ( ssh(ji,jj+1,Nbb) - ssh(ji,jj,Nbb) ) * r1_e2v(ji,jj)
#if defined key_RK3all
! ! wind stress and layer friction
zrhs_u = zrhs_u + r1_rho0 * ( z5_6*utau_b(ji,jj) + (1._wp - z5_6)*utau(ji,jj) ) / e3u(ji,jj,jk,Nbb) &
zrhs_u = zrhs_u + r1_rho0 * ( z5_6*utau_b(ji,jj) + (1._wp - z5_6)*utauU(ji,jj) ) / e3u(ji,jj,jk,Nbb) &
& - rn_rfr * uu(ji,jj,jk,Nbb)
zrhs_v = zrhs_v + r1_rho0 * ( z5_6*vtau_b(ji,jj) + (1._wp - z5_6)*vtau(ji,jj) ) / e3v(ji,jj,jk,Nbb) &
zrhs_v = zrhs_v + r1_rho0 * ( z5_6*vtau_b(ji,jj) + (1._wp - z5_6)*vtauV(ji,jj) ) / e3v(ji,jj,jk,Nbb) &
& - rn_rfr * vv(ji,jj,jk,Nbb)
#endif
! ! ==> RHS
......@@ -201,9 +201,9 @@ CONTAINS
zrhs_v = - grav * ( ssh(ji,jj+1,Nnn) - ssh(ji,jj,Nnn) ) * r1_e2v(ji,jj)
#if defined key_RK3all
! ! wind stress and layer friction
zrhs_u = zrhs_u + r1_rho0 * ( z3_4*utau_b(ji,jj) + (1._wp - z3_4)*utau(ji,jj) ) / e3u(ji,jj,jk,Nnn) &
zrhs_u = zrhs_u + r1_rho0 * ( z3_4*utau_b(ji,jj) + (1._wp - z3_4)*utauU(ji,jj) ) / e3u(ji,jj,jk,Nnn) &
& - rn_rfr * uu(ji,jj,jk,Nbb)
zrhs_v = zrhs_v + r1_rho0 * ( z3_4*vtau_b(ji,jj) + (1._wp - z3_4)*vtau(ji,jj) ) / e3v(ji,jj,jk,Nnn) &
zrhs_v = zrhs_v + r1_rho0 * ( z3_4*vtau_b(ji,jj) + (1._wp - z3_4)*vtauV(ji,jj) ) / e3v(ji,jj,jk,Nnn) &
& - rn_rfr * vv(ji,jj,jk,Nbb)
#endif
! ! ==> RHS
......@@ -265,9 +265,9 @@ CONTAINS
zrhs_u = - grav * ( ssh(ji+1,jj,Nnn) - ssh(ji,jj,Nnn) ) * r1_e1u(ji,jj)
zrhs_v = - grav * ( ssh(ji,jj+1,Nnn) - ssh(ji,jj,Nnn) ) * r1_e2v(ji,jj)
! ! wind stress and layer friction
zrhs_u = zrhs_u + z1_2rho0 * ( utau_b(ji,jj) + utau(ji,jj) ) / e3u(ji,jj,jk,Nnn) &
zrhs_u = zrhs_u + z1_2rho0 * ( utau_b(ji,jj) + utauU(ji,jj) ) / e3u(ji,jj,jk,Nnn) &
& - rn_rfr * uu(ji,jj,jk,Nbb)
zrhs_v = zrhs_v + z1_2rho0 * ( vtau_b(ji,jj) + vtau(ji,jj) ) / e3v(ji,jj,jk,Nnn) &
zrhs_v = zrhs_v + z1_2rho0 * ( vtau_b(ji,jj) + vtauV(ji,jj) ) / e3v(ji,jj,jk,Nnn) &
& - rn_rfr * vv(ji,jj,jk,Nbb)
! ! ==> RHS
uu(ji,jj,jk,Nrhs) = uu(ji,jj,jk,Nrhs) + zrhs_u
......
......@@ -46,7 +46,7 @@ CONTAINS
!!----------------------------------------------------------------------
INTEGER, INTENT(in) :: kt ! ocean time-step index
INTEGER, INTENT(in) :: Kbb, Kmm, Krhs ! ocean time level
INTEGER :: jn, jk ! dummy loop index
INTEGER :: jk ! dummy loop index
!!----------------------------------------------------------------------
!
IF( ln_timing ) CALL timing_start('trc_sms_age')
......@@ -74,7 +74,7 @@ CONTAINS
tr(:,:,jk,jp_age,Krhs) = tmask(:,:,jk) * rryear
END DO
!
IF( l_trdtrc ) CALL trd_trc( tr(:,:,:,jp_age,Krhs), jn, jptra_sms, kt, Kmm ) ! save trends
IF( l_trdtrc ) CALL trd_trc( tr(:,:,:,jp_age,Krhs), jp_age, jptra_sms, kt, Kmm ) ! save trends
!
IF( ln_timing ) CALL timing_stop('trc_sms_age')
!
......
......@@ -28,7 +28,6 @@ CONTAINS
!!---------------------------------------------------------------------
INTEGER, INTENT(in) :: Kmm ! time level indices
CHARACTER (len=20) :: cltra
INTEGER :: jn
!!---------------------------------------------------------------------
! write the tracer concentrations in the file
......
......@@ -51,6 +51,8 @@ MODULE sms_c14
REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:) :: spco2 ! Atmospheric CO2
REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:) :: tyrco2 ! Time (yr) atmospheric CO2 data
!! * Substitutions
# include "do_loop_substitute.h90"
!!----------------------------------------------------------------------
!! NEMO/TOP 4.0 , NEMO Consortium (2018)
!! $Id: sms_c14.F90 10071 2018-08-28 14:49:04Z nicolasmartin $
......@@ -64,9 +66,9 @@ CONTAINS
!! *** ROUTINE trc_sms_c14_alloc ***
!!----------------------------------------------------------------------
sms_c14_alloc = 0
ALLOCATE( exch_c14(jpi,jpj) , exch_co2(jpi,jpj) , &
& qtr_c14(jpi,jpj) , qint_c14(jpi,jpj) , &
& c14sbc(jpi,jpj) , STAT = sms_c14_alloc )
ALLOCATE( exch_c14(A2D(0)) , exch_co2(A2D(0)) , &
& qtr_c14(A2D(0)) , qint_c14(A2D(0)) , &
& c14sbc(A2D(0)) , STAT = sms_c14_alloc )
!
!
END FUNCTION sms_c14_alloc
......
......@@ -59,6 +59,13 @@ CONTAINS
!
tyrc14_now = 0._wp ! initialize
!
IF( kc14typ == 0) THEN
co2sbc=pco2at
DO_2D( 0, 0, 0, 0 )
c14sbc(ji,jj) = rc14at
END_2D
ENDIF
!
IF(kc14typ >= 1) THEN ! Transient atmospheric forcing: CO2
!
clfile = TRIM( cfileco2 )
......@@ -116,10 +123,10 @@ CONTAINS
! Linear interpolation of the C-14 source fonction
! in linear latitude bands (20N,40N) and (20S,40S)
!------------------------------------------------------
ALLOCATE( fareaz (jpi,jpj ,nc14zon) , STAT=ierr3 )
ALLOCATE( fareaz(A2D(0) ,nc14zon) , STAT=ierr3 )
IF( ierr3 /= 0 ) CALL ctl_stop( 'STOP', 'trc_atm_c14_ini: unable to allocate fareaz' )
!
DO_2D( 1, 1, 1, 1 ) ! from C14b package
DO_2D( 0, 0, 0, 0 ) ! from C14b package
IF( gphit(ji,jj) >= yn40 ) THEN
fareaz(ji,jj,1) = 0.
fareaz(ji,jj,2) = 0.
......@@ -205,9 +212,9 @@ CONTAINS
!! ** Action : atmospheric values interpolated at time-step kt
!!----------------------------------------------------------------------
INTEGER , INTENT(in ) :: kt ! ocean time-step
REAL(wp), DIMENSION(:,:), INTENT( out) :: c14sbc ! atm c14 ratio
REAL(wp), DIMENSION(A2D(0)), INTENT( out) :: c14sbc ! atm c14 ratio
REAL(wp), INTENT( out) :: co2sbc ! atm co2 p
INTEGER :: jz ! dummy loop indice
INTEGER :: ji, jj, jz ! dummy loop indice
REAL(wp) :: zdint,zint ! work
REAL(wp), DIMENSION(nc14zon) :: zonbc14 ! work
!
......@@ -215,10 +222,6 @@ CONTAINS
!
IF( ln_timing ) CALL timing_start('trc_atm_c14')
!
IF( kc14typ == 0) THEN
co2sbc=pco2at
c14sbc(:,:)=rc14at
ENDIF
!
IF(kc14typ >= 1) THEN ! Transient C14 & CO2
!
......
......@@ -80,7 +80,7 @@ CONTAINS
! CO2 solubility (Weiss, 1974; Wanninkhof, 2014)
! -------------------------------------------------------------------
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
DO_2D( 0, 0, 0, 0 )
IF( tmask(ji,jj,1) > 0. ) THEN
!
zt = MIN( 40. , ts(ji,jj,1,jp_tem,Kmm) )
......@@ -121,21 +121,21 @@ CONTAINS
!
! Flux of C-14 from air-to-sea; units: (C14/C ratio) x m/s
! already masked
qtr_c14(:,:) = exch_c14(:,:) * ( c14sbc(:,:) - tr(:,:,1,jp_c14,Kbb) )
DO_2D( 0, 0, 0, 0 )
qtr_c14(ji,jj) = exch_c14(ji,jj) * ( c14sbc(ji,jj) - tr(ji,jj,1,jp_c14,Kbb) )
END_2D
! cumulation of air-to-sea flux at each time step
qint_c14(:,:) = qint_c14(:,:) + qtr_c14(:,:) * rn_Dt
!
! Add the surface flux to the trend of jp_c14
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
DO_2D( 0, 0, 0, 0 )
tr(ji,jj,1,jp_c14,Krhs) = tr(ji,jj,1,jp_c14,Krhs) + qtr_c14(ji,jj) / e3t(ji,jj,1,Kmm)
END_2D
!
! Computation of decay effects on jp_c14
DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpk )
!
DO_3D( 0, 0, 0, 0, 1, jpkm1 )
tr(ji,jj,jk,jp_c14,Krhs) = tr(ji,jj,jk,jp_c14,Krhs) - rlam14 * tr(ji,jj,jk,jp_c14,Kbb) * tmask(ji,jj,jk)
!
END_3D
!
IF( lrst_trc ) THEN
......
......@@ -38,8 +38,8 @@ CONTAINS
CHARACTER (len=20) :: cltra ! short title for tracer
INTEGER :: ji,jj,jk,jn ! dummy loop indexes
REAL(wp) :: zage,zarea,ztemp ! temporary
REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zres, z2d ! temporary storage 2D
REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: z3d , zz3d ! temporary storage 3D
REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: z2d ! temporary storage 2D
REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: z3d ! temporary storage 3D
!!---------------------------------------------------------------------
! write the tracer concentrations in the file
......@@ -49,41 +49,35 @@ CONTAINS
! compute and write the tracer diagnostic in the file
! ---------------------------------------
IF( iom_use("qtr_c14") ) CALL iom_put( "qtr_c14" , rsiyea * qtr_c14(:,:) ) ! Radiocarbon surf flux [./m2/yr]
CALL iom_put( "qint_c14", qint_c14(:,:) ) ! cumulative flux [./m2]
IF( iom_use("DeltaC14") .OR. iom_use("C14Age") .OR. iom_use("RAge") ) THEN
!
ALLOCATE( z2d(jpi,jpj), zres(jpi,jpj) )
ALLOCATE( z3d(jpi,jpj,jpk), zz3d(jpi,jpj,jpk) )
ALLOCATE( z2d(A2D(0)), z3d(A2D(0),jpk) )
!
zage = -1._wp / rlam14 / rsiyea ! factor for radioages in year
z3d(:,:,:) = 1._wp
zz3d(:,:,:) = 0._wp
!
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) > 0._wp) THEN
z3d (ji,jj,jk) = tr(ji,jj,jk,jp_c14,Kmm)
zz3d(ji,jj,jk) = LOG( z3d(ji,jj,jk) )
z3d(ji,jj,jk) = tr(ji,jj,jk,jp_c14,Kmm)
ENDIF
END_3D
zres(:,:) = z3d(:,:,1)
CALL iom_put( "C14Age", zage * LOG( z3d(:,:,:) ) ) ! Radiocarbon age [yr]
! Reservoir age [yr]
z2d(:,:) =0._wp
jk = 1
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
ztemp = zres(ji,jj) / c14sbc(ji,jj)
IF( ztemp > 0._wp .AND. tmask(ji,jj,jk) > 0._wp ) z2d(ji,jj) = LOG( ztemp )
z2d(:,:) = 0._wp
DO_2D( 0, 0, 0, 0 )
ztemp = z3d(ji,jj,1) / c14sbc(ji,jj)
IF( ztemp > 0._wp .AND. tmask(ji,jj,1) > 0._wp ) z2d(ji,jj) = LOG( ztemp )
END_2D
CALL iom_put( "RAge" , zage * z2d(:,:) ) ! Reservoir age [yr]
!
z3d(:,:,:) = 1.d03 * ( z3d(:,:,:) - 1._wp )
CALL iom_put( "DeltaC14" , z3d(:,:,:) ) ! Delta C14 [permil]
CALL iom_put( "C14Age" , zage * zz3d(:,:,:) ) ! Radiocarbon age [yr]
CALL iom_put( "qtr_c14", rsiyea * qtr_c14(:,:) ) ! Radiocarbon surf flux [./m2/yr]
CALL iom_put( "qint_c14" , qint_c14 ) ! cumulative flux [./m2]
CALL iom_put( "RAge" , zage * z2d(:,:) ) ! Reservoir age [yr]
!
DEALLOCATE( z2d, zres, z3d, zz3d )
DEALLOCATE( z2d, z3d )
!
ENDIF
!
......@@ -91,23 +85,35 @@ CONTAINS
!
CALL iom_put( "AtmCO2", co2sbc ) ! global atmospheric CO2 [ppm]
IF( iom_use("AtmC14") ) THEN
zarea = glob_sum( 'trcwri_c14', e1e2t(:,:) ) ! global ocean surface
ztemp = glob_sum( 'trcwri_c14', c14sbc(:,:) * e1e2t(:,:) )
ztemp = ( ztemp / zarea - 1._wp ) * 1000._wp
CALL iom_put( "AtmC14" , ztemp ) ! Global atmospheric DeltaC14 [permil]
ENDIF
IF( iom_use("K_C14") ) THEN
ztemp = glob_sum ( 'trcwri_c14', exch_c14(:,:) * e1e2t(:,:) )
ztemp = rsiyea * ztemp / zarea
CALL iom_put( "K_C14" , ztemp ) ! global mean exchange velocity for C14/C ratio [m/yr]
ENDIF
IF( iom_use("K_CO2") ) THEN
IF( iom_use("AtmC14") .OR. iom_use("K_C14") .OR. iom_use("K_CO2") ) THEN
zarea = glob_sum( 'trcwri_c14', e1e2t(:,:) ) ! global ocean surface
ztemp = glob_sum ( 'trcwri_c14', exch_co2(:,:) * e1e2t(:,:) )
ztemp = 360000._wp * ztemp / zarea ! cm/h units: directly comparable with literature
CALL iom_put( "K_CO2", ztemp ) ! global mean CO2 piston velocity [cm/hr]
ENDIF
ALLOCATE( z2d(A2D(0)) )
IF( iom_use("AtmC14") ) THEN
DO_2D( 0, 0, 0, 0 )
z2d(ji,jj) = c14sbc(ji,jj) * e1e2t(ji,jj)
END_2D
ztemp = glob_sum( 'trcwri_c14', z2d(:,:) )
ztemp = ( ztemp / zarea - 1._wp ) * 1000._wp
CALL iom_put( "AtmC14" , ztemp ) ! Global atmospheric DeltaC14 [permil]
ENDIF
IF( iom_use("K_C14") ) THEN
DO_2D( 0, 0, 0, 0 )
z2d(ji,jj) = exch_c14(ji,jj) * e1e2t(ji,jj)
END_2D
ztemp = glob_sum( 'trcwri_c14', z2d(:,:) )
ztemp = rsiyea * ztemp / zarea
CALL iom_put( "K_C14" , ztemp ) ! global mean exchange velocity for C14/C ratio [m/yr]
ENDIF
IF( iom_use("K_CO2") ) THEN
DO_2D( 0, 0, 0, 0 )
z2d(ji,jj) = exch_co2(ji,jj) * e1e2t(ji,jj)
END_2D
ztemp = glob_sum( 'trcwri_c14', z2d(:,:) )
ztemp = 360000._wp * ztemp / zarea ! cm/h units: directly comparable with literature
CALL iom_put( "K_CO2", ztemp ) ! global mean CO2 piston velocity [cm/hr]
ENDIF
DEALLOCATE( z2d )
END IF
IF( iom_use("C14Inv") ) THEN
ztemp = glob_sum( 'trcwri_c14', tr(:,:,:,jp_c14,Kmm) * cvol(:,:,:) )
ztemp = atomc14 * xdicsur * ztemp
......
......@@ -131,7 +131,7 @@ CONTAINS
! Linear interpolation between 2 hemispheric function of latitud between ylats and ylatn
!---------------------------------------------------------------------------------------
zyd = ylatn - ylats
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
DO_2D( 0, 0, 0, 0 )
IF( gphit(ji,jj) >= ylatn ) THEN ; xphem(ji,jj) = 1.e0
ELSEIF( gphit(ji,jj) <= ylats ) THEN ; xphem(ji,jj) = 0.e0
ELSE ; xphem(ji,jj) = ( gphit(ji,jj) - ylats) / zyd
......
......@@ -124,9 +124,9 @@ CONTAINS
& + atm_cfc(iyear_end, jm, jl) * REAL(im2, wp) ) / 12.
END DO
! !------------!
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) ! i-j loop !
! !------------!
! !------------!
DO_2D( 0, 0, 0, 0 ) ! i-j loop !
! !------------!
! space interpolation
zpp_cfc = xphem(ji,jj) * zpatm(1,jl) &
& + ( 1.- xphem(ji,jj) ) * zpatm(2,jl)
......@@ -309,8 +309,8 @@ CONTAINS
!!----------------------------------------------------------------------
!! *** ROUTINE trc_sms_cfc_alloc ***
!!----------------------------------------------------------------------
ALLOCATE( xphem (jpi,jpj) , atm_cfc(jpyear,jphem,jp_cfc) , &
& qtr_cfc (jpi,jpj,jp_cfc) , qint_cfc(jpi,jpj,jp_cfc) , &
ALLOCATE( xphem (A2D(0)) , atm_cfc(jpyear,jphem,jp_cfc) , &
& qtr_cfc (A2D(0),jp_cfc) , qint_cfc(A2D(0),jp_cfc) , &
& soa(4,jp_cfc) , sob(3,jp_cfc) , sca(5,jp_cfc) , &
& STAT=trc_sms_cfc_alloc )
!
......
......@@ -104,7 +104,6 @@ CONTAINS
!
IF( ln_timing ) CALL timing_start('p2z_bio')
!
IF( lk_iomput ) ALLOCATE( zw2d(jpi,jpj,17), zw3d(jpi,jpj,jpk,3) )
IF( kt == nittrc000 ) THEN
IF(lwp) WRITE(numout,*)
......@@ -112,18 +111,18 @@ CONTAINS
IF(lwp) WRITE(numout,*) ' ~~~~~~~'
ENDIF
xksi(:,:) = 0.e0 ! zooplakton closure ( fbod)
IF( lk_iomput ) THEN
zw2d (:,:,:) = 0._wp
zw3d(:,:,:,:) = 0._wp
ALLOCATE( zw3d(A2D(0),jpk,3) ) ; zw3d(:,:,jpk,:) = 0._wp
ALLOCATE( zw2d(A2D(0),17) ) ; zw2d(:,:,:) = 0._wp
ENDIF
!
xksi(:,:) = 0.e0 ! zooplakton closure ( fbod)
! ! -------------------------- !
DO jk = 1, jpkbm1 ! Upper ocean (bio-layers) !
DO_3D( 0, 0, 0, 0, 1, jpkbm1 ) ! Upper ocean (bio-layers) !
! ! -------------------------- !
DO_2D( 0, 0, 0, 0 )
! trophic variables( det, zoo, phy, no3, nh4, dom)
! ------------------------------------------------
! trophic variables( det, zoo, phy, no3, nh4, dom)
! ------------------------------------------------
! negative trophic variables DO not contribute to the fluxes
zdet = MAX( 0.e0, tr(ji,jj,jk,jpdet,Kmm) )
......@@ -235,13 +234,11 @@ CONTAINS
zw3d(ji,jj,jk,3) = znh4no3 * 86400
!
ENDIF
END_2D
END DO
END_3D
! ! -------------------------- !
DO jk = jpkb, jpkm1 ! Upper ocean (bio-layers) !
DO_3D( 0, 0, 0, 0, jpkb, jpkm1 ) ! Upper ocean (bio-layers) !
! ! -------------------------- !
DO_2D( 0, 0, 0, 0 )
! remineralisation of all quantities towards nitrate
! trophic variables( det, zoo, phy, no3, nh4, dom)
......@@ -334,12 +331,9 @@ CONTAINS
zw3d(ji,jj,jk,3) = znh4no3 * 86400._wp
!
ENDIF
END_2D
END DO
END_3D
!
IF( lk_iomput ) THEN
CALL lbc_lnk( 'p2zbio', zw2d(:,:,:),'T', 1.0_wp )
CALL lbc_lnk( 'p2zbio', zw3d(:,:,:,1),'T', 1.0_wp, zw3d(:,:,:,2),'T', 1.0_wp, zw3d(:,:,:,3),'T', 1.0_wp )
! Save diagnostics
CALL iom_put( "TNO3PHY", zw2d(:,:,1) )
CALL iom_put( "TNH4PHY", zw2d(:,:,2) )
......@@ -362,6 +356,8 @@ CONTAINS
CALL iom_put( "FNH4PHY", zw3d(:,:,:,2) )
CALL iom_put( "FNH4NO3", zw3d(:,:,:,3) )
!
DEALLOCATE( zw2d, zw3d )
!
ENDIF
IF(sn_cfctl%l_prttrc) THEN ! print mean trends (used for debugging)
......@@ -370,8 +366,6 @@ CONTAINS
CALL prt_ctl(tab4d_1=tr(:,:,:,:,Krhs), mask1=tmask, clinfo=ctrcnm)
ENDIF
!
IF( lk_iomput ) DEALLOCATE( zw2d, zw3d )
!
IF( ln_timing ) CALL timing_stop('p2z_bio')
!
END SUBROUTINE p2z_bio
......
......@@ -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
......