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
No results found
Show changes
Hide whitespace changes
Inline Side-by-side
Showing
with 246 additions and 434 deletions
src/OCE/DYN/dynvor.F90
View file @ 672ff863
......@@ -115,9 +115,9 @@ CONTAINS
!!----------------------------------------------------------------------
INTEGER , INTENT( in ) :: kt ! ocean time-step index
INTEGER , INTENT( in ) :: Kmm, Krhs ! ocean time level indices
REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) :: puu, pvv ! ocean velocity field and RHS of momentum equation
REAL(dp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) :: puu, pvv ! ocean velocity field and RHS of momentum equation
!
REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: ztrdu, ztrdv
REAL(dp), ALLOCATABLE, DIMENSION(:,:,:) :: ztrdu, ztrdv
!!----------------------------------------------------------------------
!
IF( ln_timing ) CALL timing_start('dyn_vor')
......@@ -234,8 +234,8 @@ CONTAINS
INTEGER , INTENT(in ) :: kt ! ocean time-step index
INTEGER , INTENT(in ) :: Kmm ! ocean time level index
INTEGER , INTENT(in ) :: kvor ! total, planetary, relative, or metric
REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pu, pv ! now velocities
REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pu_rhs, pv_rhs ! total v-trend
REAL(dp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pu, pv ! now velocities
REAL(dp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pu_rhs, pv_rhs ! total v-trend
!
INTEGER :: ji, jj, jk ! dummy loop indices
REAL(wp) :: zx1, zy1, zx2, zy2 ! local scalars
......@@ -354,8 +354,8 @@ CONTAINS
INTEGER , INTENT(in ) :: kt ! ocean time-step index
INTEGER , INTENT(in ) :: Kmm ! ocean time level index
INTEGER , INTENT(in ) :: kvor ! total, planetary, relative, or metric
REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pu, pv ! now velocities
REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pu_rhs, pv_rhs ! total v-trend
REAL(dp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pu, pv ! now velocities
REAL(dp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pu_rhs, pv_rhs ! total v-trend
!
INTEGER :: ji, jj, jk ! dummy loop indices
REAL(wp) :: zx1, zy1, zx2, zy2, ze3f, zmsk ! local scalars
......@@ -487,8 +487,8 @@ CONTAINS
INTEGER , INTENT(in ) :: kt ! ocean time-step index
INTEGER , INTENT(in ) :: Kmm ! ocean time level index
INTEGER , INTENT(in ) :: kvor ! total, planetary, relative, or metric
REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pu, pv ! now velocities
REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pu_rhs, pv_rhs ! total v-trend
REAL(dp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pu, pv ! now velocities
REAL(dp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pu_rhs, pv_rhs ! total v-trend
!
INTEGER :: ji, jj, jk ! dummy loop indices
REAL(wp) :: zuav, zvau, ze3f, zmsk ! local scalars
......@@ -617,8 +617,8 @@ CONTAINS
INTEGER , INTENT(in ) :: kt ! ocean time-step index
INTEGER , INTENT(in ) :: Kmm ! ocean time level index
INTEGER , INTENT(in ) :: kvor ! total, planetary, relative, or metric
REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pu, pv ! now velocities
REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pu_rhs, pv_rhs ! total v-trend
REAL(dp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pu, pv ! now velocities
REAL(dp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pu_rhs, pv_rhs ! total v-trend
!
INTEGER :: ji, jj, jk ! dummy loop indices
INTEGER :: ierr ! local integer
......@@ -816,8 +816,8 @@ CONTAINS
INTEGER , INTENT(in ) :: kt ! ocean time-step index
INTEGER , INTENT(in ) :: Kmm ! ocean time level index
INTEGER , INTENT(in ) :: kvor ! total, planetary, relative, or metric
REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pu, pv ! now velocities
REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pu_rhs, pv_rhs ! total v-trend
REAL(dp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pu, pv ! now velocities
REAL(dp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pu_rhs, pv_rhs ! total v-trend
!
INTEGER :: ji, jj, jk ! dummy loop indices
INTEGER :: ierr ! local integer
......@@ -1001,7 +1001,7 @@ CONTAINS
IF(lwp) WRITE(numout,*) ' ==>>> vector form dynamics : total vorticity = Coriolis + relative vorticity'
nrvm = np_RVO ! relative vorticity
ntot = np_CRV ! relative + planetary vorticity
CASE( np_FLX_c2 , np_FLX_ubs )
CASE( np_FLX_c2 , np_FLX_ubs , np_FLX_up3 )
IF(lwp) WRITE(numout,*) ' ==>>> flux form dynamics : total vorticity = Coriolis + metric term'
nrvm = np_MET ! metric term
ntot = np_CME ! Coriolis + metric term
......
src/OCE/DYN/dynzad.F90
View file @ 672ff863
......@@ -55,13 +55,13 @@ CONTAINS
!!----------------------------------------------------------------------
INTEGER , INTENT( in ) :: kt ! ocean time-step inedx
INTEGER , INTENT( in ) :: Kmm, Krhs ! ocean time level indices
REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) :: puu, pvv ! ocean velocities and RHS of momentum equation
REAL(dp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) :: puu, pvv ! ocean velocities and RHS of momentum equation
!
INTEGER :: ji, jj, jk ! dummy loop indices
REAL(wp) :: zua, zva ! local scalars
REAL(wp), DIMENSION(A2D(nn_hls)) :: zww
REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: zwuw, zwvw
REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: ztrdu, ztrdv
REAL(dp), DIMENSION(:,:,:), ALLOCATABLE :: ztrdu, ztrdv
!!----------------------------------------------------------------------
!
IF( ln_timing ) CALL timing_start('dyn_zad')
......
src/OCE/DYN/dynzdf.F90
View file @ 672ff863
......@@ -71,7 +71,7 @@ CONTAINS
!!---------------------------------------------------------------------
INTEGER , INTENT( in ) :: kt ! ocean time-step index
INTEGER , INTENT( in ) :: Kbb, Kmm, Krhs, Kaa ! ocean time level indices
REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) :: puu, pvv ! ocean velocities and RHS of momentum equation
REAL(dp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) :: puu, pvv ! ocean velocities and RHS of momentum equation
!
INTEGER :: ji, jj, jk ! dummy loop indices
INTEGER :: iku, ikv ! local integers
......@@ -82,7 +82,7 @@ CONTAINS
REAL(wp) :: zWui, zWvi ! - -
REAL(wp) :: zWus, zWvs ! - -
REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: zwi, zwd, zws ! 3D workspace
REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: ztrdu, ztrdv ! - -
REAL(dp), DIMENSION(:,:,:), ALLOCATABLE :: ztrdu, ztrdv ! - -
!!---------------------------------------------------------------------
!
IF( ln_timing ) CALL timing_start('dyn_zdf')
......
src/OCE/DYN/sshwzv.F90
View file @ 672ff863
......@@ -75,7 +75,7 @@ CONTAINS
!!----------------------------------------------------------------------
INTEGER , INTENT(in ) :: kt ! time step
INTEGER , INTENT(in ) :: Kbb, Kmm, Kaa ! time level index
REAL(wp), DIMENSION(jpi,jpj,jpt), INTENT(inout) :: pssh ! sea-surface height
REAL(dp), DIMENSION(jpi,jpj,jpt), INTENT(inout) :: pssh ! sea-surface height
!
INTEGER :: ji, jj, jk ! dummy loop index
REAL(wp) :: zcoef ! local scalar
......@@ -95,9 +95,6 @@ CONTAINS
! !------------------------------!
! ! After Sea Surface Height !
! !------------------------------!
IF(ln_wd_il) THEN
CALL wad_lmt(pssh(:,:,Kbb), zcoef * (emp_b(:,:) + emp(:,:)), rDt, Kmm, uu, vv )
ENDIF
CALL div_hor( kt, Kbb, Kmm ) ! Horizontal divergence
!
......@@ -113,7 +110,7 @@ CONTAINS
pssh(ji,jj,Kaa) = ( pssh(ji,jj,Kbb) - rDt * ( zcoef * ( emp_b(ji,jj) + emp(ji,jj) ) + zhdiv(ji,jj) ) ) * ssmask(ji,jj)
END_2D
! pssh must be defined everywhere (true for dyn_spg_ts, not for dyn_spg_exp)
IF ( .NOT. ln_dynspg_ts .AND. nn_hls == 2 ) CALL lbc_lnk( 'sshwzv', pssh(:,:,Kaa), 'T', 1.0_wp )
IF ( .NOT. ln_dynspg_ts .AND. nn_hls == 2 ) CALL lbc_lnk( 'sshwzv', pssh(:,:,Kaa), 'T', 1.0_dp )
!
#if defined key_agrif
Kbb_a = Kbb ; Kmm_a = Kmm ; Krhs_a = Kaa
......@@ -122,7 +119,7 @@ CONTAINS
!
IF ( .NOT.ln_dynspg_ts ) THEN
IF( ln_bdy ) THEN
IF (nn_hls==1) CALL lbc_lnk( 'sshwzv', pssh(:,:,Kaa), 'T', 1.0_wp ) ! Not sure that's necessary
IF (nn_hls==1) CALL lbc_lnk( 'sshwzv', pssh(:,:,Kaa), 'T', 1.0_dp ) ! Not sure that's necessary
CALL bdy_ssh( pssh(:,:,Kaa) ) ! Duplicate sea level across open boundaries
ENDIF
ENDIF
......@@ -285,7 +282,7 @@ CONTAINS
!!----------------------------------------------------------------------
INTEGER , INTENT(in ) :: kt ! ocean time-step index
INTEGER , INTENT(in ) :: Kbb, Kmm, Kaa ! ocean time level indices
REAL(wp), DIMENSION(jpi,jpj,jpt), INTENT(inout) :: pssh ! SSH field
REAL(dp), DIMENSION(jpi,jpj,jpt), INTENT(inout) :: pssh ! SSH field
!
REAL(wp) :: zcoef ! local scalar
!!----------------------------------------------------------------------
......
src/OCE/DYN/wet_dry.F90
View file @ 672ff863
......@@ -6,7 +6,7 @@ MODULE wet_dry
!! *** MODULE wet_dry ***
!! Wetting and drying includes initialisation routine and routines to
!! compute and apply flux limiters and preserve water depth positivity
!! only effects if wetting/drying is on (ln_wd_il == .true. or ln_wd_dl==.true. )
!! only effects if wetting/drying is on ( ln_wd_dl==.true. )
!!==============================================================================
!! History : 3.6 ! 2014-09 ((H.Liu) Original code
!! ! will add the runoff and periodic BC case later
......@@ -14,8 +14,6 @@ MODULE wet_dry
!!----------------------------------------------------------------------
!! wad_init : initialisation of wetting and drying
!! wad_lmt : horizontal flux limiter and limited velocity when wetting and drying happens
!! wad_lmt_bt : same as wad_lmt for the barotropic stepping (dynspg_ts)
!!----------------------------------------------------------------------
USE oce ! ocean dynamics and tracers
USE dom_oce ! ocean space and time domain
......@@ -41,7 +39,6 @@ MODULE wet_dry
! ! (can include negative depths)
REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: wdramp, wdrampu, wdrampv !: for hpg limiting
LOGICAL, PUBLIC :: ln_wd_il !: Wetting/drying il activation switch (T:on,F:off)
LOGICAL, PUBLIC :: ln_wd_dl !: Wetting/drying dl activation switch (T:on,F:off)
REAL(wp), PUBLIC :: rn_wdmin0 !: depth at which wetting/drying starts
REAL(wp), PUBLIC :: rn_wdmin1 !: minimum water depth on dried cells
......@@ -50,17 +47,14 @@ MODULE wet_dry
REAL(wp), PUBLIC :: rn_wd_sbcdep !: Depth at which to taper sbc fluxes
REAL(wp), PUBLIC :: rn_wd_sbcfra !: Fraction of SBC at taper depth
REAL(wp), PUBLIC :: rn_wdld !: land elevation below which wetting/drying will be considered
INTEGER , PUBLIC :: nn_wdit !: maximum number of iteration for W/D limiter
LOGICAL, PUBLIC :: ln_wd_dl_bc !: DL scheme: True implies 3D velocities are set to the barotropic values at points
!: where the flow is from wet points on less than half the barotropic sub-steps
LOGICAL, PUBLIC :: ln_wd_dl_rmp !: use a ramp for the dl flux limiter between 2 rn_wdmin1 and rn_wdmin1 (rather than a cut-off at rn_wdmin1)
REAL(wp), PUBLIC :: ssh_ref !: height of z=0 with respect to the geoid;
LOGICAL, PUBLIC :: ll_wd = .FALSE. !: Wetting/drying activation switch (ln_wd_il or ln_wd_dl) <- default def if wad_init not called
LOGICAL, PUBLIC :: ll_wd = .FALSE. !: Wetting/drying activation switch (ln_wd_dl) <- default def if wad_init not called
PUBLIC wad_init ! initialisation routine called by step.F90
PUBLIC wad_lmt ! routine called by sshwzv.F90
PUBLIC wad_lmt_bt ! routine called by dynspg_ts.F90
!! * Substitutions
!!----------------------------------------------------------------------
......@@ -76,8 +70,8 @@ CONTAINS
!!----------------------------------------------------------------------
INTEGER :: ios, ierr ! Local integer
!!
NAMELIST/namwad/ ln_wd_il, ln_wd_dl , rn_wdmin0, rn_wdmin1, rn_wdmin2, rn_wdld, &
& nn_wdit , ln_wd_dl_bc, ln_wd_dl_rmp, rn_wd_sbcdep,rn_wd_sbcfra
NAMELIST/namwad/ ln_wd_dl, rn_wdmin0, rn_wdmin1, rn_wdmin2, rn_wdld, &
& ln_wd_dl_bc, ln_wd_dl_rmp, rn_wd_sbcdep,rn_wd_sbcfra
!!----------------------------------------------------------------------
!
READ ( numnam_ref, namwad, IOSTAT = ios, ERR = 905)
......@@ -92,13 +86,11 @@ CONTAINS
WRITE(numout,*) 'wad_init : Wetting and drying initialization through namelist read'
WRITE(numout,*) '~~~~~~~~'
WRITE(numout,*) ' Namelist namwad'
WRITE(numout,*) ' Logical for Iter Lim wd option ln_wd_il = ', ln_wd_il
WRITE(numout,*) ' Logical for Dir. Lim wd option ln_wd_dl = ', ln_wd_dl
WRITE(numout,*) ' Depth at which wet/drying starts rn_wdmin0 = ', rn_wdmin0
WRITE(numout,*) ' Minimum wet depth on dried cells rn_wdmin1 = ', rn_wdmin1
WRITE(numout,*) ' Tolerance of min wet depth rn_wdmin2 = ', rn_wdmin2
WRITE(numout,*) ' land elevation threshold rn_wdld = ', rn_wdld
WRITE(numout,*) ' Max iteration for W/D limiter nn_wdit = ', nn_wdit
WRITE(numout,*) ' T => baroclinic u,v=0 at dry pts: ln_wd_dl_bc = ', ln_wd_dl_bc
WRITE(numout,*) ' use a ramp for rwd limiter: ln_wd_dl_rwd_rmp = ', ln_wd_dl_rmp
WRITE(numout,*) ' cut off depth sbc for wd rn_wd_sbcdep = ', rn_wd_sbcdep
......@@ -110,285 +102,15 @@ CONTAINS
ENDIF
r_rn_wdmin1 = 1 / rn_wdmin1
IF( ln_wd_il .OR. ln_wd_dl ) THEN
IF( ln_wd_dl ) THEN
ll_wd = .TRUE.
ALLOCATE( wdmask(jpi,jpj), STAT=ierr )
ALLOCATE( wdramp(jpi,jpj), wdrampu(jpi,jpj), wdrampv(jpi,jpj), STAT=ierr )
IF( ierr /= 0 ) CALL ctl_stop('STOP', 'wad_init : Array allocation error')
ENDIF
IF( ln_tile .AND. ln_wd_il ) CALL ctl_warn('Tiling has not been tested with ln_wd_il = T')
!
END SUBROUTINE wad_init
SUBROUTINE wad_lmt( psshb1, psshemp, z2dt, Kmm, puu, pvv )
!!----------------------------------------------------------------------
!! *** ROUTINE wad_lmt ***
!!
!! ** Purpose : generate flux limiters for wetting/drying
!!
!! ** Method : - Prevent negative depth occurring (Not ready for Agrif)
!!
!! ** Action : - calculate flux limiter and W/D flag
!!----------------------------------------------------------------------
REAL(wp), DIMENSION(:,:) , INTENT(inout) :: psshb1
REAL(wp), DIMENSION(:,:) , INTENT(in ) :: psshemp
REAL(wp) , INTENT(in ) :: z2dt
INTEGER , INTENT(in ) :: Kmm ! time level index
REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) :: puu, pvv ! velocity arrays
!
INTEGER :: ji, jj, jk, jk1 ! dummy loop indices
INTEGER :: jflag ! local scalar
REAL(wp) :: zcoef, zdep1, zdep2 ! local scalars
REAL(wp) :: zzflxp, zzflxn ! local scalars
REAL(wp) :: zdepwd ! local scalar, always wet cell depth
REAL(wp) :: ztmp ! local scalars
REAL(wp), DIMENSION(jpi,jpj) :: zwdlmtu, zwdlmtv ! W/D flux limiters
REAL(wp), DIMENSION(jpi,jpj) :: zflxp , zflxn ! local 2D workspace
REAL(wp), DIMENSION(jpi,jpj) :: zflxu , zflxv ! local 2D workspace
REAL(wp), DIMENSION(jpi,jpj) :: zflxu1 , zflxv1 ! local 2D workspace
!!----------------------------------------------------------------------
IF( ln_timing ) CALL timing_start('wad_lmt') !
!
DO jk = 1, jpkm1
puu(:,:,jk,Kmm) = puu(:,:,jk,Kmm) * zwdlmtu(:,:)
pvv(:,:,jk,Kmm) = pvv(:,:,jk,Kmm) * zwdlmtv(:,:)
END DO
jflag = 0
zdepwd = 50._wp ! maximum depth on which that W/D could possibly happen
!
zflxp(:,:) = 0._wp
zflxn(:,:) = 0._wp
zflxu(:,:) = 0._wp
zflxv(:,:) = 0._wp
!
zwdlmtu(:,:) = 1._wp
zwdlmtv(:,:) = 1._wp
!
DO jk = 1, jpkm1 ! Horizontal Flux in u and v direction
zflxu(:,:) = zflxu(:,:) + e3u(:,:,jk,Kmm) * puu(:,:,jk,Kmm) * umask(:,:,jk)
zflxv(:,:) = zflxv(:,:) + e3v(:,:,jk,Kmm) * pvv(:,:,jk,Kmm) * vmask(:,:,jk)
END DO
zflxu(:,:) = zflxu(:,:) * e2u(:,:)
zflxv(:,:) = zflxv(:,:) * e1v(:,:)
!
wdmask(:,:) = 1._wp
DO_2D( 0, 1, 0, 1 )
!
IF( tmask(ji,jj,1) < 0.5_wp ) CYCLE ! we don't care about land cells
IF( ht_0(ji,jj) - ssh_ref > zdepwd ) CYCLE ! and cells which are unlikely to dry
!
zflxp(ji,jj) = MAX( zflxu(ji,jj) , 0._wp ) - MIN( zflxu(ji-1,jj ) , 0._wp ) &
& + MAX( zflxv(ji,jj) , 0._wp ) - MIN( zflxv(ji, jj-1) , 0._wp )
zflxn(ji,jj) = MIN( zflxu(ji,jj) , 0._wp ) - MAX( zflxu(ji-1,jj ) , 0._wp ) &
& + MIN( zflxv(ji,jj) , 0._wp ) - MAX( zflxv(ji, jj-1) , 0._wp )
!
zdep2 = ht_0(ji,jj) + psshb1(ji,jj) - rn_wdmin1
IF( zdep2 <= 0._wp ) THEN ! add more safty, but not necessary
psshb1(ji,jj) = rn_wdmin1 - ht_0(ji,jj)
IF(zflxu(ji, jj) > 0._wp) zwdlmtu(ji ,jj) = 0._wp
IF(zflxu(ji-1,jj) < 0._wp) zwdlmtu(ji-1,jj) = 0._wp
IF(zflxv(ji, jj) > 0._wp) zwdlmtv(ji ,jj) = 0._wp
IF(zflxv(ji,jj-1) < 0._wp) zwdlmtv(ji,jj-1) = 0._wp
wdmask(ji,jj) = 0._wp
END IF
END_2D
!
! ! HPG limiter from jholt
wdramp(:,:) = min((ht_0(:,:) + psshb1(:,:) - rn_wdmin1)/(rn_wdmin0 - rn_wdmin1),1.0_wp)
!jth assume don't need a lbc_lnk here
DO_2D( 1, 0, 1, 0 )
wdrampu(ji,jj) = MIN( wdramp(ji,jj) , wdramp(ji+1,jj) )
wdrampv(ji,jj) = MIN( wdramp(ji,jj) , wdramp(ji,jj+1) )
END_2D
! ! end HPG limiter
!
!
DO jk1 = 1, nn_wdit + 1 !== start limiter iterations ==!
!
zflxu1(:,:) = zflxu(:,:) * zwdlmtu(:,:)
zflxv1(:,:) = zflxv(:,:) * zwdlmtv(:,:)
jflag = 0 ! flag indicating if any further iterations are needed
!
DO_2D( 0, 1, 0, 1 )
IF( tmask(ji, jj, 1) < 0.5_wp ) CYCLE
IF( ht_0(ji,jj) > zdepwd ) CYCLE
!
ztmp = e1e2t(ji,jj)
!
zzflxp = MAX( zflxu1(ji,jj) , 0._wp ) - MIN( zflxu1(ji-1,jj ) , 0._wp) &
& + MAX( zflxv1(ji,jj) , 0._wp ) - MIN( zflxv1(ji, jj-1) , 0._wp)
zzflxn = MIN( zflxu1(ji,jj) , 0._wp ) - MAX( zflxu1(ji-1,jj ) , 0._wp) &
& + MIN( zflxv1(ji,jj) , 0._wp ) - MAX( zflxv1(ji, jj-1) , 0._wp)
!
zdep1 = (zzflxp + zzflxn) * z2dt / ztmp
zdep2 = ht_0(ji,jj) + psshb1(ji,jj) - rn_wdmin1 - z2dt * psshemp(ji,jj)
!
IF( zdep1 > zdep2 ) THEN
wdmask(ji, jj) = 0._wp
zcoef = ( ( zdep2 - rn_wdmin2 ) * ztmp - zzflxn * z2dt ) / ( zflxp(ji,jj) * z2dt )
!zcoef = ( ( zdep2 - rn_wdmin2 ) * ztmp - zzflxn * z2dt ) / ( zzflxp * z2dt )
! flag if the limiter has been used but stop flagging if the only
! changes have zeroed the coefficient since further iterations will
! not change anything
IF( zcoef > 0._wp ) THEN ; jflag = 1
ELSE ; zcoef = 0._wp
ENDIF
IF( jk1 > nn_wdit ) zcoef = 0._wp
IF( zflxu1(ji ,jj ) > 0._wp ) zwdlmtu(ji ,jj ) = zcoef
IF( zflxu1(ji-1,jj ) < 0._wp ) zwdlmtu(ji-1,jj ) = zcoef
IF( zflxv1(ji ,jj ) > 0._wp ) zwdlmtv(ji ,jj ) = zcoef
IF( zflxv1(ji ,jj-1) < 0._wp ) zwdlmtv(ji ,jj-1) = zcoef
ENDIF
END_2D
CALL lbc_lnk( 'wet_dry', zwdlmtu, 'U', 1.0_wp, zwdlmtv, 'V', 1.0_wp )
!
CALL mpp_max('wet_dry', jflag) !max over the global domain
!
IF( jflag == 0 ) EXIT
!
END DO ! jk1 loop
!
DO jk = 1, jpkm1
puu(:,:,jk,Kmm) = puu(:,:,jk,Kmm) * zwdlmtu(:,:)
pvv(:,:,jk,Kmm) = pvv(:,:,jk,Kmm) * zwdlmtv(:,:)
END DO
uu_b(:,:,Kmm) = uu_b(:,:,Kmm) * zwdlmtu(:, :)
vv_b(:,:,Kmm) = vv_b(:,:,Kmm) * zwdlmtv(:, :)
!
!!gm TO BE SUPPRESSED ? these lbc_lnk are useless since zwdlmtu and zwdlmtv are defined everywhere !
CALL lbc_lnk( 'wet_dry', puu(:,:,:,Kmm) , 'U', -1.0_wp, pvv(:,:,:,Kmm) , 'V', -1.0_wp )
CALL lbc_lnk( 'wet_dry', uu_b(:,:,Kmm), 'U', -1.0_wp, vv_b(:,:,Kmm), 'V', -1.0_wp )
!!gm
!
IF(jflag == 1 .AND. lwp) WRITE(numout,*) 'Need more iterations in wad_lmt!!!'
!
!IF( ln_rnf ) CALL sbc_rnf_div( hdiv ) ! runoffs (update hdiv field)
!
IF( ln_timing ) CALL timing_stop('wad_lmt') !
!
END SUBROUTINE wad_lmt
SUBROUTINE wad_lmt_bt( zflxu, zflxv, sshn_e, zssh_frc, rDt_e )
!!----------------------------------------------------------------------
!! *** ROUTINE wad_lmt ***
!!
!! ** Purpose : limiting flux in the barotropic stepping (dynspg_ts)
!!
!! ** Method : - Prevent negative depth occurring (Not ready for Agrif)
!!
!! ** Action : - calculate flux limiter and W/D flag
!!----------------------------------------------------------------------
REAL(wp) , INTENT(in ) :: rDt_e ! ocean time-step index
REAL(wp), DIMENSION(:,:), INTENT(inout) :: zflxu, zflxv, sshn_e, zssh_frc
!
INTEGER :: ji, jj, jk, jk1 ! dummy loop indices
INTEGER :: jflag ! local integer
REAL(wp) :: z2dt
REAL(wp) :: zcoef, zdep1, zdep2 ! local scalars
REAL(wp) :: zzflxp, zzflxn ! local scalars
REAL(wp) :: zdepwd ! local scalar, always wet cell depth
REAL(wp) :: ztmp ! local scalars
REAL(wp), DIMENSION(jpi,jpj) :: zwdlmtu, zwdlmtv !: W/D flux limiters
REAL(wp), DIMENSION(jpi,jpj) :: zflxp, zflxn ! local 2D workspace
REAL(wp), DIMENSION(jpi,jpj) :: zflxu1, zflxv1 ! local 2D workspace
!!----------------------------------------------------------------------
IF( ln_timing ) CALL timing_start('wad_lmt_bt') !
!
jflag = 0
zdepwd = 50._wp ! maximum depth that ocean cells can have W/D processes
!
z2dt = rDt_e
!
zflxp(:,:) = 0._wp
zflxn(:,:) = 0._wp
zwdlmtu(:,:) = 1._wp
zwdlmtv(:,:) = 1._wp
!
DO_2D( 0, 1, 0, 1 ) ! Horizontal Flux in u and v direction
!
IF( tmask(ji, jj, 1 ) < 0.5_wp) CYCLE ! we don't care about land cells
IF( ht_0(ji,jj) > zdepwd ) CYCLE ! and cells which are unlikely to dry
!
zflxp(ji,jj) = MAX( zflxu(ji,jj) , 0._wp ) - MIN( zflxu(ji-1,jj ) , 0._wp ) &
& + MAX( zflxv(ji,jj) , 0._wp ) - MIN( zflxv(ji, jj-1) , 0._wp )
zflxn(ji,jj) = MIN( zflxu(ji,jj) , 0._wp ) - MAX( zflxu(ji-1,jj ) , 0._wp ) &
& + MIN( zflxv(ji,jj) , 0._wp ) - MAX( zflxv(ji, jj-1) , 0._wp )
!
zdep2 = ht_0(ji,jj) + sshn_e(ji,jj) - rn_wdmin1
IF( zdep2 <= 0._wp ) THEN !add more safety, but not necessary
sshn_e(ji,jj) = rn_wdmin1 - ht_0(ji,jj)
IF( zflxu(ji ,jj ) > 0._wp) zwdlmtu(ji ,jj ) = 0._wp
IF( zflxu(ji-1,jj ) < 0._wp) zwdlmtu(ji-1,jj ) = 0._wp
IF( zflxv(ji ,jj ) > 0._wp) zwdlmtv(ji ,jj ) = 0._wp
IF( zflxv(ji ,jj-1) < 0._wp) zwdlmtv(ji ,jj-1) = 0._wp
ENDIF
END_2D
!
DO jk1 = 1, nn_wdit + 1 !! start limiter iterations
!
zflxu1(:,:) = zflxu(:,:) * zwdlmtu(:,:)
zflxv1(:,:) = zflxv(:,:) * zwdlmtv(:,:)
jflag = 0 ! flag indicating if any further iterations are needed
!
DO_2D( 0, 1, 0, 1 )
!
IF( tmask(ji, jj, 1 ) < 0.5_wp ) CYCLE
IF( ht_0(ji,jj) > zdepwd ) CYCLE
!
ztmp = e1e2t(ji,jj)
!
zzflxp = max(zflxu1(ji,jj), 0._wp) - min(zflxu1(ji-1,jj), 0._wp) &
& + max(zflxv1(ji,jj), 0._wp) - min(zflxv1(ji, jj-1), 0._wp)
zzflxn = min(zflxu1(ji,jj), 0._wp) - max(zflxu1(ji-1,jj), 0._wp) &
& + min(zflxv1(ji,jj), 0._wp) - max(zflxv1(ji, jj-1), 0._wp)
zdep1 = (zzflxp + zzflxn) * z2dt / ztmp
zdep2 = ht_0(ji,jj) + sshn_e(ji,jj) - rn_wdmin1 - z2dt * zssh_frc(ji,jj)
IF(zdep1 > zdep2) THEN
zcoef = ( ( zdep2 - rn_wdmin2 ) * ztmp - zzflxn * z2dt ) / ( zflxp(ji,jj) * z2dt )
!zcoef = ( ( zdep2 - rn_wdmin2 ) * ztmp - zzflxn * z2dt ) / ( zzflxp * z2dt )
! flag if the limiter has been used but stop flagging if the only
! changes have zeroed the coefficient since further iterations will
! not change anything
IF( zcoef > 0._wp ) THEN
jflag = 1
ELSE
zcoef = 0._wp
ENDIF
IF(jk1 > nn_wdit) zcoef = 0._wp
IF(zflxu1(ji, jj) > 0._wp) zwdlmtu(ji ,jj) = zcoef
IF(zflxu1(ji-1,jj) < 0._wp) zwdlmtu(ji-1,jj) = zcoef
IF(zflxv1(ji, jj) > 0._wp) zwdlmtv(ji ,jj) = zcoef
IF(zflxv1(ji,jj-1) < 0._wp) zwdlmtv(ji,jj-1) = zcoef
END IF
END_2D
!
CALL lbc_lnk( 'wet_dry', zwdlmtu, 'U', 1.0_wp, zwdlmtv, 'V', 1.0_wp )
!
CALL mpp_max('wet_dry', jflag) !max over the global domain
!
IF(jflag == 0) EXIT
!
END DO ! jk1 loop
!
zflxu(:,:) = zflxu(:,:) * zwdlmtu(:, :)
zflxv(:,:) = zflxv(:,:) * zwdlmtv(:, :)
!
!!gm THIS lbc_lnk is useless since it is already done at the end of the jk1-loop
CALL lbc_lnk( 'wet_dry', zflxu, 'U', -1.0_wp, zflxv, 'V', -1.0_wp )
!!gm end
!
IF( jflag == 1 .AND. lwp ) WRITE(numout,*) 'Need more iterations in wad_lmt_bt!!!'
!
!IF( ln_rnf ) CALL sbc_rnf_div( hdiv ) ! runoffs (update hdiv field)
!
IF( ln_timing ) CALL timing_stop('wad_lmt_bt') !
!
END SUBROUTINE wad_lmt_bt
!!==============================================================================
END MODULE wet_dry
src/OCE/FLO/floblk.F90
View file @ 672ff863
......@@ -383,4 +383,4 @@ CONTAINS
END SUBROUTINE flo_blk
!!======================================================================
END MODULE floblk
END MODULE floblk
src/OCE/ICB/icbini.F90
View file @ 672ff863
......@@ -522,6 +522,16 @@ CONTAINS
CALL ctl_stop( 'icb_nam: a negative rn_distribution value encountered ==>> change your namelist namberg' )
ENDIF
!
! ensure that nn_verbose_write is a multiple of nn_fsbc
IF (MOD(nn_verbose_write, nn_fsbc) /= 0) THEN
CALL ctl_stop( 'icb_nam: nn_verbose_write is not a multiple of nn_fsbc')
END IF
!
! ensure that nn_sample_rate is a multiple of nn_fsbc
IF (MOD(nn_sample_rate, nn_fsbc) /= 0) THEN
CALL ctl_stop( 'icb_nam: nn_sample_rate is not a multiple of nn_fsbc')
END IF
!
END SUBROUTINE icb_nam
!!======================================================================
......
src/OCE/ICB/icbrst.F90
View file @ 672ff863
......@@ -24,6 +24,7 @@ MODULE icbrst
USE lib_mpp ! NEMO MPI library, lk_mpp in particular
USE netcdf ! netcdf routines for IO
USE iom
USE ioipsl, ONLY : ju2ymds ! for calendar
USE icb_oce ! define iceberg arrays
USE icbutl ! iceberg utility routines
......@@ -190,9 +191,12 @@ CONTAINS
INTEGER :: jn ! dummy loop index
INTEGER :: idg ! number of digits
INTEGER :: ix_dim, iy_dim, ik_dim, in_dim
INTEGER :: iyear, imonth, iday
REAL (wp) :: zsec
REAL (wp) :: zfjulday
CHARACTER(len=256) :: cl_path
CHARACTER(len=256) :: cl_filename
CHARACTER(len=8 ) :: cl_kt
CHARACTER(LEN=20) :: cl_kt ! ocean time-step deine as a character
CHARACTER(LEN=12 ) :: clfmt ! writing format
TYPE(iceberg), POINTER :: this
TYPE(point) , POINTER :: pt
......@@ -212,8 +216,17 @@ CONTAINS
IF( cl_path(LEN_TRIM(cl_path):) /= '/' ) cl_path = TRIM(cl_path) // '/'
!
! file name
WRITE(cl_kt, '(i8.8)') kt
cl_filename = TRIM(cexper)//"_"//cl_kt//"_"//TRIM(cn_icbrst_out)
IF ( ln_rstdate ) THEN
zfjulday = fjulday + rdt / rday
IF( ABS(zfjulday - REAL(NINT(zfjulday),wp)) < 0.1 / rday ) zfjulday = REAL(NINT(zfjulday),wp) ! avoid truncation error
CALL ju2ymds( zfjulday, iyear, imonth, iday, zsec )
WRITE(cl_kt, '(i4.4,2i2.2)') iyear, imonth, iday
ELSE
IF( kt > 999999999 ) THEN ; WRITE(cl_kt, * ) kt
ELSE ; WRITE(cl_kt, '(i8.8)') kt
ENDIF
ENDIF
cl_filename = TRIM(cexper)//"_"//TRIM(cl_kt)//"_"//TRIM(cn_icbrst_out)
IF( lk_mpp ) THEN
idg = MAX( INT(LOG10(REAL(MAX(1,jpnij-1),wp))) + 1, 4 ) ! how many digits to we need to write? min=4, max=9
WRITE(clfmt, "('(a,a,i', i1, '.', i1, ',a)')") idg, idg ! '(a,a,ix.x,a)'
......
src/OCE/ICB/icbstp.F90
View file @ 672ff863
......@@ -68,78 +68,79 @@ CONTAINS
IF( ln_timing ) CALL timing_start('icb_stp')
! !== start of timestep housekeeping ==!
!
nktberg = kt
!
!CALL test_icb_utl_getkb
!CALL ctl_stop('end test icb')
!
IF( nn_test_icebergs < 0 .OR. ln_use_calving ) THEN !* read calving data
IF( MOD( kt-1, nn_fsbc ) == 0 ) THEN
!
CALL fld_read ( kt, 1, sf_icb )
src_calving (:,:) = sf_icb(1)%fnow(:,:,1) ! calving in km^3/year (water equivalent)
src_calving_hflx(:,:) = 0._wp ! NO heat flux for now
nktberg = kt
!
ENDIF
!
berg_grid%floating_melt(:,:) = 0._wp
!
! !* anything that needs to be reset to zero each timestep
CALL icb_dia_step() ! for budgets is dealt with here
!
! !* write out time
ll_verbose = .FALSE.
IF( nn_verbose_write > 0 .AND. MOD( kt-1 , nn_verbose_write ) == 0 ) ll_verbose = ( nn_verbose_level > 0 )
!
IF( ll_verbose ) WRITE(numicb,9100) nktberg, ndastp, nsec_day
9100 FORMAT('kt= ',i8, ' day= ',i8,' secs=',i8)
!
! !* copy nemo forcing arrays into iceberg versions with extra halo
CALL icb_utl_copy( Kmm ) ! only necessary for variables not on T points
!
!
! !== process icebergs ==!
! !
CALL icb_clv_flx( kt ) ! Accumulate ice from calving
! !
CALL icb_clv( kt ) ! Calve excess stored ice into icebergs
! !
!
! !== For each berg, evolve ==!
!
IF( ASSOCIATED(first_berg) ) CALL icb_dyn( kt ) ! ice berg dynamics
IF( lk_mpp ) THEN ; CALL icb_lbc_mpp() ! Send bergs to other PEs
ELSE ; CALL icb_lbc() ! Deal with any cyclic boundaries in non-mpp case
ENDIF
IF( ASSOCIATED(first_berg) ) CALL icb_thm( kt ) ! Ice berg thermodynamics (melting) + rolling
!
!
! !== diagnostics and output ==!
!
! !* For each berg, record trajectory (when needed)
ll_sample_traj = .FALSE.
IF( nn_sample_rate > 0 .AND. MOD(kt-1,nn_sample_rate) == 0 ) ll_sample_traj = .TRUE.
IF( ll_sample_traj .AND. ASSOCIATED(first_berg) ) CALL icb_trj_write( kt )
! !* Gridded diagnostics
! ! To get these iom_put's and those preceding to actually do something
! ! use key_xios in cpp file and create content for XML file
!
CALL iom_put( "calving" , berg_grid%calving (:,:) ) ! 'calving mass input'
CALL iom_put( "berg_floating_melt", berg_grid%floating_melt(:,:) ) ! 'Melt rate of icebergs + bits' , 'kg/m2/s'
CALL iom_put( "berg_stored_ice" , berg_grid%stored_ice (:,:,:) ) ! 'Accumulated ice mass by class', 'kg'
!
CALL icb_dia_put() !* store mean budgets
!
! !* Dump icebergs to screen
IF( nn_verbose_level >= 2 ) CALL icb_utl_print( 'icb_stp, status', kt )
!
! !* Diagnose budgets
ll_budget = .FALSE.
IF( nn_verbose_write > 0 .AND. MOD(kt-1,nn_verbose_write) == 0 ) ll_budget = ln_bergdia
CALL icb_dia( ll_budget )
IF( nn_test_icebergs < 0 .OR. ln_use_calving ) THEN !* read calving data
!
CALL fld_read ( kt, 1, sf_icb )
src_calving (:,:) = sf_icb(1)%fnow(:,:,1) ! calving in km^3/year (water equivalent)
src_calving_hflx(:,:) = 0._wp ! NO heat flux for now
!
ENDIF
!
berg_grid%floating_melt(:,:) = 0._wp
!
! !* anything that needs to be reset to zero each timestep
CALL icb_dia_step() ! for budgets is dealt with here
!
! !* write out time
ll_verbose = .FALSE.
IF( nn_verbose_write > 0 .AND. MOD( kt-1 , nn_verbose_write ) == 0 ) ll_verbose = ( nn_verbose_level > 0 )
!
IF( ll_verbose ) WRITE(numicb,9100) nktberg, ndastp, nsec_day
9100 FORMAT('kt= ',i8, ' day= ',i8,' secs=',i8)
!
! !* copy nemo forcing arrays into iceberg versions with extra halo
CALL icb_utl_copy( Kmm ) ! only necessary for variables not on T points
!
!
! !== process icebergs ==!
! !
CALL icb_clv_flx( kt ) ! Accumulate ice from calving
! !
CALL icb_clv( kt ) ! Calve excess stored ice into icebergs
! !
!
! !== For each berg, evolve ==!
!
IF( ASSOCIATED(first_berg) ) CALL icb_dyn( kt ) ! ice berg dynamics
IF( lk_mpp ) THEN ; CALL icb_lbc_mpp() ! Send bergs to other PEs
ELSE ; CALL icb_lbc() ! Deal with any cyclic boundaries in non-mpp case
ENDIF
IF( ASSOCIATED(first_berg) ) CALL icb_thm( kt ) ! Ice berg thermodynamics (melting) + rolling
!
!
! !== diagnostics and output ==!
!
! !* For each berg, record trajectory (when needed)
ll_sample_traj = .FALSE.
IF( nn_sample_rate > 0 .AND. MOD(kt-1,nn_sample_rate) == 0 ) ll_sample_traj = .TRUE.
IF( ll_sample_traj .AND. ASSOCIATED(first_berg) ) CALL icb_trj_write( kt )
! !* Gridded diagnostics
! ! To get these iom_put's and those preceding to actually do something
! ! use key_xios in cpp file and create content for XML file
!
CALL iom_put( "calving" , berg_grid%calving (:,:) ) ! 'calving mass input'
CALL iom_put( "berg_floating_melt", berg_grid%floating_melt(:,:) ) ! 'Melt rate of icebergs + bits' , 'kg/m2/s'
CALL iom_put( "berg_stored_ice" , berg_grid%stored_ice (:,:,:) ) ! 'Accumulated ice mass by class', 'kg'
!
CALL icb_dia_put() !* store mean budgets
!
! !* Dump icebergs to screen
IF( nn_verbose_level >= 2 ) CALL icb_utl_print( 'icb_stp, status', kt )
!
! !* Diagnose budgets
ll_budget = .FALSE.
IF( nn_verbose_write > 0 .AND. MOD(kt-1,nn_verbose_write) == 0 ) ll_budget = ln_bergdia
CALL icb_dia( ll_budget )
!
END IF
!
IF( lrst_oce ) THEN !* restart
CALL icb_rst_write( kt )
......
src/OCE/ICB/icbtrj.F90
View file @ 672ff863
......@@ -62,7 +62,7 @@ CONTAINS
!
INTEGER :: iret, iyear, imonth, iday
INTEGER :: idg ! number of digits
REAL(wp) :: zfjulday, zsec
REAL(dp) :: zfjulday, zsec
CHARACTER(len=80) :: cl_filename
CHARACTER(LEN=12) :: clfmt ! writing format
CHARACTER(LEN=8 ) :: cldate_ini, cldate_end
......
src/OCE/ICB/icbutl.F90
View file @ 672ff863
......@@ -57,6 +57,7 @@ MODULE icbutl
PUBLIC icb_utl_heat ! routine called in icbdia module
!! * Substitutions
# include "single_precision_substitute.h90"
# include "domzgr_substitute.h90"
!!----------------------------------------------------------------------
!! NEMO/OCE 4.0 , NEMO Consortium (2018)
......@@ -181,8 +182,8 @@ CONTAINS
CALL icb_utl_pos( pi, pj, 'F', iiF, ijF, zwF, zmskF )
!
! metrics and coordinates
IF ( PRESENT(pe1 ) ) pe1 = icb_utl_bilin_e( e1t, e1u, e1v, e1f, pi, pj ) ! scale factors
IF ( PRESENT(pe2 ) ) pe2 = icb_utl_bilin_e( e2t, e2u, e2v, e2f, pi, pj )
IF ( PRESENT(pe1 ) ) pe1 = icb_utl_bilin_e( e1t, CASTSP(e1u), e1v, e1f, pi, pj ) ! scale factors
IF ( PRESENT(pe2 ) ) pe2 = icb_utl_bilin_e( e2t, e2u, CASTSP(e2v), e2f, pi, pj )
IF ( PRESENT(plon) ) plon= icb_utl_bilin_h( rlon_e, iiT, ijT, zwT, .true. )
IF ( PRESENT(plat) ) plat= icb_utl_bilin_h( rlat_e, iiT, ijT, zwT, .false. )
!
......@@ -214,17 +215,17 @@ CONTAINS
!
! Estimate SSH gradient in i- and j-direction (centred evaluation)
IF ( PRESENT(pssh_i) .AND. PRESENT(pssh_j) ) THEN
CALL icb_utl_pos( pi+0.1, pj , 'T', iiTp, ijTp, zwTp, zmskTp )
CALL icb_utl_pos( pi-0.1, pj , 'T', iiTm, ijTm, zwTm, zmskTm )
CALL icb_utl_pos( pi+0.1_wp, pj , 'T', iiTp, ijTp, zwTp, zmskTp )
CALL icb_utl_pos( pi-0.1_wp, pj , 'T', iiTm, ijTm, zwTm, zmskTm )
!
IF ( .NOT. PRESENT(pe1) ) pe1 = icb_utl_bilin_e( e1t, e1u, e1v, e1f, pi, pj )
IF ( .NOT. PRESENT(pe1) ) pe1 = icb_utl_bilin_e( e1t, CASTSP(e1u), e1v, e1f, pi, pj )
pssh_i = ( icb_utl_bilin_h( ssh_e, iiTp, ijTp, zwTp*zmskTp, .false. ) - &
& icb_utl_bilin_h( ssh_e, iiTm, ijTm, zwTm*zmskTm, .false. ) ) / ( 0.2_wp * pe1 )
!
CALL icb_utl_pos( pi , pj+0.1, 'T', iiTp, ijTp, zwTp, zmskTp )
CALL icb_utl_pos( pi , pj-0.1, 'T', iiTm, ijTm, zwTm, zmskTm )
CALL icb_utl_pos( pi , pj+0.1_wp, 'T', iiTp, ijTp, zwTp, zmskTp )
CALL icb_utl_pos( pi , pj-0.1_wp, 'T', iiTm, ijTm, zwTm, zmskTm )
!
IF ( .NOT. PRESENT(pe2) ) pe2 = icb_utl_bilin_e( e2t, e2u, e2v, e2f, pi, pj )
IF ( .NOT. PRESENT(pe2) ) pe2 = icb_utl_bilin_e( e2t, e2u, CASTSP(e2v), e2f, pi, pj )
pssh_j = ( icb_utl_bilin_h( ssh_e, iiTp, ijTp, zwTp*zmskTp, .false. ) - &
& icb_utl_bilin_h( ssh_e, iiTm, ijTm, zwTm*zmskTm, .false. ) ) / ( 0.2_wp * pe2 )
END IF
......@@ -437,7 +438,8 @@ CONTAINS
!! ** Method : interpolation done using the 4 nearest grid points among
!! t-, u-, v-, and f-points.
!!----------------------------------------------------------------------
REAL(wp), DIMENSION(:,:), INTENT(in) :: pet, peu, pev, pef ! horizontal scale factor to be interpolated at t-,u-,v- & f-pts
REAL(wp), DIMENSION(:,:), INTENT(in) :: peu, pev! horizontal scale factor to be interpolated at t-,u-,v- & f-pts
REAL(dp), DIMENSION(:,:), INTENT(in) :: pet, pef! horizontal scale factor to be interpolated at t-,u-,v- & f-pts
REAL(wp) , INTENT(IN) :: pi , pj ! iceberg position
!
! weights corresponding to corner points of a T cell quadrant
......@@ -521,7 +523,7 @@ CONTAINS
kb = kb + 1
END DO
kb = MIN(kb - 1,jpk)
END SUBROUTINE
END SUBROUTINE icb_utl_getkb
SUBROUTINE icb_utl_zavg(pzavg, pdat, pe3, pD, kb )
!!----------------------------------------------------------------------
......@@ -546,7 +548,7 @@ CONTAINS
! if kb is limited by mbkt => bottom value is used between bathy and icb tail
! if kb not limited by mbkt => ocean value over mask is used (ie 0.0 for u, v)
pzavg = ( pzavg + (pD - zdep)*pdat(kb)) / pD
END SUBROUTINE
END SUBROUTINE icb_utl_zavg
SUBROUTINE icb_utl_add( bergvals, ptvals )
!!----------------------------------------------------------------------
......@@ -750,7 +752,7 @@ CONTAINS
!!
!!----------------------------------------------------------------------
CHARACTER(len=*) :: cd_label
INTEGER :: kt ! timestep number
INTEGER, INTENT(IN) :: kt ! timestep number
!
INTEGER :: ibergs, inbergs
TYPE(iceberg), POINTER :: this
......@@ -916,7 +918,8 @@ CONTAINS
!! ** Comments : not called, if needed a CALL test_icb_utl_getkb need to be added in icb_step
!!----------------------------------------------------------------------
INTEGER :: ikb
REAL(wp) :: zD, zout
REAL(wp) :: zout
REAL(wp) :: zD
REAL(wp), DIMENSION(jpk) :: ze3, zin
WRITE(numout,*) 'Test icb_utl_getkb : '
zD = 0.0 ; ze3= 20.0
......
src/OCE/IOM/in_out_manager.F90
View file @ 672ff863
......@@ -27,6 +27,7 @@ MODULE in_out_manager
CHARACTER(lc) :: cn_ocerst_outdir !: restart output directory
LOGICAL :: ln_rstart !: start from (F) rest or (T) a restart file
LOGICAL :: ln_rst_list !: output restarts at list of times (T) or by frequency (F)
LOGICAL :: ln_rst_eos !: check equation of state used for the restart is consistent with model
INTEGER :: nn_rstctl !: control of the time step (0, 1 or 2)
INTEGER :: nn_rstssh = 0 !: hand made initilization of ssh or not (1/0)
INTEGER :: nn_it000 !: index of the first time step
......@@ -39,6 +40,7 @@ MODULE in_out_manager
INTEGER :: nn_stock !: restart file frequency
INTEGER, DIMENSION(10) :: nn_stocklist !: restart dump times
LOGICAL :: ln_mskland !: mask land points in NetCDF outputs (costly: + ~15%)
LOGICAL :: ln_rstdate !: T=> stamp output restart files with date instead of timestep
LOGICAL :: ln_cfmeta !: output additional data to netCDF files required for compliance with the CF metadata standard
LOGICAL :: ln_clobber !: clobber (overwrite) an existing file
INTEGER :: nn_chunksz !: chunksize (bytes) for NetCDF file (works only with iom_nf90 routines)
......
src/OCE/IOM/iom.F90
View file @ 672ff863
......@@ -263,7 +263,9 @@ CONTAINS
# if defined key_si3
CALL iom_set_axis_attr( "ncatice", (/ (REAL(ji,wp), ji=1,jpl) /) )
! SIMIP diagnostics (4 main arctic straits)
CALL iom_set_axis_attr( "nstrait", (/ (REAL(ji,wp), ji=1,4) /) )
! remove from si3 case to general case to use dct diagnostics over
! oce+ice
! CALL iom_set_axis_attr( "nstrait", (/ (REAL(ji,wp), ji=1,4) /) )
# endif
#if defined key_top
IF( ALLOCATED(profsed) ) CALL iom_set_axis_attr( "profsed", paxis = profsed )
......@@ -276,6 +278,12 @@ CONTAINS
INQUIRE( FILE = 'subbasins.nc', EXIST = ll_exist )
nbasin = 1 + 4 * COUNT( (/ll_exist/) )
CALL iom_set_axis_attr( "basin" , (/ (REAL(ji,wp), ji=1,nbasin) /) )
! Transport diagnostics diadct - max number of section 150 in diadct -> ! pb to ensure consistency here
IF( ln_diadct ) THEN
CALL iom_set_axis_attr( "nstrait" , (/ (REAL(ji,wp), ji=1,nb_sec) /) )
ELSE
CALL iom_set_axis_attr( "nstrait", (/ (REAL(ji,wp), ji=1,4) /) )
ENDIF
ENDIF
!
! automatic definitions of some of the xml attributs
......@@ -1942,7 +1950,9 @@ CONTAINS
IF( iom_use(cdname) ) THEN
#if defined key_xios
IF( is_tile(pfield2d) == 1 ) THEN
#if ! defined key_xios3
CALL xios_send_field( cdname, pfield2d, ntile - 1 )
#endif
ELSE IF( .NOT. l_istiled .OR. ntile == nijtile ) THEN
CALL xios_send_field( cdname, pfield2d )
ENDIF
......@@ -1958,7 +1968,9 @@ CONTAINS
IF( iom_use(cdname) ) THEN
#if defined key_xios
IF( is_tile(pfield2d) == 1 ) THEN
#if ! defined key_xios3
CALL xios_send_field( cdname, pfield2d, ntile - 1 )
#endif
ELSE IF( .NOT. l_istiled .OR. ntile == nijtile ) THEN
CALL xios_send_field( cdname, pfield2d )
ENDIF
......@@ -1974,7 +1986,9 @@ CONTAINS
IF( iom_use(cdname) ) THEN
#if defined key_xios
IF( is_tile(pfield3d) == 1 ) THEN
#if ! defined key_xios3
CALL xios_send_field( cdname, pfield3d, ntile - 1 )
#endif
ELSE IF( .NOT. l_istiled .OR. ntile == nijtile ) THEN
CALL xios_send_field( cdname, pfield3d )
ENDIF
......@@ -1990,7 +2004,9 @@ CONTAINS
IF( iom_use(cdname) ) THEN
#if defined key_xios
IF( is_tile(pfield3d) == 1 ) THEN
#if ! defined key_xios3
CALL xios_send_field( cdname, pfield3d, ntile - 1 )
#endif
ELSE IF( .NOT. l_istiled .OR. ntile == nijtile ) THEN
CALL xios_send_field( cdname, pfield3d )
ENDIF
......@@ -2006,7 +2022,9 @@ CONTAINS
IF( iom_use(cdname) ) THEN
#if defined key_xios
IF( is_tile(pfield4d) == 1 ) THEN
#if ! defined key_xios3
CALL xios_send_field( cdname, pfield4d, ntile - 1 )
#endif
ELSE IF( .NOT. l_istiled .OR. ntile == nijtile ) THEN
CALL xios_send_field( cdname, pfield4d )
ENDIF
......@@ -2022,7 +2040,9 @@ CONTAINS
IF( iom_use(cdname) ) THEN
#if defined key_xios
IF( is_tile(pfield4d) == 1 ) THEN
#if ! defined key_xios3
CALL xios_send_field( cdname, pfield4d, ntile - 1 )
#endif
ELSE IF( .NOT. l_istiled .OR. ntile == nijtile ) THEN
CALL xios_send_field( cdname, pfield4d )
ENDIF
......@@ -2058,18 +2078,22 @@ CONTAINS
IF( xios_is_valid_domain (cdid) ) THEN
CALL xios_set_domain_attr ( cdid, ni_glo=ni_glo, nj_glo=nj_glo, ibegin=ibegin, jbegin=jbegin, ni=ni, nj=nj, &
& data_dim=data_dim, data_ibegin=data_ibegin, data_ni=data_ni, data_jbegin=data_jbegin, data_nj=data_nj , &
#if ! defined key_xios3
& ntiles=ntiles, tile_ibegin=tile_ibegin, tile_jbegin=tile_jbegin, tile_ni=tile_ni, tile_nj=tile_nj, &
& tile_data_ibegin=tile_data_ibegin, tile_data_jbegin=tile_data_jbegin, &
& tile_data_ni=tile_data_ni, tile_data_nj=tile_data_nj, &
#endif
& lonvalue_1D=lonvalue, latvalue_1D=latvalue, mask_1D=mask, nvertex=nvertex, bounds_lon_1D=bounds_lon, &
& bounds_lat_1D=bounds_lat, area=area, type='curvilinear')
ENDIF
IF( xios_is_valid_domaingroup(cdid) ) THEN
CALL xios_set_domaingroup_attr( cdid, ni_glo=ni_glo, nj_glo=nj_glo, ibegin=ibegin, jbegin=jbegin, ni=ni, nj=nj, &
& data_dim=data_dim, data_ibegin=data_ibegin, data_ni=data_ni, data_jbegin=data_jbegin, data_nj=data_nj , &
#if ! defined key_xios3
& ntiles=ntiles, tile_ibegin=tile_ibegin, tile_jbegin=tile_jbegin, tile_ni=tile_ni, tile_nj=tile_nj, &
& tile_data_ibegin=tile_data_ibegin, tile_data_jbegin=tile_data_jbegin, &
& tile_data_ni=tile_data_ni, tile_data_nj=tile_data_nj, &
#endif
& lonvalue_1D=lonvalue, latvalue_1D=latvalue, mask_1D=mask, nvertex=nvertex, bounds_lon_1D=bounds_lon, &
& bounds_lat_1D=bounds_lat, area=area, type='curvilinear' )
ENDIF
......@@ -2250,9 +2274,11 @@ CONTAINS
LOGICAL, INTENT(IN) :: ldxios, ldrxios
!!----------------------------------------------------------------------
!
! nn_hls halo points
CALL iom_set_domain_attr("grid_"//cdgrd, ni_glo=Ni0glo,nj_glo=Nj0glo,ibegin=mig0(Nis0)-1,jbegin=mjg0(Njs0)-1,ni=Ni_0,nj=Nj_0)
CALL iom_set_domain_attr("grid_"//cdgrd, data_dim=2, data_ibegin = -nn_hls, data_ni=jpi, data_jbegin = -nn_hls, data_nj=jpj)
! Inner domain only
CALL iom_set_domain_attr("grid_"//cdgrd//"_inner", ni_glo = Ni0glo, nj_glo = Nj0glo, &
& ibegin = mig0(Nis0) - 1, jbegin = mjg0(Njs0) - 1, ni = Ni_0, nj = Nj_0)
CALL iom_set_domain_attr("grid_"//cdgrd//"_inner", data_dim=2, data_ibegin = 0, data_ni=Ni_0, data_jbegin = 0, data_nj=Nj_0)
......@@ -2264,14 +2290,16 @@ CONTAINS
idb(jn) = -nn_hls ! Tile data offset (halo size)
END DO
! Tile_[ij]begin are defined with respect to the processor data domain, so data_[ij]begin is added
! Data includes all halo points
CALL iom_set_domain_attr("grid_"//cdgrd, ntiles=nijtile, &
& tile_ibegin=ntsi_a(1:nijtile) + idb(:) - 1, tile_jbegin=ntsj_a(1:nijtile) + idb(:) - 1, &
& tile_ibegin=ntsi_a(1:nijtile) - nn_hls - 1, tile_jbegin=ntsj_a(1:nijtile) - nn_hls - 1, &
& tile_ni=ini(:), tile_nj=inj(:), &
& tile_data_ibegin=idb(:), tile_data_jbegin=idb(:), &
& tile_data_ni=ini(:) - 2 * idb(:), tile_data_nj=inj(:) - 2 * idb(:))
! Data contains no halo points
idb(:) = 0
CALL iom_set_domain_attr("grid_"//cdgrd//"_inner", ntiles=nijtile, &
& tile_ibegin=ntsi_a(1:nijtile) + idb(:) - 1, tile_jbegin=ntsj_a(1:nijtile) + idb(:) - 1, &
& tile_ibegin=ntsi_a(1:nijtile) - nn_hls - 1, tile_jbegin=ntsj_a(1:nijtile) - nn_hls - 1, &
& tile_ni=ini(:), tile_nj=inj(:), &
& tile_data_ibegin=idb(:), tile_data_jbegin=idb(:), &
& tile_data_ni=ini(:) - 2 * idb(:), tile_data_nj=inj(:) - 2 * idb(:))
......@@ -2440,6 +2468,9 @@ CONTAINS
f_op%timestep = nn_fsbc ; f_of%timestep = 0 ; CALL iom_set_field_attr('SBC' , freq_op=f_op, freq_offset=f_of)
f_op%timestep = nn_fsbc ; f_of%timestep = 0 ; CALL iom_set_field_attr('SBC_scalar' , freq_op=f_op, freq_offset=f_of)
f_op%timestep = nn_fsbc ; f_of%timestep = 0 ; CALL iom_set_field_attr('ABL' , freq_op=f_op, freq_offset=f_of)
IF( ln_diadct ) THEN
f_op%timestep = nn_dctwri ; f_of%timestep = 0 ; CALL iom_set_field_attr('DCT' , freq_op=f_op, freq_offset=f_of)
ENDIF
! output file names (attribut: name)
DO ji = 1, 9
......@@ -2568,7 +2599,7 @@ CONTAINS
INTEGER :: jn, iln
INTEGER :: itrlen
INTEGER :: iyear, imonth, iday, isec
REAL(wp) :: zsec
REAL(dp) :: zsec
LOGICAL :: llexist
TYPE(xios_duration) :: output_freq
!!----------------------------------------------------------------------
......@@ -2668,14 +2699,14 @@ CONTAINS
!!
!! ** Purpose : send back the date corresponding to the given julian day
!!----------------------------------------------------------------------
REAL(wp), INTENT(in ) :: pjday ! julian day
REAL(dp), INTENT(in ) :: pjday ! julian day
LOGICAL , INTENT(in ), OPTIONAL :: ld24 ! true to force 24:00 instead of 00:00
LOGICAL , INTENT(in ), OPTIONAL :: ldfull ! true to get the compleate date: yyyymmdd_hh:mm:ss
!
CHARACTER(LEN=20) :: iom_sdate
CHARACTER(LEN=50) :: clfmt ! format used to write the date
INTEGER :: iyear, imonth, iday, ihour, iminute, isec
REAL(wp) :: zsec
REAL(dp) :: zsec
LOGICAL :: ll24, llfull
!!----------------------------------------------------------------------
!
......
src/OCE/IOM/iom_def.F90
View file @ 672ff863
......@@ -25,7 +25,7 @@ MODULE iom_def
INTEGER, PARAMETER, PUBLIC :: jp_i2 = 203 !: write INTEGER(2)
INTEGER, PARAMETER, PUBLIC :: jp_i1 = 204 !: write INTEGER(1)
INTEGER, PARAMETER, PUBLIC :: jpmax_files = 100 !: maximum number of simultaneously opened file
INTEGER, PARAMETER, PUBLIC :: jpmax_files = 200 !: maximum number of simultaneously opened file
INTEGER, PARAMETER, PUBLIC :: jpmax_vars = 1200 !: maximum number of variables in one file
INTEGER, PARAMETER, PUBLIC :: jpmax_dims = 4 !: maximum number of dimensions for one variable
INTEGER, PARAMETER, PUBLIC :: jpmax_digits = 9 !: maximum number of digits for the cpu number in the file name
......
src/OCE/IOM/prtctl.F90
View file @ 672ff863
......@@ -38,11 +38,11 @@ CONTAINS
SUBROUTINE prt_ctl (tab2d_1, tab3d_1, tab4d_1, tab2d_2, tab3d_2, mask1, mask2, &
& clinfo, clinfo1, clinfo2, clinfo3, kdim )
!!
REAL(wp), DIMENSION(:,:) , INTENT(in), OPTIONAL :: tab2d_1
REAL(wp), DIMENSION(:,:,:) , INTENT(in), OPTIONAL :: tab3d_1
REAL(wp), DIMENSION(:,:,:,:), INTENT(in), OPTIONAL :: tab4d_1
REAL(wp), DIMENSION(:,:) , INTENT(in), OPTIONAL :: tab2d_2
REAL(wp), DIMENSION(:,:,:) , INTENT(in), OPTIONAL :: tab3d_2
REAL(dp), DIMENSION(:,:) , INTENT(in), OPTIONAL :: tab2d_1
REAL(dp), DIMENSION(:,:,:) , INTENT(in), OPTIONAL :: tab3d_1
REAL(dp), DIMENSION(:,:,:,:), INTENT(in), OPTIONAL :: tab4d_1
REAL(dp), DIMENSION(:,:) , INTENT(in), OPTIONAL :: tab2d_2
REAL(dp), DIMENSION(:,:,:) , INTENT(in), OPTIONAL :: tab3d_2
REAL(wp), DIMENSION(:,:,:) , INTENT(in), OPTIONAL :: mask1
REAL(wp), DIMENSION(:,:,:) , INTENT(in), OPTIONAL :: mask2
CHARACTER(len=*), DIMENSION(:) , INTENT(in), OPTIONAL :: clinfo ! information about the tab3d array
......@@ -53,29 +53,29 @@ CONTAINS
!
IF( PRESENT(tab2d_2) ) THEN
CALL prt_ctl_t(ktab2d_1 = is_tile(tab2d_1), ktab3d_1 = 0, ktab4d_1 = 0, ktab2d_2 = is_tile(tab2d_2), ktab3d_2 = 0, &
& tab2d_1 = REAL(tab2d_1, 2*wp), tab2d_2 = REAL(tab2d_2, 2*wp), &
& tab2d_1 = REAL(tab2d_1, dp), tab2d_2 = REAL(tab2d_2, dp), &
& mask1 = mask1, mask2 = mask2, &
& clinfo = clinfo, clinfo1 = clinfo1, clinfo2 = clinfo2, clinfo3 = clinfo3 )
ELSEIF( PRESENT(tab3d_2) ) THEN
CALL prt_ctl_t(ktab2d_1 = 0, ktab3d_1 = is_tile(tab3d_1), ktab4d_1 = 0, ktab2d_2 = 0, ktab3d_2 = is_tile(tab3d_2), &
& tab3d_1 = REAL(tab3d_1, 2*wp), tab3d_2 = REAL(tab3d_2, 2*wp), &
& tab3d_1 = REAL(tab3d_1, dp), tab3d_2 = REAL(tab3d_2, dp), &
& mask1 = mask1, mask2 = mask2, &
& clinfo = clinfo, clinfo1 = clinfo1, clinfo2 = clinfo2, clinfo3 = clinfo3, kdim = kdim )
ELSEIF( PRESENT(tab2d_1) ) THEN
CALL prt_ctl_t(ktab2d_1 = is_tile(tab2d_1), ktab3d_1 = 0, ktab4d_1 = 0, ktab2d_2 = 0, ktab3d_2 = 0, &
& tab2d_1 = REAL(tab2d_1,2*wp), &
& tab2d_1 = REAL(tab2d_1,dp), &
& mask1 = mask1, &
& clinfo = clinfo, clinfo1 = clinfo1, clinfo3 = clinfo3 )
ELSEIF( PRESENT(tab3d_1) ) THEN
CALL prt_ctl_t(ktab2d_1 = 0, ktab3d_1 = is_tile(tab3d_1), ktab4d_1 = 0, ktab2d_2 = 0, ktab3d_2 = 0, &
& tab3d_1 = REAL(tab3d_1, 2*wp), &
& tab3d_1 = REAL(tab3d_1, dp), &
& mask1 = mask1, &
& clinfo = clinfo, clinfo1 = clinfo1, clinfo3 = clinfo3, kdim = kdim )
ELSEIF( PRESENT(tab4d_1) ) THEN
CALL prt_ctl_t(ktab2d_1 = 0, ktab3d_1 = 0, ktab4d_1 = is_tile(tab4d_1), ktab2d_2 = 0, ktab3d_2 = 0, &
& tab4d_1 = REAL(tab4d_1, 2*wp), &
& tab4d_1 = REAL(tab4d_1, dp), &
& mask1 = mask1, &
& clinfo = clinfo, clinfo1 = clinfo1, clinfo3 = clinfo3, kdim = kdim )
& clinfo = clinfo, clinfo1 = clinfo1, clinfo3 = clinfo3, kdim = kdim )
ENDIF
END SUBROUTINE prt_ctl
......@@ -119,11 +119,11 @@ CONTAINS
!! clinfo3 : additional information
!!----------------------------------------------------------------------
INTEGER , INTENT(in) :: ktab2d_1, ktab3d_1, ktab4d_1, ktab2d_2, ktab3d_2
REAL(2*wp), DIMENSION(A2D_T(ktab2d_1)) , INTENT(in), OPTIONAL :: tab2d_1
REAL(2*wp), DIMENSION(A2D_T(ktab3d_1),:) , INTENT(in), OPTIONAL :: tab3d_1
REAL(2*wp), DIMENSION(A2D_T(ktab4d_1),:,:), INTENT(in), OPTIONAL :: tab4d_1
REAL(2*wp), DIMENSION(A2D_T(ktab2d_2)) , INTENT(in), OPTIONAL :: tab2d_2
REAL(2*wp), DIMENSION(A2D_T(ktab3d_2),:) , INTENT(in), OPTIONAL :: tab3d_2
REAL(dp), DIMENSION(A2D_T(ktab2d_1)) , INTENT(in), OPTIONAL :: tab2d_1
REAL(dp), DIMENSION(A2D_T(ktab3d_1),:) , INTENT(in), OPTIONAL :: tab3d_1
REAL(dp), DIMENSION(A2D_T(ktab4d_1),:,:), INTENT(in), OPTIONAL :: tab4d_1
REAL(dp), DIMENSION(A2D_T(ktab2d_2)) , INTENT(in), OPTIONAL :: tab2d_2
REAL(dp), DIMENSION(A2D_T(ktab3d_2),:) , INTENT(in), OPTIONAL :: tab3d_2
REAL(wp), DIMENSION(:,:,:) , INTENT(in), OPTIONAL :: mask1
REAL(wp), DIMENSION(:,:,:) , INTENT(in), OPTIONAL :: mask2
CHARACTER(len=*), DIMENSION(:) , INTENT(in), OPTIONAL :: clinfo ! information about the tab3d array
......@@ -137,7 +137,7 @@ CONTAINS
INTEGER :: jn, jl, kdir
INTEGER :: iis, iie, jjs, jje
INTEGER :: itra, inum
REAL(2*wp) :: zsum1, zsum2, zvctl1, zvctl2
REAL(wp) :: zsum1, zsum2, zvctl1, zvctl2
!!----------------------------------------------------------------------
!
! Arrays, scalars initialization
......
src/OCE/IOM/restart.F90
View file @ 672ff863
......@@ -39,6 +39,7 @@ MODULE restart
!
USE in_out_manager ! I/O manager
USE iom ! I/O module
USE ioipsl, ONLY : ju2ymds ! for calendar
USE lib_mpp ! distribued memory computing library
IMPLICIT NONE
......@@ -72,6 +73,9 @@ CONTAINS
!!----------------------------------------------------------------------
INTEGER, INTENT(in) :: kt ! ocean time-step
!!
INTEGER :: iyear, imonth, iday
REAL (wp) :: zsec
REAL (wp) :: zfjulday
CHARACTER(LEN=20) :: clkt ! ocean time-step deine as a character
CHARACTER(LEN=50) :: clname ! ocean output restart file name
CHARACTER(lc) :: clpath ! full path to ocean output restart file
......@@ -103,8 +107,16 @@ CONTAINS
IF( kt == nitrst - 1 .OR. nn_stock == 1 .OR. ( kt == nitend .AND. .NOT. lrst_oce ) ) THEN
IF( nitrst <= nitend .AND. nitrst > 0 ) THEN
! beware of the format used to write kt (default is i8.8, that should be large enough...)
IF( nitrst > 999999999 ) THEN ; WRITE(clkt, * ) nitrst
ELSE ; WRITE(clkt, '(i8.8)') nitrst
IF ( ln_rstdate ) THEN
zfjulday = fjulday + rdt / rday
IF( ABS(zfjulday - REAL(NINT(zfjulday),wp)) < 0.1 / rday ) zfjulday = REAL(NINT(zfjulday),wp) ! avoid truncation error
CALL ju2ymds( zfjulday, iyear, imonth, iday, zsec )
WRITE(clkt, '(i4.4,2i2.2)') iyear, imonth, iday
ELSE
IF( nitrst > 999999999 ) THEN ; WRITE(clkt, * ) nitrst
ELSE ; WRITE(clkt, '(i8.8)') nitrst
ENDIF
ENDIF
! create the file
clname = TRIM(cexper)//"_"//TRIM(ADJUSTL(clkt))//"_"//TRIM(cn_ocerst_out)
......@@ -183,6 +195,9 @@ CONTAINS
#endif
ENDIF
CALL iom_rstput( kt, nitrst, numrow, 'neos' , REAL(neos)) ! equation of state
!CALL iom_rstput( kt, nitrst, numrow, 'neos' , neos , ktype = jp_i1, ldxios = lwxios) ! equation of state
IF( ln_diurnal ) CALL iom_rstput( kt, nitrst, numrow, 'Dsst', x_dsst )
IF( kt == nitrst ) THEN
IF( .NOT.lwxios ) THEN
......@@ -267,6 +282,8 @@ CONTAINS
!!----------------------------------------------------------------------
INTEGER , INTENT(in) :: Kbb, Kmm ! ocean time level indices
INTEGER :: jk
INTEGER :: id1
REAL(wp) :: zeos
REAL(wp), DIMENSION(jpi, jpj, jpk) :: w3d
REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: zgdept ! 3D workspace for QCO
!!----------------------------------------------------------------------
......@@ -294,11 +311,23 @@ CONTAINS
#else
! !* Read Kmm fields (MLF only)
IF(lwp) WRITE(numout,*) ' Kmm u, v and T-S fields read in the restart file'
CALL iom_get( numror, jpdom_auto, 'un', uu(:,:,: ,Kmm), cd_type = 'U', psgn = -1._wp )
CALL iom_get( numror, jpdom_auto, 'vn', vv(:,:,: ,Kmm), cd_type = 'V', psgn = -1._wp )
CALL iom_get( numror, jpdom_auto, 'un', uu(:,:,: ,Kmm), cd_type = 'U', psgn = -1._dp )
CALL iom_get( numror, jpdom_auto, 'vn', vv(:,:,: ,Kmm), cd_type = 'V', psgn = -1._dp )
CALL iom_get( numror, jpdom_auto, 'tn', ts(:,:,:,jp_tem,Kmm) )
CALL iom_get( numror, jpdom_auto, 'sn', ts(:,:,:,jp_sal,Kmm) )
!
IF ( ln_rst_eos ) THEN
! Check equation of state used is consistent with the restart
IF( iom_varid( numror, 'neos') == -1 ) THEN
CALL ctl_stop( 'restart, rst_read: variable neos not found. STOP check that the equations of state in the restart file and in the namelist nameos are consistent and use ln_rst_eos=F')
ELSE
CALL iom_get( numror, 'neos' , zeos )
IF ( INT(zeos) /= neos ) CALL ctl_stop( 'restart, rst_read: equation of state used in restart file differs from namelist nameos')
ENDIF
ENDIF
IF( l_1st_euler ) THEN !* Euler restart (MLF only)
IF(lwp) WRITE(numout,*) ' Kbb u, v and T-S fields set to Kmm values'
uu(:,:,: ,Kbb) = uu(:,:,: ,Kmm) ! all before fields set to now values
......@@ -307,8 +336,8 @@ CONTAINS
!
ELSE !* Leap frog restart (MLF only)
IF(lwp) WRITE(numout,*) ' Kbb u, v and T-S fields read in the restart file'
CALL iom_get( numror, jpdom_auto, 'ub', uu(:,:,: ,Kbb), cd_type = 'U', psgn = -1._wp )
CALL iom_get( numror, jpdom_auto, 'vb', vv(:,:,: ,Kbb), cd_type = 'V', psgn = -1._wp )
CALL iom_get( numror, jpdom_auto, 'ub', uu(:,:,: ,Kbb), cd_type = 'U', psgn = -1._dp )
CALL iom_get( numror, jpdom_auto, 'vb', vv(:,:,: ,Kbb), cd_type = 'V', psgn = -1._dp )
CALL iom_get( numror, jpdom_auto, 'tb', ts(:,:,:,jp_tem,Kbb) )
CALL iom_get( numror, jpdom_auto, 'sb', ts(:,:,:,jp_sal,Kbb) )
ENDIF
......
src/OCE/ISF/isf_oce.F90
View file @ 672ff863
......@@ -12,7 +12,7 @@ MODULE isf_oce
!!----------------------------------------------------------------------
!! isf : define and allocate ice shelf variables
!!----------------------------------------------------------------------
USE par_kind
USE par_oce , ONLY: jpi, jpj, jpk
USE in_out_manager, ONLY: wp, jpts ! I/O manager
USE lib_mpp , ONLY: ctl_stop, mpp_sum ! MPP library
......
src/OCE/ISF/isfcavmlt.F90
View file @ 672ff863
......@@ -32,6 +32,7 @@ MODULE isfcavmlt
PUBLIC isfcav_mlt
!! * Substitutions
# include "single_precision_substitute.h90"
# include "do_loop_substitute.h90"
!!----------------------------------------------------------------------
!! NEMO/OCE 4.0 , NEMO Consortium (2018)
......@@ -112,7 +113,7 @@ CONTAINS
!!--------------------------------------------------------------------
!
! Compute freezing temperature
CALL eos_fzp( pstbl(:,:), ztfrz(:,:), risfdep(:,:) )
CALL eos_fzp( pstbl(:,:), ztfrz(:,:), CASTDP(risfdep(:,:)) )
!
! read input file of fwf (from isf to oce; ie melt)
CALL fld_read ( kt, 1, sf_isfcav_fwf )
......@@ -157,7 +158,7 @@ CONTAINS
!!--------------------------------------------------------------------
!
! Calculate freezing temperature
CALL eos_fzp( pstbl(:,:), ztfrz(:,:), risfdep(:,:) )
CALL eos_fzp( pstbl(:,:), ztfrz(:,:), CASTDP(risfdep(:,:)) )
!
! thermal driving
zthd (:,:) = ( pttbl(:,:) - ztfrz(:,:) ) * mskisf_cav(:,:)
......@@ -281,7 +282,7 @@ CONTAINS
!!--------------------------------------------------------------------
!
! Calculate freezing temperature
CALL eos_fzp( pstbl(:,:), ztfrz(:,:), risfdep(:,:) )
CALL eos_fzp( pstbl(:,:), ztfrz(:,:), CASTDP(risfdep(:,:)) )
!
! read input file of fwf from isf to oce
CALL fld_read ( kt, 1, sf_isfcav_fwf )
......
src/OCE/ISF/isfcpl.F90
View file @ 672ff863
......@@ -48,6 +48,7 @@ MODULE isfcpl
!
!! * Substitutions
# include "do_loop_substitute.h90"
# include "single_precision_substitute.h90"
# include "domzgr_substitute.h90"
!!----------------------------------------------------------------------
!! NEMO/OCE 4.0 , NEMO Consortium (2018)
......@@ -200,7 +201,8 @@ CONTAINS
zssmask_b(ji,jj) = 1._wp
ENDIF
END_2D
CALL lbc_lnk( 'isfcpl', ssh(:,:,Kmm), 'T', 1.0_wp, zssmask_b(:,:), 'T', 1.0_wp )
CALL lbc_lnk( 'isfcpl', ssh(:,:,Kmm), 'T', 1.0_dp)
CALL lbc_lnk( 'isfcpl', zssmask_b(:,:), 'T', 1.0_wp )
!
zssh(:,:) = ssh(:,:,Kmm)
zssmask0(:,:) = zssmask_b(:,:)
......@@ -213,7 +215,7 @@ CONTAINS
!
ssh(:,:,Kbb) = ssh(:,:,Kmm)
!
IF ( ln_isfdebug ) CALL debug('isfcpl_ssh: sshn',ssh(:,:,Kmm))
IF ( ln_isfdebug ) CALL debug('isfcpl_ssh: sshn',CASTSP(ssh(:,:,Kmm)))
!
! recompute the vertical scale factor, depth and water thickness
IF(lwp) write(numout,*) 'isfcpl_ssh : recompute scale factor from ssh (new wet cell,Kmm)'
......@@ -358,7 +360,8 @@ CONTAINS
END_2D
END DO
!
CALL lbc_lnk( 'isfcpl', ts(:,:,:,jp_tem,Kmm), 'T', 1.0_wp, ts(:,:,:,jp_sal,Kmm), 'T', 1.0_wp, ztmask1, 'T', 1.0_wp)
CALL lbc_lnk( 'isfcpl', ts(:,:,:,jp_tem,Kmm), 'T', 1.0_dp, ts(:,:,:,jp_sal,Kmm), 'T', 1.0_dp)
CALL lbc_lnk( 'isfcpl', ztmask1, 'T', 1.0_wp)
!
! update temperature and salinity and mask
zts0(:,:,:,:) = ts(:,:,:,:,Kmm)
......
src/OCE/ISF/isfdiags.F90
View file @ 672ff863
......@@ -72,7 +72,7 @@ CONTAINS
IF ( iom_use( TRIM(cvarqlat3d) ) ) CALL isf_diags_2dto3d( Kmm, ktop, kbot, phtbl, pfrac, TRIM(cvarqlat3d) , pqoce(:,:))
IF ( iom_use( TRIM(cvarqhc3d) ) ) CALL isf_diags_2dto3d( Kmm, ktop, kbot, phtbl, pfrac, TRIM(cvarqhc3d) , pqhc (:,:))
!
END SUBROUTINE
END SUBROUTINE isf_diags_flx
SUBROUTINE isf_diags_2dto3d(Kmm, ktop, kbot, phtbl, pfrac, cdvar, pvar2d)
!!---------------------------------------------------------------------
......