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
tests/ISOMIP+/MY_SRC/istate.F90
View file @ 4389904b
......@@ -32,6 +32,7 @@ MODULE istate
USE in_out_manager ! I/O manager
USE iom ! I/O library
USE lib_mpp ! MPP library
USE lbclnk ! lateal boundary condition / mpp exchanges
USE restart ! restart
#if defined key_agrif
......@@ -49,7 +50,7 @@ MODULE istate
# include "domzgr_substitute.h90"
!!----------------------------------------------------------------------
!! NEMO/OCE 4.0 , NEMO Consortium (2018)
!! $Id: istate.F90 11423 2019-08-08 14:02:49Z mathiot $
!! $Id: istate.F90 15581 2021-12-07 13:08:22Z techene $
!! Software governed by the CeCILL license (see ./LICENSE)
!!----------------------------------------------------------------------
CONTAINS
......@@ -59,6 +60,8 @@ CONTAINS
!! *** ROUTINE istate_init ***
!!
!! ** Purpose : Initialization of the dynamics and tracer fields.
!!
!! ** Method :
!!----------------------------------------------------------------------
INTEGER, INTENT( in ) :: Kbb, Kmm, Kaa ! ocean time level indices
!
......@@ -86,7 +89,7 @@ CONTAINS
#endif
#if defined key_agrif
IF ( (.NOT.Agrif_root()).AND.ln_init_chfrpar ) THEN
IF ( .NOT.Agrif_root() .AND. ln_init_chfrpar ) THEN
numror = 0 ! define numror = 0 -> no restart file to read
ln_1st_euler = .true. ! Set time-step indicator at nit000 (euler forward)
CALL day_init
......@@ -125,37 +128,50 @@ CONTAINS
DO jk = 1, jpk
zgdept(:,:,jk) = gdept(:,:,jk,Kbb)
END DO
CALL usr_def_istate( zgdept, tmask, ts(:,:,:,:,Kbb), uu(:,:,:,Kbb), vv(:,:,:,Kbb) )
CALL usr_def_istate( zgdept, tmask, ts(:,:,:,:,Kbb), uu(:,:,:,Kbb), vv(:,:,:,Kbb) )
! make sure that periodicities are properly applied
CALL lbc_lnk( 'istate', ts(:,:,:,jp_tem,Kbb), 'T', 1._wp, ts(:,:,:,jp_sal,Kbb), 'T', 1._wp, &
& uu(:,:,:, Kbb), 'U', -1._wp, vv(:,:,:, Kbb), 'V', -1._wp )
ENDIF
ts (:,:,:,:,Kmm) = ts (:,:,:,:,Kbb) ! set now values from to before ones
uu (:,:,:,Kmm) = uu (:,:,:,Kbb)
vv (:,:,:,Kmm) = vv (:,:,:,Kbb)
ENDIF
#if defined key_agrif
ENDIF
#endif
!
! Initialize "now" and "before" barotropic velocities:
! Do it whatever the free surface method, these arrays being eventually used
! Initialize "now" barotropic velocities:
! Do it whatever the free surface method, these arrays being used eventually
!
!!gm the use of umask & vmask is not necessary below as uu(:,:,:,Kmm), vv(:,:,:,Kmm), uu(:,:,:,Kbb), vv(:,:,:,Kbb) are always masked
#if ! defined key_RK3
uu_b(:,:,Kmm) = 0._wp ; vv_b(:,:,Kmm) = 0._wp
uu_b(:,:,Kbb) = 0._wp ; vv_b(:,:,Kbb) = 0._wp
!
!!gm the use of umsak & vmask is not necessary below as uu(:,:,:,Kmm), vv(:,:,:,Kmm), uu(:,:,:,Kbb), vv(:,:,:,Kbb) are always masked
DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1 )
uu_b(ji,jj,Kmm) = uu_b(ji,jj,Kmm) + e3u(ji,jj,jk,Kmm) * uu(ji,jj,jk,Kmm) * umask(ji,jj,jk)
vv_b(ji,jj,Kmm) = vv_b(ji,jj,Kmm) + e3v(ji,jj,jk,Kmm) * vv(ji,jj,jk,Kmm) * vmask(ji,jj,jk)
!
uu_b(ji,jj,Kbb) = uu_b(ji,jj,Kbb) + e3u(ji,jj,jk,Kbb) * uu(ji,jj,jk,Kbb) * umask(ji,jj,jk)
vv_b(ji,jj,Kbb) = vv_b(ji,jj,Kbb) + e3v(ji,jj,jk,Kbb) * vv(ji,jj,jk,Kbb) * vmask(ji,jj,jk)
END_3D
!
uu_b(:,:,Kmm) = uu_b(:,:,Kmm) * r1_hu(:,:,Kmm)
vv_b(:,:,Kmm) = vv_b(:,:,Kmm) * r1_hv(:,:,Kmm)
#endif
!
uu_b(:,:,Kbb) = uu_b(:,:,Kbb) * r1_hu(:,:,Kbb)
vv_b(:,:,Kbb) = vv_b(:,:,Kbb) * r1_hv(:,:,Kbb)
#if defined key_RK3
IF( .NOT. ln_rstart ) THEN
#endif
! Initialize "before" barotropic velocities. "now" values are always set but
! "before" values may have been read from a restart to ensure restartability.
! In the non-restart or non-RK3 cases they need to be initialised here:
uu_b(:,:,Kbb) = 0._wp ; vv_b(:,:,Kbb) = 0._wp
DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1 )
uu_b(ji,jj,Kbb) = uu_b(ji,jj,Kbb) + e3u(ji,jj,jk,Kbb) * uu(ji,jj,jk,Kbb) * umask(ji,jj,jk)
vv_b(ji,jj,Kbb) = vv_b(ji,jj,Kbb) + e3v(ji,jj,jk,Kbb) * vv(ji,jj,jk,Kbb) * vmask(ji,jj,jk)
END_3D
uu_b(:,:,Kbb) = uu_b(:,:,Kbb) * r1_hu(:,:,Kbb)
vv_b(:,:,Kbb) = vv_b(:,:,Kbb) * r1_hv(:,:,Kbb)
!
#if defined key_RK3
ENDIF
#endif
!
END SUBROUTINE istate_init
......
tests/ISOMIP+/MY_SRC/sbcfwb.F90
View file @ 4389904b
......@@ -35,8 +35,9 @@ MODULE sbcfwb
PUBLIC sbc_fwb ! routine called by step
REAL(wp) :: rn_fwb0 ! initial freshwater adjustment flux [kg/m2/s] (nn_fwb = 2 only)
REAL(wp) :: a_fwb ! annual domain averaged freshwater budget from the
! previous year
REAL(wp) :: a_fwb ! annual domain averaged freshwater budget from the previous year
REAL(wp) :: a_fwb_b ! annual domain averaged freshwater budget from the year before or at initial state
REAL(wp) :: a_fwb_ini ! initial domain averaged freshwater budget
REAL(wp) :: area ! global mean ocean surface (interior domain)
!!----------------------------------------------------------------------
......@@ -128,68 +129,65 @@ CONTAINS
IF( iom_use('vflx_fwb_cea') ) CALL iom_put( 'vflx_fwb_cea', z_fwfprv(1) * tmask(:,:,1) )
ENDIF
!
CASE ( 4 ) !== global mean fwf set to zero (ISOMIP case) ==!
!
IF( MOD( kt-1, kn_fsbc ) == 0 ) THEN
z_fwf = glob_sum( 'sbcfwb', e1e2t(:,:) * ( emp(:,:) - rnf(:,:) - fwfisf_cav(:,:) - fwfisf_par(:,:) - snwice_fmass(:,:) ) )
!
! correction for ice sheet coupling testing (ie remove the excess through the surface)
! test impact on the melt as conservation correction made in depth
! test conservation level as sbcfwb is conserving
! avoid the model to blow up for large ssh drop (isomip OCEAN3 with melt switch off and uniform T/S)
IF (ln_isfcpl .AND. ln_isfcpl_cons) THEN
z_fwf = z_fwf + glob_sum( 'sbcfwb', e1e2t(:,:) * risfcpl_cons_ssh(:,:) * rho0 )
END IF
!
z_fwf = z_fwf / area
zcoef = z_fwf * rcp
emp(:,:) = emp(:,:) - z_fwf * tmask(:,:,1) ! (Eq. 34 AD2015)
qns(:,:) = qns(:,:) + zcoef * sst_m(:,:) * tmask(:,:,1) ! (Eq. 35 AD2015) ! use sst_m to avoid generation of any bouyancy fluxes
sfx(:,:) = sfx(:,:) + z_fwf * sss_m(:,:) * tmask(:,:,1) ! (Eq. 36 AD2015) ! use sss_m to avoid generation of any bouyancy fluxes
ENDIF
!
CASE ( 2 ) !== fw adjustment based on fw budget at the end of the previous year ==!
!
IF( kt == nit000 ) THEN ! initialisation
! ! set the fw adjustment (a_fwb)
IF ( ln_rstart .AND. iom_varid( numror, 'a_fwb', ldstop = .FALSE. ) > 0 ) THEN ! as read from restart file
IF(lwp) WRITE(numout,*) 'sbc_fwb : reading FW-budget adjustment from restart file'
CALL iom_get( numror, 'a_fwb', a_fwb )
ELSE ! as specified in namelist
a_fwb = rn_fwb0
! simulation is supposed to start 1st of January
IF( kt == nit000 ) THEN ! initialisation
! ! set the fw adjustment (a_fwb)
IF ( ln_rstart .AND. iom_varid( numror, 'a_fwb_b', ldstop = .FALSE. ) > 0 & ! as read from restart file
& .AND. iom_varid( numror, 'a_fwb', ldstop = .FALSE. ) > 0 ) THEN
IF(lwp) WRITE(numout,*) 'sbc_fwb : reading freshwater-budget from restart file'
CALL iom_get( numror, 'a_fwb_b', a_fwb_b )
CALL iom_get( numror, 'a_fwb' , a_fwb )
!
a_fwb_ini = a_fwb_b
ELSE ! as specified in namelist
IF(lwp) WRITE(numout,*) 'sbc_fwb : setting freshwater-budget from namelist rn_fwb0'
a_fwb = rn_fwb0
a_fwb_b = 0._wp ! used only the first year then it is replaced by a_fwb_ini
!
a_fwb_ini = glob_sum( 'sbcfwb', e1e2t(:,:) * ( ssh(:,:,Kmm) + snwice_mass(:,:) * r1_rho0 ) ) &
& * rho0 / ( area * rday * REAL(nyear_len(1), wp) )
END IF
!
IF(lwp)WRITE(numout,*)
IF(lwp)WRITE(numout,*)'sbc_fwb : initial freshwater-budget adjustment = ', a_fwb, 'kg/m2/s'
IF(lwp) WRITE(numout,*)
IF(lwp) WRITE(numout,*)'sbc_fwb : freshwater-budget at the end of previous year = ', a_fwb , 'kg/m2/s'
IF(lwp) WRITE(numout,*)' freshwater-budget at initial state = ', a_fwb_ini, 'kg/m2/s'
!
ELSE
! at the end of year n:
ikty = nyear_len(1) * 86400 / NINT(rn_Dt)
IF( MOD( kt, ikty ) == 0 ) THEN ! Update a_fwb at the last time step of a year
! It should be the first time step of a year MOD(kt-1,ikty) but then the restart would be wrong
! Hence, we make a small error here but the code is restartable
a_fwb_b = a_fwb_ini
! mean sea level taking into account ice+snow
a_fwb = glob_sum( 'sbcfwb', e1e2t(:,:) * ( ssh(:,:,Kmm) + snwice_mass(:,:) * r1_rho0 ) )
a_fwb = a_fwb * rho0 / ( area * rday * REAL(nyear_len(1), wp) ) ! convert in kg/m2/s
ENDIF
!
ENDIF
! ! Update a_fwb if new year start
ikty = 365 * 86400 / rn_Dt !!bug use of 365 days leap year or 360d year !!!!!!!
IF( MOD( kt, ikty ) == 0 ) THEN
! mean sea level taking into account the ice+snow
! sum over the global domain
a_fwb = glob_sum( 'sbcfwb', e1e2t(:,:) * ( ssh(:,:,Kmm) + snwice_mass(:,:) * r1_rho0 ) )
a_fwb = a_fwb * 1.e+3 / ( area * rday * 365. ) ! convert in Kg/m3/s = mm/s
!!gm ! !!bug 365d year
ENDIF
!
IF( MOD( kt-1, kn_fsbc ) == 0 ) THEN ! correct the freshwater fluxes
zcoef = a_fwb * rcp
emp(:,:) = emp(:,:) + a_fwb * tmask(:,:,1)
qns(:,:) = qns(:,:) - zcoef * sst_m(:,:) * tmask(:,:,1) ! account for change to the heat budget due to fw correction
!
IF( MOD( kt-1, kn_fsbc ) == 0 ) THEN ! correct the freshwater fluxes using previous year budget minus initial state
zcoef = ( a_fwb - a_fwb_b )
emp(:,:) = emp(:,:) + zcoef * tmask(:,:,1)
qns(:,:) = qns(:,:) - zcoef * rcp * sst_m(:,:) * tmask(:,:,1) ! account for change to the heat budget due to fw correction
! outputs
IF( iom_use('hflx_fwb_cea') ) CALL iom_put( 'hflx_fwb_cea', -zcoef * sst_m(:,:) * tmask(:,:,1) )
IF( iom_use('vflx_fwb_cea') ) CALL iom_put( 'vflx_fwb_cea', -a_fwb * tmask(:,:,1) )
IF( iom_use('hflx_fwb_cea') ) CALL iom_put( 'hflx_fwb_cea', -zcoef * rcp * sst_m(:,:) * tmask(:,:,1) )
IF( iom_use('vflx_fwb_cea') ) CALL iom_put( 'vflx_fwb_cea', -zcoef * tmask(:,:,1) )
ENDIF
! Output restart information
IF( lrst_oce ) THEN
IF(lwp) WRITE(numout,*)
IF(lwp) WRITE(numout,*) 'sbc_fwb : writing FW-budget adjustment to ocean restart file at it = ', kt
IF(lwp) WRITE(numout,*) '~~~~'
CALL iom_rstput( kt, nitrst, numrow, 'a_fwb', a_fwb )
CALL iom_rstput( kt, nitrst, numrow, 'a_fwb_b', a_fwb_b )
CALL iom_rstput( kt, nitrst, numrow, 'a_fwb', a_fwb )
END IF
!
IF( kt == nitend .AND. lwp ) WRITE(numout,*) 'sbc_fwb : final freshwater-budget adjustment = ', a_fwb, 'kg/m2/s'
IF( kt == nitend .AND. lwp ) THEN
WRITE(numout,*) 'sbc_fwb : freshwater-budget at the end of simulation (year now) = ', a_fwb , 'kg/m2/s'
WRITE(numout,*) ' freshwater-budget at initial state = ', a_fwb_b, 'kg/m2/s'
ENDIF
!
CASE ( 3 ) !== global fwf set to zero and spread out over erp area ==!
!
......@@ -248,6 +246,26 @@ CONTAINS
ENDIF
DEALLOCATE( ztmsk_neg , ztmsk_pos , ztmsk_tospread , z_wgt , zerp_cor )
!
CASE ( 4 ) !== global mean fwf set to zero (ISOMIP case) ==!
!
IF( MOD( kt-1, kn_fsbc ) == 0 ) THEN
z_fwf = glob_sum( 'sbcfwb', e1e2t(:,:) * ( emp(:,:) - rnf(:,:) - fwfisf_cav(:,:) - fwfisf_par(:,:) - snwice_fmass(:,:) ) )
!
! correction for ice sheet coupling testing (ie remove the excess through the surface)
! test impact on the melt as conservation correction made in depth
! test conservation level as sbcfwb is conserving
! avoid the model to blow up for large ssh drop (isomip OCEAN3 with melt switch off and uniform T/S)
IF (ln_isfcpl .AND. ln_isfcpl_cons) THEN
z_fwf = z_fwf + glob_sum( 'sbcfwb', e1e2t(:,:) * risfcpl_cons_ssh(:,:) * rho0 )
END IF
!
z_fwf = z_fwf / area
zcoef = z_fwf * rcp
emp(:,:) = emp(:,:) - z_fwf * tmask(:,:,1) ! (Eq. 34 AD2015)
qns(:,:) = qns(:,:) + zcoef * sst_m(:,:) * tmask(:,:,1) ! (Eq. 35 AD2015) ! use sst_m to avoid generation of any bouyancy fluxes
sfx(:,:) = sfx(:,:) + z_fwf * sss_m(:,:) * tmask(:,:,1) ! (Eq. 36 AD2015) ! use sss_m to avoid generation of any bouyancy fluxes
ENDIF
!
CASE DEFAULT !== you should never be there ==!
CALL ctl_stop( 'sbc_fwb : wrong nn_fwb value for the FreshWater Budget correction, choose either 1, 2 or 3' )
!
......
tests/ISOMIP+/MY_SRC/tradmp.F90
View file @ 4389904b
......@@ -23,7 +23,6 @@ MODULE tradmp
!!----------------------------------------------------------------------
USE oce ! ocean: variables
USE dom_oce ! ocean: domain variables
USE c1d ! 1D vertical configuration
USE trd_oce ! trends: ocean variables
USE trdtra ! trends manager: tracers
USE zdf_oce ! ocean: vertical physics
......@@ -55,7 +54,7 @@ MODULE tradmp
# include "domzgr_substitute.h90"
!!----------------------------------------------------------------------
!! NEMO/OCE 4.0 , NEMO Consortium (2018)
!! $Id: tradmp.F90 10425 2018-12-19 21:54:16Z smasson $
!! $Id: tradmp.F90 15574 2021-12-03 19:32:50Z techene $
!! Software governed by the CeCILL license (see ./LICENSE)
!!----------------------------------------------------------------------
CONTAINS
......@@ -96,15 +95,19 @@ CONTAINS
!
INTEGER :: ji, jj, jk, jn ! dummy loop indices
REAL(wp), DIMENSION(A2D(nn_hls),jpk,jpts) :: zts_dta
REAL(wp), DIMENSION(jpi,jpj,jpk) :: ze3t
REAL(wp), DIMENSION(:,:,:) , ALLOCATABLE :: zwrk
REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: ztrdts
!!----------------------------------------------------------------------
!
IF( ln_timing ) CALL timing_start('tra_dmp')
!
IF( l_trdtra .OR. iom_use('hflx_dmp_cea') .OR. iom_use('sflx_dmp_cea') ) THEN !* Save ta and sa trends
ALLOCATE( ztrdts(jpi,jpj,jpk,jpts) )
ztrdts(:,:,:,:) = pts(:,:,:,:,Krhs)
ALLOCATE( ztrdts(A2D(nn_hls),jpk,jpts) )
DO jn = 1, jpts
DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpk )
ztrdts(ji,jj,jk,jn) = pts(ji,jj,jk,jn,Krhs)
END_3D
END DO
ENDIF
! !== input T-S data at kt ==!
CALL dta_tsd( kt, 'dmp', zts_dta ) ! read and interpolates T-S data at kt
......@@ -142,16 +145,25 @@ CONTAINS
END SELECT
!
! outputs (clem trunk)
DO jk = 1, jpk
ze3t(:,:,jk) = e3t(:,:,jk,Kmm)
END DO
!
IF( iom_use('hflx_dmp_cea') ) &
& CALL iom_put('hflx_dmp_cea', &
& SUM( ( pts(:,:,:,jp_tem,Krhs) - ztrdts(:,:,:,jp_tem) ) * ze3t(:,:,:), dim=3 ) * rcp * rho0 ) ! W/m2
IF( iom_use('sflx_dmp_cea') ) &
& CALL iom_put('sflx_dmp_cea', &
& SUM( ( pts(:,:,:,jp_sal,Krhs) - ztrdts(:,:,:,jp_sal) ) * ze3t(:,:,:), dim=3 ) * rho0 ) ! g/m2/s
IF( iom_use('hflx_dmp_cea') .OR. iom_use('sflx_dmp_cea') ) THEN
ALLOCATE( zwrk(A2D(nn_hls),jpk) ) ! Needed to handle expressions containing e3t when using key_qco or key_linssh
zwrk(:,:,:) = 0._wp
IF( iom_use('hflx_dmp_cea') ) THEN
DO_3D( 0, 0, 0, 0, 1, jpk )
zwrk(ji,jj,jk) = ( pts(ji,jj,jk,jp_tem,Krhs) - ztrdts(ji,jj,jk,jp_tem) ) * e3t(ji,jj,jk,Kmm)
END_3D
CALL iom_put('hflx_dmp_cea', SUM( zwrk(:,:,:), dim=3 ) * rcp * rho0 ) ! W/m2
ENDIF
IF( iom_use('sflx_dmp_cea') ) THEN
DO_3D( 0, 0, 0, 0, 1, jpk )
zwrk(ji,jj,jk) = ( pts(ji,jj,jk,jp_sal,Krhs) - ztrdts(ji,jj,jk,jp_sal) ) * e3t(ji,jj,jk,Kmm)
END_3D
CALL iom_put('sflx_dmp_cea', SUM( zwrk(:,:,:), dim=3 ) * rho0 ) ! g/m2/s
ENDIF
DEALLOCATE( zwrk )
ENDIF
!
IF( l_trdtra ) THEN ! trend diagnostic
ztrdts(:,:,:,:) = pts(:,:,:,:,Krhs) - ztrdts(:,:,:,:)
......