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
src/OCE/ZDF/zdftke.F90
View file @ 4a9f4b18
......@@ -84,6 +84,7 @@ MODULE zdftke
REAL(wp) :: rn_bshear ! background shear (>0) currently a numerical threshold (do not change it)
INTEGER :: nn_etau ! type of depth penetration of surface tke (=0/1/2/3)
INTEGER :: nn_htau ! type of tke profile of penetration (=0/1)
REAL(wp) :: rn_htau_scaling ! a acaling factor to apply to the penetration of TKE
INTEGER :: nn_bc_surf! surface condition (0/1=Dir/Neum) ! Only applicable for wave coupling
INTEGER :: nn_bc_bot ! surface condition (0/1=Dir/Neum) ! Only applicable for wave coupling
REAL(wp) :: rn_efr ! fraction of TKE surface value which penetrates in the ocean
......@@ -96,7 +97,7 @@ MODULE zdftke
REAL(wp) :: rhftau_add = 1.e-3_wp ! add offset applied to HF part of taum (nn_etau=3)
REAL(wp) :: rhftau_scl = 1.0_wp ! scale factor applied to HF part of taum (nn_etau=3)
REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: htau ! depth of tke penetration (nn_htau)
REAL(wp), ALLOCATABLE, SAVE, PUBLIC, DIMENSION(:,:) :: htau ! depth of tke penetration (nn_htau)
REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: dissl ! now mixing lenght of dissipation
REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: apdlr ! now mixing lenght of dissipation
......@@ -733,6 +734,7 @@ CONTAINS
& rn_emin0, rn_bshear, nn_mxl , ln_mxl0 , &
& rn_mxl0 , nn_mxlice, rn_mxlice, &
& nn_pdl , ln_lc , rn_lc , &
& rn_htau_scaling , &
& nn_etau , nn_htau , rn_efr , nn_eice , &
& nn_bc_surf, nn_bc_bot, ln_mxhsw
!!----------------------------------------------------------------------
......@@ -789,6 +791,7 @@ CONTAINS
ENDIF
WRITE(numout,*) ' test param. to add tke induced by wind nn_etau = ', nn_etau
WRITE(numout,*) ' type of tke penetration profile nn_htau = ', nn_htau
WRITE(numout,*) ' scaling factor for tke penetration depth rn_htau_scaling = ', rn_htau_scaling
WRITE(numout,*) ' fraction of TKE that penetrates rn_efr = ', rn_efr
WRITE(numout,*) ' langmuir & surface wave breaking under ice nn_eice = ', nn_eice
SELECT CASE( nn_eice )
......@@ -819,7 +822,7 @@ CONTAINS
! !* Check of some namelist values
IF( nn_mxl < 0 .OR. nn_mxl > 3 ) CALL ctl_stop( 'bad flag: nn_mxl is 0, 1, 2 or 3' )
IF( nn_pdl < 0 .OR. nn_pdl > 1 ) CALL ctl_stop( 'bad flag: nn_pdl is 0 or 1' )
IF( nn_htau < 0 .OR. nn_htau > 1 ) CALL ctl_stop( 'bad flag: nn_htau is 0 or 1' )
IF( ( nn_htau < 0 .OR. nn_htau > 1 ) .AND. nn_htau .NE. 4 .AND. nn_htau .NE. 5 ) CALL ctl_stop( 'bad flag: nn_htau is 0, 1 , 4 or 5 ' )
IF( nn_etau == 3 .AND. .NOT. ln_cpl ) CALL ctl_stop( 'nn_etau == 3 : HF taum only known in coupled mode' )
!
IF( ln_mxl0 ) THEN
......@@ -831,9 +834,24 @@ CONTAINS
IF( nn_etau /= 0 ) THEN
SELECT CASE( nn_htau ) ! Choice of the depth of penetration
CASE( 0 ) ! constant depth penetration (here 10 meters)
htau(:,:) = 10._wp
htau(:,:) = rn_htau_scaling*10._wp
CASE( 1 ) ! F(latitude) : 0.5m to 30m poleward of 40 degrees
htau(:,:) = MAX( 0.5_wp, MIN( 30._wp, 45._wp* ABS( SIN( rpi/180._wp * gphit(:,:) ) ) ) )
CASE( 4 ) ! F(latitude) : 0.5m to 10m/30m poleward of 13/40 degrees north/south
DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
IF( gphit(ji,jj) <= 0._wp ) THEN
htau(ji,jj) = MAX( 0.5_wp, MIN( 30._wp, 45._wp* rn_htau_scaling*ABS( SIN( rpi/180._wp * gphit(ji,jj) ) ) ) )
ELSE
htau(ji,jj) = MAX( 0.5_wp, MIN( 10._wp, 45._wp* rn_htau_scaling*ABS( SIN( rpi/180._wp * gphit(ji,jj) ) ) ) )
ENDIF
END_2D
CASE( 5 ) ! Variation on case 4 with a steeper ramp further south in Southern Hemisphere
DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
htau(ji,jj) = MAX( 0.5_wp, MIN( 10._wp, 45._wp* rn_htau_scaling*ABS( SIN( rpi/180._wp * gphit(ji,jj) ) ) ) )
IF( gphit(ji,jj) <= -40._wp ) THEN
htau(ji,jj) = htau(ji,jj) + MIN( 20._wp, 135._wp * ABS( SIN( rpi/180._wp * (gphit(ji,jj) + 40.0) ) ) )
ENDIF
END_2D
END SELECT
ENDIF
! !* read or initialize all required files
......
src/OCE/ZDF/zdftmx.F90 0 → 100644
View file @ 4a9f4b18
MODULE zdftmx
!!========================================================================
!! *** MODULE zdftmx ***
!! Ocean physics: vertical tidal mixing coefficient
!!========================================================================
!! History : 1.0 ! 2004-04 (L. Bessieres, G. Madec) Original code
!! - ! 2006-08 (A. Koch-Larrouy) Indonesian strait
!! 3.3 ! 2010-10 (C. Ethe, G. Madec) reorganisation of initialisation phase
!!----------------------------------------------------------------------
!!----------------------------------------------------------------------
!! 'key_zdftmx' Tidal vertical mixing
!!----------------------------------------------------------------------
!! zdf_tmx : global momentum & tracer Kz with tidal induced Kz
!! tmx_itf : Indonesian momentum & tracer Kz with tidal induced Kz
!!----------------------------------------------------------------------
USE oce ! ocean dynamics and tracers variables
USE dom_oce ! ocean space and time domain variables
USE zdf_oce ! ocean vertical physics variables
USE lbclnk ! ocean lateral boundary conditions (or mpp link)
USE eosbn2 ! ocean equation of state
USE phycst ! physical constants
USE prtctl ! Print control
USE in_out_manager ! I/O manager
USE iom ! I/O Manager
USE lib_mpp ! MPP library
USE timing ! Timing
USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)
IMPLICIT NONE
PRIVATE
PUBLIC zdf_tmx ! called in step module
PUBLIC zdf_tmx_init ! called in opa module
PUBLIC zdf_tmx_alloc ! called in nemogcm module
LOGICAL, PUBLIC, PARAMETER :: lk_zdftmx = .TRUE. !: tidal mixing flag
! !!* Namelist namzdf_tmx : tidal mixing *
REAL(wp) :: rn_htmx ! vertical decay scale for turbulence (meters)
REAL(wp) :: rn_n2min ! threshold of the Brunt-Vaisala frequency (s-1)
REAL(wp) :: rn_tfe ! tidal dissipation efficiency (St Laurent et al. 2002)
REAL(wp) :: rn_me ! mixing efficiency (Osborn 1980)
LOGICAL :: ln_tmx_itf ! Indonesian Through Flow (ITF): Koch-Larrouy et al. (2007) parameterization
REAL(wp) :: rn_tfe_itf ! ITF tidal dissipation efficiency (St Laurent et al. 2002)
REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: en_tmx ! energy available for tidal mixing (W/m2)
REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: mask_itf ! mask to use over Indonesian area
REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: az_tmx ! coefficient used to evaluate the tidal induced Kz
!! * Substitutions
# include "do_loop_substitute.h90"
# include "domzgr_substitute.h90"
!!----------------------------------------------------------------------
!! NEMO/OPA 4.0 , NEMO Consortium (2011)
!! $Id: zdftmx.F90 8788 2017-11-22 18:01:02Z davestorkey $
!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
!!----------------------------------------------------------------------
CONTAINS
INTEGER FUNCTION zdf_tmx_alloc()
!!----------------------------------------------------------------------
!! *** FUNCTION zdf_tmx_alloc ***
!!----------------------------------------------------------------------
ALLOCATE(en_tmx(jpi,jpj), mask_itf(jpi,jpj), az_tmx(jpi,jpj,jpk), STAT=zdf_tmx_alloc )
!
IF( lk_mpp ) CALL mpp_sum ('zdftmx', zdf_tmx_alloc )
IF( zdf_tmx_alloc /= 0 ) CALL ctl_warn('zdf_tmx_alloc: failed to allocate arrays')
END FUNCTION zdf_tmx_alloc
SUBROUTINE zdf_tmx( kt, Kmm, p_avm, p_avt, p_avs)
!!----------------------------------------------------------------------
!! *** ROUTINE zdf_tmx ***
!!
!! ** Purpose : add to the vertical mixing coefficients the effect of
!! tidal mixing (Simmons et al 2004).
!!
!! ** Method : - tidal-induced vertical mixing is given by:
!! Kz_tides = az_tmx / max( rn_n2min, N^2 )
!! where az_tmx is a coefficient that specified the 3D space
!! distribution of the faction of tidal energy taht is used
!! for mixing. Its expression is set in zdf_tmx_init routine,
!! following Simmons et al. 2004.
!! NB: a specific bounding procedure is performed on av_tide
!! so that the input tidal energy is actually almost used. The
!! basic maximum value is 60 cm2/s, but values of 300 cm2/s
!! can be reached in area where bottom stratification is too
!! weak.
!!
!! - update av_tide in the Indonesian Through Flow area
!! following Koch-Larrouy et al. (2007) parameterisation
!! (see tmx_itf routine).
!!
!! - update the model vertical eddy viscosity and diffusivity:
!! avt = avt + av_tides
!! avm = avm + av_tides
!!
!! ** Action : avt, avm increased by tidal mixing
!!
!! References : Simmons et al. 2004, Ocean Modelling, 6, 3-4, 245-263.
!! Koch-Larrouy et al. 2007, GRL.
!!----------------------------------------------------------------------
INTEGER, INTENT(in) :: kt , Kmm ! ocean time-step , time level
REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: p_avm ! momentum Kz (w-points)
REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: p_avt, p_avs ! tracer Kz (w-points)
!!
INTEGER :: ji, jj, jk ! dummy loop indices
REAL(wp) :: ztpc ! scalar workspace
REAL(wp), DIMENSION(jpi,jpj) :: zkz
REAL(wp), DIMENSION(jpi,jpj,jpk) :: zav_tide
!!----------------------------------------------------------------------
!
! ! ----------------------- !
! ! Standard tidal mixing ! (compute zav_tide)
! ! ----------------------- !
! !* First estimation (with n2 bound by rn_n2min) bounded by 60 cm2/s
zav_tide(:,:,:) = MIN( 60.e-4, az_tmx(:,:,:) / MAX( rn_n2min, rn2(:,:,:) ) )
zkz(:,:) = 0.e0 !* Associated potential energy consummed over the whole water column
DO jk = 2, jpkm1
zkz(:,:) = zkz(:,:) + e3w(:,:,jk,Kmm) * MAX( 0.e0, rn2(:,:,jk) ) * rho0 * zav_tide(:,:,jk) * wmask(:,:,jk)
END DO
DO jj = 1, jpj !* Here zkz should be equal to en_tmx ==> multiply by en_tmx/zkz to recover en_tmx
DO ji = 1, jpi
IF( zkz(ji,jj) /= 0.e0 ) zkz(ji,jj) = en_tmx(ji,jj) / zkz(ji,jj)
END DO
END DO
DO jk = 2, jpkm1 !* Mutiply by zkz to recover en_tmx, BUT bound by 30/6 ==> zav_tide bound by 300 cm2/s
DO jj = 1, jpj !* Here zkz should be equal to en_tmx ==> multiply by en_tmx/zkz to recover en_tmx
DO ji = 1, jpi
zav_tide(ji,jj,jk) = zav_tide(ji,jj,jk) * MIN( zkz(ji,jj), 30./6. ) * wmask(ji,jj,jk) !kz max = 300 cm2/s
END DO
END DO
END DO
IF( kt == nit000 ) THEN !* check at first time-step: diagnose the energy consumed by zav_tide
ztpc = 0.e0
DO jk= 1, jpk
DO jj= 1, jpj
DO ji= 1, jpi
ztpc = ztpc + e3w(ji,jj,jk,Kmm) * e1t(ji,jj) * e2t(ji,jj) &
& * MAX( 0.e0, rn2(ji,jj,jk) ) * zav_tide(ji,jj,jk) * tmask(ji,jj,jk) * tmask_i(ji,jj)
END DO
END DO
END DO
ztpc= rho0 / ( rn_tfe * rn_me ) * ztpc
IF(lwp) WRITE(numout,*)
IF(lwp) WRITE(numout,*) ' N Total power consumption by av_tide : ztpc = ', ztpc * 1.e-12 ,'TW'
ENDIF
! ! ----------------------- !
! ! ITF tidal mixing ! (update zav_tide)
! ! ----------------------- !
IF( ln_tmx_itf ) CALL tmx_itf( kt, zav_tide , Kmm)
! ! ----------------------- !
! ! Update mixing coefs !
! ! ----------------------- !
DO jk = 2, jpkm1 !* update momentum & tracer diffusivity with tidal mixing
DO jj = 1, jpj !* Here zkz should be equal to en_tmx ==> multiply by en_tmx/zkz to recover en_tmx
DO ji = 1, jpi
p_avt(ji,jj,jk) = p_avt(ji,jj,jk) + zav_tide(ji,jj,jk) * wmask(ji,jj,jk)
p_avs(ji,jj,jk) = p_avs(ji,jj,jk) + zav_tide(ji,jj,jk) * wmask(ji,jj,jk)
p_avm(ji,jj,jk) = p_avm(ji,jj,jk) + zav_tide(ji,jj,jk) * wmask(ji,jj,jk)
END DO
END DO
END DO
! !* output tidal mixing coefficient
CALL iom_put( "av_tmx", zav_tide )
IF(sn_cfctl%l_prtctl) CALL prt_ctl(tab3d_1=zav_tide , clinfo1=' tmx - av_tide: ', tab3d_2=p_avt, clinfo2=' p_avt: ', kdim=jpk)
!
END SUBROUTINE zdf_tmx
SUBROUTINE tmx_itf( kt, pav , Kmm)
!!----------------------------------------------------------------------
!! *** ROUTINE tmx_itf ***
!!
!! ** Purpose : modify the vertical eddy diffusivity coefficients
!! (pav) in the Indonesian Through Flow area (ITF).
!!
!! ** Method : - Following Koch-Larrouy et al. (2007), in the ITF defined
!! by msk_itf (read in a file, see tmx_init), the tidal
!! mixing coefficient is computed with :
!! * q=1 (i.e. all the tidal energy remains trapped in
!! the area and thus is used for mixing)
!! * the vertical distribution of the tifal energy is a
!! proportional to N above the thermocline (d(N^2)/dz > 0)
!! and to N^2 below the thermocline (d(N^2)/dz < 0)
!!
!! ** Action : av_tide updated in the ITF area (msk_itf)
!!
!! References : Koch-Larrouy et al. 2007, GRL
!!----------------------------------------------------------------------
INTEGER , INTENT(in ) :: kt, Kmm ! ocean time-step
REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) :: pav ! Tidal mixing coef.
!!
INTEGER :: ji, jj, jk ! dummy loop indices
REAL(wp) :: zcoef, ztpc ! temporary scalar
REAL(wp), DIMENSION(jpi,jpj) :: zkz ! 2D workspace
REAL(wp), DIMENSION(jpi,jpj) :: zsum1 , zsum2 , zsum ! - -
REAL(wp), DIMENSION(jpi,jpj,jpk) :: zempba_3d_1, zempba_3d_2 ! 3D workspace
REAL(wp), DIMENSION(jpi,jpj,jpk) :: zempba_3d , zdn2dz ! - -
REAL(wp), DIMENSION(jpi,jpj,jpk) :: zavt_itf ! - -
!!----------------------------------------------------------------------
!
! ! compute the form function using N2 at each time step
zdn2dz (:,:,jpk) = 0.e0
zempba_3d_1(:,:,jpk) = 0.e0
zempba_3d_2(:,:,jpk) = 0.e0
DO jk = 1, jpkm1
zdn2dz (:,:,jk) = rn2(:,:,jk) - rn2(:,:,jk+1) ! Vertical profile of dN2/dz
!CDIR NOVERRCHK
zempba_3d_1(:,:,jk) = SQRT( MAX( 0.e0, rn2(:,:,jk) ) ) ! - - of N
zempba_3d_2(:,:,jk) = MAX( 0.e0, rn2(:,:,jk) ) ! - - of N^2
END DO
!
zsum (:,:) = 0.e0
zsum1(:,:) = 0.e0
zsum2(:,:) = 0.e0
DO jk= 2, jpk
zsum1(:,:) = zsum1(:,:) + zempba_3d_1(:,:,jk) * e3w(:,:,jk,Kmm) * tmask(:,:,jk) * tmask(:,:,jk-1)
zsum2(:,:) = zsum2(:,:) + zempba_3d_2(:,:,jk) * e3w(:,:,jk,Kmm) * tmask(:,:,jk) * tmask(:,:,jk-1)
END DO
DO jj = 1, jpj
DO ji = 1, jpi
IF( zsum1(ji,jj) /= 0.e0 ) zsum1(ji,jj) = 1.e0 / zsum1(ji,jj)
IF( zsum2(ji,jj) /= 0.e0 ) zsum2(ji,jj) = 1.e0 / zsum2(ji,jj)
END DO
END DO
DO jk= 1, jpk
DO jj = 1, jpj
DO ji = 1, jpi
zcoef = 0.5 - SIGN( 0.5, zdn2dz(ji,jj,jk) ) ! =0 if dN2/dz > 0, =1 otherwise
ztpc = zempba_3d_1(ji,jj,jk) * zsum1(ji,jj) * zcoef &
& + zempba_3d_2(ji,jj,jk) * zsum2(ji,jj) * ( 1. - zcoef )
!
zempba_3d(ji,jj,jk) = ztpc
zsum (ji,jj) = zsum(ji,jj) + ztpc * e3w(ji,jj,jk,Kmm)
END DO
END DO
END DO
DO jj = 1, jpj
DO ji = 1, jpi
IF( zsum(ji,jj) > 0.e0 ) zsum(ji,jj) = 1.e0 / zsum(ji,jj)
END DO
END DO
! ! first estimation bounded by 10 cm2/s (with n2 bounded by rn_n2min)
zcoef = rn_tfe_itf / ( rn_tfe * rho0 )
DO jk = 1, jpk
zavt_itf(:,:,jk) = MIN( 10.e-4, zcoef * en_tmx(:,:) * zsum(:,:) * zempba_3d(:,:,jk) &
& / MAX( rn_n2min, rn2(:,:,jk) ) * tmask(:,:,jk) )
END DO
zkz(:,:) = 0.e0 ! Associated potential energy consummed over the whole water column
DO jk = 2, jpkm1
zkz(:,:) = zkz(:,:) + e3w(:,:,jk,Kmm) * MAX( 0.e0, rn2(:,:,jk) ) * rho0 * zavt_itf(:,:,jk) * tmask(:,:,jk) * tmask(:,:,jk-1)
END DO
DO jj = 1, jpj ! Here zkz should be equal to en_tmx ==> multiply by en_tmx/zkz to recover en_tmx
DO ji = 1, jpi
IF( zkz(ji,jj) /= 0.e0 ) zkz(ji,jj) = en_tmx(ji,jj) * rn_tfe_itf / rn_tfe / zkz(ji,jj)
END DO
END DO
DO jk = 2, jpkm1 ! Mutiply by zkz to recover en_tmx, BUT bound by 30/6 ==> zavt_itf bound by 300 cm2/s
zavt_itf(:,:,jk) = zavt_itf(:,:,jk) * MIN( zkz(:,:), 120./10. ) * tmask(:,:,jk) * tmask(:,:,jk-1) ! kz max = 120 cm2/s
END DO
IF( kt == nit000 ) THEN ! diagnose the nergy consumed by zavt_itf
ztpc = 0.e0
DO jk= 1, jpk
DO jj= 1, jpj
DO ji= 1, jpi
ztpc = ztpc + e1t(ji,jj) * e2t(ji,jj) * e3w(ji,jj,jk,Kmm) * MAX( 0.e0, rn2(ji,jj,jk) ) &
& * zavt_itf(ji,jj,jk) * tmask(ji,jj,jk) * tmask_i(ji,jj)
END DO
END DO
END DO
ztpc= rho0 * ztpc / ( rn_me * rn_tfe_itf )
IF(lwp) WRITE(numout,*) ' N Total power consumption by zavt_itf: ztpc = ', ztpc * 1.e-12 ,'TW'
ENDIF
! ! Update pav with the ITF mixing coefficient
DO jk = 2, jpkm1
pav(:,:,jk) = pav (:,:,jk) * ( 1.e0 - mask_itf(:,:) ) &
& + zavt_itf(:,:,jk) * mask_itf(:,:)
END DO
!
END SUBROUTINE tmx_itf
SUBROUTINE zdf_tmx_init
!!----------------------------------------------------------------------
!! *** ROUTINE zdf_tmx_init ***
!!
!! ** Purpose : Initialization of the vertical tidal mixing, Reading
!! of M2 and K1 tidal energy in nc files
!!
!! ** Method : - Read the namtmx namelist and check the parameters
!!
!! - Read the input data in NetCDF files :
!! M2 and K1 tidal energy. The total tidal energy, en_tmx,
!! is the sum of M2, K1 and S2 energy where S2 is assumed
!! to be: S2=(1/2)^2 * M2
!! mask_itf, a mask array that determine where substituing
!! the standard Simmons et al. (2005) formulation with the
!! one of Koch_Larrouy et al. (2007).
!!
!! - Compute az_tmx, a 3D coefficient that allows to compute
!! the standard tidal-induced vertical mixing as follows:
!! Kz_tides = az_tmx / max( rn_n2min, N^2 )
!! with az_tmx a bottom intensified coefficient is given by:
!! az_tmx(z) = en_tmx / ( rho0 * rn_htmx ) * EXP( -(H-z)/rn_htmx )
!! / ( 1. - EXP( - H /rn_htmx ) )
!! where rn_htmx the characteristic length scale of the bottom
!! intensification, en_tmx the tidal energy, and H the ocean depth
!!
!! ** input : - Namlist namtmx
!! - NetCDF file : M2_ORCA2.nc, K1_ORCA2.nc, and mask_itf.nc
!!
!! ** Action : - Increase by 1 the nstop flag is setting problem encounter
!! - defined az_tmx used to compute tidal-induced mixing
!!
!! References : Simmons et al. 2004, Ocean Modelling, 6, 3-4, 245-263.
!! Koch-Larrouy et al. 2007, GRL.
!!----------------------------------------------------------------------
INTEGER :: ji, jj, jk ! dummy loop indices
INTEGER :: inum ! local integer
INTEGER :: ios
REAL(wp) :: ztpc, ze_z ! local scalars
REAL(wp), DIMENSION(jpi,jpj) :: zem2, zek1 ! read M2 and K1 tidal energy
REAL(wp), DIMENSION(jpi,jpj) :: zkz ! total M2, K1 and S2 tidal energy
REAL(wp), DIMENSION(jpi,jpj) :: zfact ! used for vertical structure function
REAL(wp), DIMENSION(jpi,jpj) :: zhdep ! Ocean depth
REAL(wp), DIMENSION(jpi,jpj,jpk) :: zpc ! power consumption
REAL(wp), DIMENSION(jpi,jpj,jpk) :: zav_tide ! tidal mixing coefficient
!!
NAMELIST/namzdf_tmx/ rn_htmx, rn_n2min, rn_tfe, rn_me, ln_tmx_itf, rn_tfe_itf
!!----------------------------------------------------------------------
!
! Namelist namzdf_tmx in reference namelist : Tidal Mixing
READ ( numnam_ref, namzdf_tmx, IOSTAT = ios, ERR = 901)
901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf_tmx in reference namelist' )
! Namelist namzdf_tmx in configuration namelist : Tidal Mixing
READ ( numnam_cfg, namzdf_tmx, IOSTAT = ios, ERR = 902 )
902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf_tmx in configuration namelist' )
IF(lwm) WRITE ( numond, namzdf_tmx )
IF(lwp) THEN ! Control print
WRITE(numout,*)
WRITE(numout,*) 'zdf_tmx_init : tidal mixing'
WRITE(numout,*) '~~~~~~~~~~~~'
WRITE(numout,*) ' Namelist namzdf_tmx : set tidal mixing parameters'
WRITE(numout,*) ' Vertical decay scale for turbulence = ', rn_htmx
WRITE(numout,*) ' Brunt-Vaisala frequency threshold = ', rn_n2min
WRITE(numout,*) ' Tidal dissipation efficiency = ', rn_tfe
WRITE(numout,*) ' Mixing efficiency = ', rn_me
WRITE(numout,*) ' ITF specific parameterisation = ', ln_tmx_itf
WRITE(numout,*) ' ITF tidal dissipation efficiency = ', rn_tfe_itf
ENDIF
! ! allocate tmx arrays
IF( zdf_tmx_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'zdf_tmx_init : unable to allocate tmx arrays' )
IF( ln_tmx_itf ) THEN ! read the Indonesian Through Flow mask
CALL iom_open('mask_itf',inum)
CALL iom_get (inum, jpdom_global, 'tmaskitf',mask_itf,1) !
CALL iom_close(inum)
ENDIF
! read M2 tidal energy flux : W/m2 ( zem2 < 0 )
CALL iom_open('M2rowdrg',inum)
CALL iom_get (inum, jpdom_global, 'field',zem2,1) !
CALL iom_close(inum)
! read K1 tidal energy flux : W/m2 ( zek1 < 0 )
CALL iom_open('K1rowdrg',inum)
CALL iom_get (inum, jpdom_global, 'field',zek1,1) !
CALL iom_close(inum)
! Total tidal energy ( M2, S2 and K1 with S2=(1/2)^2 * M2 )
! only the energy available for mixing is taken into account,
! (mixing efficiency tidal dissipation efficiency)
en_tmx(:,:) = - rn_tfe * rn_me * ( zem2(:,:) * 1.25 + zek1(:,:) ) * ssmask(:,:)
!============
!TG: Bug for VVL? Should this section be moved out of _init and be updated at every timestep?
! Vertical structure (az_tmx)
DO jj = 1, jpj ! part independent of the level
DO ji = 1, jpi
zhdep(ji,jj) = gdepw_0(ji,jj,mbkt(ji,jj)+1) ! depth of the ocean
zfact(ji,jj) = rho0 * rn_htmx * ( 1. - EXP( -zhdep(ji,jj) / rn_htmx ) )
IF( zfact(ji,jj) /= 0 ) zfact(ji,jj) = en_tmx(ji,jj) / zfact(ji,jj)
END DO
END DO
DO jk= 1, jpk ! complete with the level-dependent part
DO jj = 1, jpj
DO ji = 1, jpi
az_tmx(ji,jj,jk) = zfact(ji,jj) * EXP( -( zhdep(ji,jj)-gdepw_0(ji,jj,jk) ) / rn_htmx ) * tmask(ji,jj,jk)
END DO
END DO
END DO
!===========
IF( lwp ) THEN
! Control print
! Total power consumption due to vertical mixing
! zpc = rho0 * 1/rn_me * rn2 * zav_tide
zav_tide(:,:,:) = 0.e0
DO jk = 2, jpkm1
zav_tide(:,:,jk) = az_tmx(:,:,jk) / MAX( rn_n2min, rn2(:,:,jk) )
END DO
ztpc = 0.e0
zpc(:,:,:) = MAX(rn_n2min,rn2(:,:,:)) * zav_tide(:,:,:)
DO jk= 2, jpkm1
DO jj = 1, jpj
DO ji = 1, jpi
ztpc = ztpc + e3w_0(ji,jj,jk) * e1t(ji,jj) * e2t(ji,jj) * zpc(ji,jj,jk) * wmask(ji,jj,jk) * tmask_i(ji,jj)
END DO
END DO
END DO
ztpc= rho0 * 1/(rn_tfe * rn_me) * ztpc
WRITE(numout,*)
WRITE(numout,*) ' Total power consumption of the tidally driven part of Kz : ztpc = ', ztpc * 1.e-12 ,'TW'
! control print 2
zav_tide(:,:,:) = MIN( zav_tide(:,:,:), 60.e-4 )
zkz(:,:) = 0.e0
DO jk = 2, jpkm1
DO jj = 1, jpj
DO ji = 1, jpi
zkz(ji,jj) = zkz(ji,jj) + e3w_0(ji,jj,jk) * MAX(0.e0, rn2(ji,jj,jk)) * rho0 * zav_tide(ji,jj,jk) * wmask(ji,jj,jk)
END DO
END DO
END DO
! Here zkz should be equal to en_tmx ==> multiply by en_tmx/zkz
DO jj = 1, jpj
DO ji = 1, jpi
IF( zkz(ji,jj) /= 0.e0 ) THEN
zkz(ji,jj) = en_tmx(ji,jj) / zkz(ji,jj)
ENDIF
END DO
END DO
ztpc = 1.e50
DO jj = 1, jpj
DO ji = 1, jpi
IF( zkz(ji,jj) /= 0.e0 ) THEN
ztpc = Min( zkz(ji,jj), ztpc)
ENDIF
END DO
END DO
WRITE(numout,*) ' Min de zkz ', ztpc, ' Max = ', maxval(zkz(:,:) )
DO jk = 2, jpkm1
DO jj = 1, jpj
DO ji = 1, jpi
zav_tide(ji,jj,jk) = zav_tide(ji,jj,jk) * MIN( zkz(ji,jj), 30./6. ) * wmask(ji,jj,jk) !kz max = 300 cm2/s
END DO
END DO
END DO
ztpc = 0.e0
zpc(:,:,:) = Max(0.e0,rn2(:,:,:)) * zav_tide(:,:,:)
DO jk= 1, jpk
DO jj = 1, jpj
DO ji = 1, jpi
ztpc = ztpc + e3w_0(ji,jj,jk) * e1t(ji,jj) * e2t(ji,jj) * zpc(ji,jj,jk) * wmask(ji,jj,jk) * tmask_i(ji,jj)
END DO
END DO
END DO
ztpc= rho0 * 1/(rn_tfe * rn_me) * ztpc
WRITE(numout,*) ' 2 Total power consumption of the tidally driven part of Kz : ztpc = ', ztpc * 1.e-12 ,'TW'
DO jk = 1, jpk
ze_z = SUM( e1t(:,:) * e2t(:,:) * zav_tide(:,:,jk) * tmask_i(:,:) ) &
& / MAX( 1.e-20, SUM( e1t(:,:) * e2t(:,:) * wmask (:,:,jk) * tmask_i(:,:) ) )
ztpc = 1.E50
DO jj = 1, jpj
DO ji = 1, jpi
IF( zav_tide(ji,jj,jk) /= 0.e0 ) ztpc =Min( ztpc, zav_tide(ji,jj,jk) )
END DO
END DO
WRITE(numout,*) ' N2 min - jk= ', jk,' ', ze_z * 1.e4,' cm2/s min= ',ztpc*1.e4, &
& 'max= ', MAXVAL(zav_tide(:,:,jk) )*1.e4, ' cm2/s'
END DO
WRITE(numout,*) ' e_tide : ', SUM( e1t*e2t*en_tmx ) / ( rn_tfe * rn_me ) * 1.e-12, 'TW'
WRITE(numout,*)
WRITE(numout,*) ' Initial profile of tidal vertical mixing'
DO jk = 1, jpk
DO jj = 1,jpj
DO ji = 1,jpi
zkz(ji,jj) = az_tmx(ji,jj,jk) /MAX( rn_n2min, rn2(ji,jj,jk) )
END DO
END DO
ze_z = SUM( e1t(:,:) * e2t(:,:) * zkz(:,:) * tmask_i(:,:) ) &
& / MAX( 1.e-20, SUM( e1t(:,:) * e2t(:,:) * wmask (:,:,jk) * tmask_i(:,:) ) )
WRITE(numout,*) ' jk= ', jk,' ', ze_z * 1.e4,' cm2/s'
END DO
DO jk = 1, jpk
zkz(:,:) = az_tmx(:,:,jk) /rn_n2min
ze_z = SUM( e1t(:,:) * e2t(:,:) * zkz(:,:) * tmask_i(:,:) ) &
& / MAX( 1.e-20, SUM( e1t(:,:) * e2t(:,:) * wmask (:,:,jk) * tmask_i(:,:) ) )
WRITE(numout,*)
WRITE(numout,*) ' N2 min - jk= ', jk,' ', ze_z * 1.e4,' cm2/s min= ',MINVAL(zkz)*1.e4, &
& 'max= ', MAXVAL(zkz)*1.e4, ' cm2/s'
END DO
!
ENDIF
!
END SUBROUTINE zdf_tmx_init
!!======================================================================
END MODULE zdftmx
src/OCE/nemogcm.F90
View file @ 4a9f4b18
......@@ -73,6 +73,7 @@ MODULE nemogcm
USE lib_mpp ! distributed memory computing
USE mppini ! shared/distributed memory setting (mpp_init routine)
USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)
USE sbccpl
USE halo_mng ! halo manager
IMPLICIT NONE
......@@ -173,6 +174,9 @@ CONTAINS
IF ( istp == nitend ) elapsed_time = zstptiming - elapsed_time
ENDIF
!
IF (lk_oasis) THEN
CALL sbc_cpl_snd( istp, Nbb, Nnn ) ! Coupling to atmos
ENDIF
# if defined key_qco || defined key_linssh
CALL stp_MLF( istp )
# else
......@@ -262,7 +266,7 @@ CONTAINS
#if defined key_xios
IF( Agrif_Root() ) THEN
IF( lk_oasis ) THEN
CALL cpl_init( "oceanx", ilocal_comm ) ! nemo local communicator given by oasis
CALL cpl_init( "toyoce", ilocal_comm ) ! nemo local communicator given by oasis
CALL xios_initialize( "not used" , local_comm =ilocal_comm ) ! send nemo communicator to xios
ELSE
CALL xios_initialize( "for_xios_mpi_id", return_comm=ilocal_comm ) ! nemo local communicator given by xios
......@@ -272,7 +276,7 @@ CONTAINS
#else
IF( lk_oasis ) THEN
IF( Agrif_Root() ) THEN
CALL cpl_init( "oceanx", ilocal_comm ) ! nemo local communicator given by oasis
CALL cpl_init( "toyoce", ilocal_comm ) ! nemo local communicator given by oasis
ENDIF
CALL mpp_start( ilocal_comm )
ELSE
......@@ -496,6 +500,11 @@ CONTAINS
!
IF(lwp) WRITE(numout,cform_aaa) ! Flag AAAAAAA
!
IF (nstop > 0) THEN
CALL CTL_STOP('STOP','Critical errors in NEMO initialisation')
END IF
IF( ln_timing ) CALL timing_stop( 'nemo_init')
!
END SUBROUTINE nemo_init
......
src/OCE/par_oce.F90
View file @ 4a9f4b18
......@@ -99,6 +99,12 @@ MODULE par_oce
INTEGER, PUBLIC :: Ni_0, Nj_0 !: local domain size without halo
INTEGER, PUBLIC :: Ni0glo, Nj0glo !: global domain size without halo
INTEGER, PUBLIC :: Nis0_ext !: start I-index with wrap/N-fold
INTEGER, PUBLIC :: Nie0_ext !: end I-index with wrap/N-fold
INTEGER, PUBLIC :: Njs0_ext !: start J-index with wrap/N-fold
INTEGER, PUBLIC :: Nje0_ext !: end J-index with wrap/N-fold
INTEGER, PUBLIC :: Ni_0_ext, Nj_0_ext !: local domain size with wrap/N-fold
INTEGER, PUBLIC :: Ni0glo_ext, Nj0glo_ext !: global domain size with wrap/N-fold
!!----------------------------------------------------------------------
!! NEMO/OCE 4.0 , NEMO Consortium (2018)
!! $Id: par_oce.F90 15119 2021-07-13 14:43:22Z jchanut $
......
src/OCE/step.F90
View file @ 4a9f4b18
......@@ -172,6 +172,7 @@ CONTAINS
CALL bn2 ( ts(:,:,:,:,Nbb), rab_b, rn2b, Nnn ) ! before Brunt-Vaisala frequency
CALL bn2 ( ts(:,:,:,:,Nnn), rab_n, rn2, Nnn ) ! now Brunt-Vaisala frequency
! VERTICAL PHYSICS
! lbc_lnk needed for zdf_sh2 when using nn_hls = 2, moved here to allow tiling in zdf_phy
IF( nn_hls == 2 .AND. l_zdfsh2 ) CALL lbc_lnk( 'stp', avm_k, 'W', 1.0_wp )
......@@ -296,6 +297,7 @@ CONTAINS
IF( ln_diadct ) CALL dia_dct ( kstp, Nnn ) ! Transports
CALL dia_ar5 ( kstp, Nnn ) ! ar5 diag
CALL dia_ptr ( kstp, Nnn ) ! Poleward adv/ldf TRansports diagnostics
CALL dia_prod ( kstp, Nnn ) ! ocean model: products
CALL dia_wri ( kstp, Nnn ) ! ocean model: outputs
IF( ln_crs ) CALL crs_fld ( kstp, Nnn ) ! ocean model: online field coarsening & output
IF( lk_diadetide ) CALL dia_detide( kstp ) ! Weights computation for daily detiding of model diagnostics
......@@ -429,7 +431,7 @@ CONTAINS
!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
! Coupled mode
!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
IF( lk_oasis .AND. nstop == 0 ) CALL sbc_cpl_snd( kstp, Nbb, Nnn ) ! coupled mode : field exchanges
!!IF( lk_oasis .AND. nstop == 0 ) CALL sbc_cpl_snd( kstp, Nbb, Nnn ) ! coupled mode : field exchanges
!
#if defined key_xios
!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
......
src/OCE/step_oce.F90
View file @ 4a9f4b18
......@@ -85,6 +85,7 @@ MODULE step_oce
USE diahsb ! heat, salt and volume budgets (dia_hsb routine)
USE diacfl ! CFL diagnostics (dia_cfl routine)
USE diaobs ! Observation operator (dia_obs routine)
USE diaprod
USE diadetide ! Weights computation for daily detiding of model diagnostics
USE diamlr ! IOM context management for multiple-linear-regression analysis
USE flo_oce ! floats variables
......
src/OCE/stpmlf.F90
View file @ 4a9f4b18
......@@ -443,7 +443,7 @@ CONTAINS
!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
! Coupled mode
!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
IF( lk_oasis .AND. nstop == 0 ) CALL sbc_cpl_snd( kstp, Nbb, Nnn ) ! coupled mode : field exchanges
!!IF( lk_oasis .AND. nstop == 0 ) CALL sbc_cpl_snd( kstp, Nbb, Nnn ) ! coupled mode : field exchanges
!
#if defined key_xios
!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
......@@ -561,6 +561,7 @@ CONTAINS
CALL lbc_lnk( 'finalize_lbc', r3u(:,:,Kaa), 'U', 1._wp, r3v(:,:,Kaa), 'V', 1._wp, &
& r3u_f(:,:), 'U', 1._wp, r3v_f(:,:), 'V', 1._wp )
ENDIF
! !* BDY open boundaries
IF( ln_bdy ) THEN
CALL bdy_tra( kt, Kbb, pts, Kaa )
......
src/SAS/sbcssm.F90
View file @ 4a9f4b18
......@@ -79,9 +79,19 @@ CONTAINS
INTEGER :: ji, jj ! dummy loop indices
REAL(wp) :: ztinta ! ratio applied to after records when doing time interpolation
REAL(wp) :: ztintb ! ratio applied to before records when doing time interpolation
CHARACTER(len=4),SAVE :: stype
!!----------------------------------------------------------------------
!
IF( ln_timing ) CALL timing_start( 'sbc_ssm')
IF( kt == nit000 ) THEN
IF( ln_TEOS10 ) THEN
stype='abs' ! teos-10: using absolute salinity (sst is converted to potential temperature for the surface module)
ELSE IF( ln_EOS80 ) THEN
stype='pra' ! eos-80: using practical salinity
ELSE IF ( ln_SEOS) THEN
stype='seos' ! seos using Simplified Equation of state (sst is converted to potential temperature for the surface module)
ENDIF
ENDIF
IF ( l_sasread ) THEN
IF( nfld_3d > 0 ) CALL fld_read( kt, 1, sf_ssm_3d ) !== read data at kt time step ==!
......@@ -156,8 +166,8 @@ CONTAINS
IF( l_initdone ) THEN ! Mean value at each nn_fsbc time-step !
CALL iom_put( 'ssu_m', ssu_m )
CALL iom_put( 'ssv_m', ssv_m )
CALL iom_put( 'sst_m', sst_m )
CALL iom_put( 'sss_m', sss_m )
CALL iom_put( 'sst_m_pot', sst_m )
CALL iom_put( 'sss_m_'//stype, sss_m )
CALL iom_put( 'ssh_m', ssh_m )
CALL iom_put( 'e3t_m', e3t_m )
IF( ln_read_frq ) CALL iom_put( 'frq_m', frq_m )
......