MODULE domvvl !!====================================================================== !! *** MODULE domvvl *** !! Ocean : !!====================================================================== !! History : 5.0 ! 2018-07 (G. Madec) Flux Form with Kinetic Energy conservation !! ==>>> here z* and s* only (no z-tilde) ! 1- remove z-tilde ==>>> pure z-star (or s-star) ! 2- remove dom_vvl_interpol !!---------------------------------------------------------------------- !!---------------------------------------------------------------------- !! dom_vvl_init : define initial vertical scale factors, depths and column thickness !! dom_vvl_sf_nxt : Compute next vertical scale factors !! dom_vvl_sf_swp : Swap vertical scale factors and update the vertical grid !! dom_vvl_rst : read/write restart file !! dom_vvl_ctl : Check the vvl options !!---------------------------------------------------------------------- USE oce ! ocean dynamics and tracers USE phycst ! physical constant USE dom_oce ! ocean space and time domain USE sbc_oce ! ocean surface boundary condition USE wet_dry ! wetting and drying USE usrdef_istate ! user defined initial state (wad only) USE restart ! ocean restart ! USE in_out_manager ! I/O manager USE iom ! I/O manager library USE lib_mpp ! distributed memory computing library USE lbclnk ! ocean lateral boundary conditions (or mpp link) USE timing ! Timing IMPLICIT NONE PRIVATE PUBLIC dom_vvl_init ! called by domain.F90 PUBLIC dom_vvl_sf_nxt ! called by step.F90 PUBLIC dom_vvl_sf_swp ! called by step.F90 PUBLIC ssh2e3_before ! ... PUBLIC ssh2e3_now ! ... ! !!* Namelist nam_vvl LOGICAL , PUBLIC :: ln_vvl_zstar = .FALSE. ! zstar vertical coordinate LOGICAL , PUBLIC :: ln_vvl_ztilde = .FALSE. ! ztilde vertical coordinate LOGICAL , PUBLIC :: ln_vvl_layer = .FALSE. ! level vertical coordinate LOGICAL , PUBLIC :: ln_vvl_ztilde_as_zstar = .FALSE. ! ztilde vertical coordinate LOGICAL , PUBLIC :: ln_vvl_zstar_at_eqtor = .FALSE. ! ztilde vertical coordinate LOGICAL , PUBLIC :: ln_vvl_kepe = .FALSE. ! kinetic/potential energy transfer ! ! conservation: not used yet REAL(wp) :: rn_ahe3 ! thickness diffusion coefficient REAL(wp) :: rn_rst_e3t ! ztilde to zstar restoration timescale [days] REAL(wp) :: rn_lf_cutoff ! cutoff frequency for low-pass filter [days] REAL(wp) :: rn_zdef_max ! maximum fractional e3t deformation LOGICAL , PUBLIC :: ln_vvl_dbg = .FALSE. ! debug control prints REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: un_td, vn_td ! thickness diffusion transport REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: hdiv_lf ! low frequency part of hz divergence REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: te3t_b, te3t_n ! baroclinic scale factors REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: te3t_a, dte3t_a ! baroclinic scale factors REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:) :: frq_rst_e3t ! retoring period for scale factors REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:) :: frq_rst_hdv ! retoring period for low freq. divergence !!gm add !!gm !! * Substitutions # include "vectopt_loop_substitute.h90" !!---------------------------------------------------------------------- !! NEMO/OCE 4.0 , NEMO Consortium (2018) !! $Id$ !! Software governed by the CeCILL licence (./LICENSE) !!---------------------------------------------------------------------- CONTAINS INTEGER FUNCTION dom_vvl_alloc() !!---------------------------------------------------------------------- !! *** FUNCTION dom_vvl_alloc *** !!---------------------------------------------------------------------- ! IF( ln_vvl_zstar ) dom_vvl_alloc = 0 IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN ALLOCATE( te3t_b(jpi,jpj,jpk) , te3t_n(jpi,jpj,jpk) , te3t_a(jpi,jpj,jpk) , & & dte3t_a(jpi,jpj,jpk) , un_td (jpi,jpj,jpk) , vn_td (jpi,jpj,jpk) , & & STAT = dom_vvl_alloc ) IF( lk_mpp ) CALL mpp_sum ( dom_vvl_alloc ) IF( dom_vvl_alloc /= 0 ) CALL ctl_warn('dom_vvl_alloc: failed to allocate arrays') un_td = 0._wp vn_td = 0._wp ENDIF IF( ln_vvl_ztilde ) THEN ALLOCATE( frq_rst_e3t(jpi,jpj) , frq_rst_hdv(jpi,jpj) , hdiv_lf(jpi,jpj,jpk) , STAT= dom_vvl_alloc ) IF( lk_mpp ) CALL mpp_sum ( dom_vvl_alloc ) IF( dom_vvl_alloc /= 0 ) CALL ctl_warn('dom_vvl_alloc: failed to allocate arrays') ENDIF ! END FUNCTION dom_vvl_alloc SUBROUTINE dom_vvl_init !!---------------------------------------------------------------------- !! *** ROUTINE dom_vvl_init *** !! !! ** Purpose : Initialization of all scale factors, depths !! and water column heights !! !! ** Method : - use restart file and/or initialize !! - interpolate scale factors !! !! ** Action : - e3t_(n/b) and te3t_(n/b) !! - Regrid: e3(u/v)_n !! e3(u/v)_b !! e3w_n !! e3(u/v)w_b !! e3(u/v)w_n !! gdept_n, gdepw_n and gde3w_n !! - h(t/u/v)_0 !! - frq_rst_e3t and frq_rst_hdv !! !! Reference : Leclair, M., and G. Madec, 2011, Ocean Modelling. !!---------------------------------------------------------------------- !!---------------------------------------------------------------------- ! IF(lwp) WRITE(numout,*) IF(lwp) WRITE(numout,*) 'dom_vvl_init : Variable volume activated' IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~' ! CALL dom_vvl_ctl ! choose vertical coordinate (z_star, z_tilde or layer) ! ! ! Allocate module arrays IF( dom_vvl_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'dom_vvl_init : unable to allocate arrays' ) ! ! ! Read or initialize e3t_(b/n), ssh(b/n) CALL dom_vvl_rst( nit000, 'READ' ) ! ! !== Set of all other vertical mesh fields ==! (now and before) ! ! !* BEFORE fields : CALL ssh2e3_before ! set: hu , hv , r1_hu, r1_hv ! ! e3t, e3w, e3u, e3uw, e3v, e3vw ! ! ! set one for all last level to the e3._0 value e3t_b(:,:,jpk) = e3t_0(:,:,jpk) ; e3u_b(:,:,jpk) = e3w_0(:,:,jpk) ; e3v_b(:,:,jpk) = e3v_0(:,:,jpk) e3w_b(:,:,jpk) = e3w_0(:,:,jpk) ; e3uw_b(:,:,jpk) = e3uw_0(:,:,jpk) ; e3vw_b(:,:,jpk) = e3vw_0(:,:,jpk) ! ! !* NOW fields : CALL ssh2e3_now ! set: ht , hu , hv , r1_hu, r1_hv ! ! e3t, e3w, e3u, e3uw, e3v, e3vw, e3f ! ! gdept_n, gdepw_n, gde3w_n ! ! ! set one for all last level to the e3._0 value e3t_n(:,:,jpk) = e3t_0(:,:,jpk) ; e3u_n(:,:,jpk) = e3w_0(:,:,jpk) ; e3v_n(:,:,jpk) = e3v_0(:,:,jpk) e3w_n(:,:,jpk) = e3w_0(:,:,jpk) ; e3uw_n(:,:,jpk) = e3uw_0(:,:,jpk) ; e3vw_n(:,:,jpk) = e3vw_0(:,:,jpk) e3f_n(:,:,jpk) = e3f_0(:,:,jpk) ! ! !* AFTER fields : (last level for OPA, 3D required for AGRIF initialisation) e3t_a(:,:,:) = e3t_n(:,:,:) ; e3u_a(:,:,:) = e3u_n(:,:,:) ; e3v_a(:,:,:) = e3v_n(:,:,:) !!gm !!gm ===>>>> below: some issues to think about !!! !!gm !!gm fmask definition checked (0 or 1 nothing else) !! in z-tilde or ALE e3._0 should be the time varying fields step forward with an euler scheme !!gm e3w_b & gdept_b are not used except its update in agrif !! and used to compute before slope of surface in dynldf_iso ==>>> remove it !!!! !! NB: in triads on TRA, gdept_n is used !!!! BUG? !!gm e3f_n almost not used ===>>>> verify whether it can be removed or not... !! verify the use of wumask & wvmask mau be replaced by a multiplication by umask(k)*umask(k+1) ??? !! !!gm ISF case to be checked by Pierre Mathiot !! !!gm setting of e3._a for agrif.... TO BE CHECKED ! IF(lwxios) THEN ! define variables in restart file when writing with XIOS CALL iom_set_rstw_var_active('e3t_b') CALL iom_set_rstw_var_active('e3t_n') ENDIF ! END SUBROUTINE dom_vvl_init SUBROUTINE dom_vvl_sf_nxt( kt, kcall ) !!---------------------------------------------------------------------- !! *** ROUTINE dom_vvl_sf_nxt *** !! !! ** Purpose : - compute the after scale factors used in tra_zdf, dynnxt, !! tranxt and dynspg routines !! !! ** Method : - z_star case: Repartition of ssh INCREMENT proportionnaly to the level thickness. !! - z_tilde_case: after scale factor increment = !! high frequency part of horizontal divergence !! + retsoring towards the background grid !! + thickness difusion !! Then repartition of ssh INCREMENT proportionnaly !! to the "baroclinic" level thickness. !! !! ** Action : - hdiv_lf : restoring towards full baroclinic divergence in z_tilde case !! - te3t_a: after increment of vertical scale factor !! in z_tilde case !! - e3(t/u/v)_a !! !! Reference : Leclair, M., and Madec, G. 2011, Ocean Modelling. !!---------------------------------------------------------------------- INTEGER, INTENT( in ) :: kt ! time step INTEGER, INTENT( in ), OPTIONAL :: kcall ! optional argument indicating call sequence ! INTEGER :: ji, jj, jk ! dummy loop indices REAL(wp), DIMENSION(jpi,jpj) :: zssht_h, zsshu_h, zsshv_h !!---------------------------------------------------------------------- ! IF( ln_linssh ) RETURN ! No calculation in linear free surface ! IF( ln_timing ) CALL timing_start('dom_vvl_sf_nxt') ! IF( kt == nit000 ) THEN IF(lwp) WRITE(numout,*) IF(lwp) WRITE(numout,*) 'dom_vvl_sf_nxt : compute after scale factors' IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~' ENDIF ! ! ------------------! ! ! z_star coordinate ! (and barotropic z-tilde part) ! ! ------------------! ! ! !== after ssh ==! (u- and v-points) DO jj = 2, jpjm1 ; DO ji = 2, jpim1 zsshu_h(ji,jj) = 0.5_wp * ( ssha(ji,jj) + ssha(ji+1,jj) ) * ssumask(ji,jj) zsshv_h(ji,jj) = 0.5_wp * ( ssha(ji,jj) + ssha(ji,jj+1) ) * ssvmask(ji,jj) END DO ; END DO CALL lbc_lnk_multi( zsshu_h(:,:), 'U', 1._wp , zsshv_h(:,:), 'V', 1._wp ) ! ! !== after depths and its inverse ==! hu_a(:,:) = hu_0(:,:) + zsshu_h(:,:) hv_a(:,:) = hv_0(:,:) + zsshv_h(:,:) r1_hu_a(:,:) = ssumask(:,:) / ( hu_a(:,:) + 1._wp - ssumask(:,:) ) r1_hv_a(:,:) = ssvmask(:,:) / ( hv_a(:,:) + 1._wp - ssvmask(:,:) ) ! ! !== after scale factors ==! (e3t , e3u , e3v) zssht_h(:,:) = ssha (:,:) * r1_ht_0(:,:) ! t-point zsshu_h(:,:) = zsshu_h(:,:) * r1_hu_0(:,:) ! u-point zsshv_h(:,:) = zsshv_h(:,:) * r1_hv_0(:,:) ! v-point DO jk = 1, jpkm1 e3t_a(:,:,jk) = e3t_0(:,:,jk) * ( 1._wp + zssht_h(:,:) * tmask(:,:,jk) ) e3u_a(:,:,jk) = e3u_0(:,:,jk) * ( 1._wp + zsshu_h(:,:) * umask(:,:,jk) ) e3v_a(:,:,jk) = e3v_0(:,:,jk) * ( 1._wp + zsshv_h(:,:) * vmask(:,:,jk) ) END DO ! IF( ln_timing ) CALL timing_stop('dom_vvl_sf_nxt') ! END SUBROUTINE dom_vvl_sf_nxt SUBROUTINE dom_vvl_sf_swp( kt ) !!---------------------------------------------------------------------- !! *** ROUTINE dom_vvl_sf_swp *** !! !! ** Purpose : compute time filter and swap of scale factors !! compute all depths and related variables for next time step !! write outputs and restart file !! !! ** Method : - swap of e3t with trick for volume/tracer conservation !! - reconstruct scale factor at other grid points (interpolate) !! - recompute depths and water height fields !! !! ** Action : - e3t_(b/n), te3t_(b/n) and e3(u/v)_n ready for next time step !! - Recompute: !! e3(u/v)_b !! e3w_n !! e3(u/v)w_b !! e3(u/v)w_n !! gdept_n, gdepw_n and gde3w_n !! h(u/v) and h(u/v)r !! !! Reference : Leclair, M., and G. Madec, 2009, Ocean Modelling. !! Leclair, M., and G. Madec, 2011, Ocean Modelling. !!---------------------------------------------------------------------- INTEGER, INTENT( in ) :: kt ! time step ! INTEGER :: ji, jj, jk ! dummy loop indices REAL(wp), DIMENSION(jpi,jpj) :: zssht_h, zsshu_h, zsshv_h, zsshf_h !!---------------------------------------------------------------------- ! IF( ln_linssh ) RETURN ! No calculation in linear free surface ! IF( ln_timing ) CALL timing_start('dom_vvl_sf_swp') ! IF( kt == nit000 ) THEN IF(lwp) WRITE(numout,*) IF(lwp) WRITE(numout,*) 'dom_vvl_sf_swp : - time filter and swap of scale factors' IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~ - interpolate scale factors and compute depths for next time step' ENDIF ! ! Swap and Compute all missing vertical scale factor and depths ! ============================================================= ! - ML - e3u_b and e3v_b are allready computed in dynnxt ! - JC - hu_b, hv_b, hur_b, hvr_b also ! ! - GM - to be updated : e3f_n, e3w_n , e3uw_n , e3vw_n ! e3w_b , e3uw_b , e3vw_b ! gdept_n , gdepw_n , gde3w_n ! ht_n ! to be swap : hu_a , hv_a , r1_hu_a , r1_hv_a ! ! Local depth and Inverse of the local depth of the water ! ------------------------------------------------------- ! ! swap of depth and scale factors ! ! =============================== DO jk = 1, jpkm1 ! swap n--> a gdept_b(:,:,jk) = gdept_n(:,:,jk) ! depth at t and w gdepw_b(:,:,jk) = gdepw_n(:,:,jk) e3t_n (:,:,jk) = e3t_a (:,:,jk) ! e3t, e3u, e3v e3u_n (:,:,jk) = e3u_a (:,:,jk) e3v_n (:,:,jk) = e3v_a (:,:,jk) END DO ht_n(:,:) = ht_0(:,:) + sshn(:,:) ! ocean thickness ! hu_n(:,:) = hu_a(:,:) ; r1_hu_n(:,:) = r1_hu_a(:,:) hv_n(:,:) = hv_a(:,:) ; r1_hv_n(:,:) = r1_hv_a(:,:) ! ! !== before : ! !* ssh at u- and v-points) DO jj = 2, jpjm1 ; DO ji = 2, jpim1 zsshu_h(ji,jj) = 0.5_wp * ( sshb(ji ,jj) + sshb(ji+1,jj ) ) * ssumask(ji,jj) zsshv_h(ji,jj) = 0.5_wp * ( sshb(ji ,jj) + sshb(ji ,jj+1) ) * ssvmask(ji,jj) END DO ; END DO CALL lbc_lnk_multi( zsshu_h(:,:),'U', 1._wp , zsshv_h(:,:),'V', 1._wp ) ! ! !* e3w_b , e3uw_b , e3vw_b zssht_h(:,:) = sshb (:,:) * r1_ht_0(:,:) ! w-point zsshu_h(:,:) = zsshu_h(:,:) * r1_hu_0(:,:) ! uw-point zsshv_h(:,:) = zsshv_h(:,:) * r1_hv_0(:,:) ! vw-point DO jk = 1, jpkm1 e3w_b(:,:,jk) = e3w_0(:,:,jk) * ( 1._wp + zssht_h(:,:) * tmask(:,:,jk) ) e3uw_b(:,:,jk) = e3uw_0(:,:,jk) * ( 1._wp + zsshu_h(:,:) * umask(:,:,jk) ) e3vw_b(:,:,jk) = e3vw_0(:,:,jk) * ( 1._wp + zsshv_h(:,:) * vmask(:,:,jk) ) END DO ! ! !== now : ! !* ssh at u- and v-points) DO jj = 1, jpjm1 ; DO ji = 1, jpim1 ! start from 1 for f-point zsshu_h(ji,jj) = 0.50_wp * ( sshn(ji ,jj) + sshn(ji+1,jj ) ) * ssumask(ji,jj) zsshv_h(ji,jj) = 0.50_wp * ( sshn(ji ,jj) + sshn(ji ,jj+1) ) * ssvmask(ji,jj) zsshf_h(ji,jj) = 0.25_wp * ( sshn(ji ,jj) + sshn(ji ,jj+1) & & + sshn(ji+1,jj) + sshn(ji+1,jj+1) ) * ssfmask(ji,jj) END DO ; END DO CALL lbc_lnk_multi( zsshu_h(:,:),'U', 1._wp , zsshv_h(:,:),'V', 1._wp , zsshf_h(:,:),'F', 1._wp ) ! ! !* e3w_n , e3uw_n , e3vw_n, e3f_n zssht_h(:,:) = sshn (:,:) * r1_ht_0(:,:) ! t- & w-point zsshu_h(:,:) = zsshu_h(:,:) * r1_hu_0(:,:) ! uw-point zsshv_h(:,:) = zsshv_h(:,:) * r1_hv_0(:,:) ! vw-point zsshf_h(:,:) = zsshf_h(:,:) * r1_hf_0(:,:) ! f-point DO jk = 1, jpkm1 e3w_n(:,:,jk) = e3w_0(:,:,jk) * ( 1._wp + zssht_h(:,:) * tmask(:,:,jk) ) e3uw_n(:,:,jk) = e3uw_0(:,:,jk) * ( 1._wp + zsshu_h(:,:) * wumask(:,:,jk) ) e3vw_n(:,:,jk) = e3vw_0(:,:,jk) * ( 1._wp + zsshv_h(:,:) * wvmask(:,:,jk) ) e3f_n(:,:,jk) = e3f_0(:,:,jk) * ( 1._wp + zsshf_h(:,:) * fmask(:,:,jk) ) END DO ! zssht_h(:,:) = 1._wp + sshn (:,:) * r1_ht_0(:,:) ! t-point ! IF( ln_isfcav ) THEN ! ISF cavities : ssh scaling not applied over the iceshelf thickness DO jk = 1, jpkm1 gdept_n(:,:,jk) = ( gdept_0(:,:,jk) - risfdep(:,:) ) * zssht_h(:,:) + risfdep(:,:) gdepw_n(:,:,jk) = ( gdepw_0(:,:,jk) - risfdep(:,:) ) * zssht_h(:,:) + risfdep(:,:) gde3w_n(:,:,jk) = gdept_n(:,:,jk) - sshn (:,:) END DO ELSE ! no ISF cavities DO jk = 1, jpkm1 gdept_n(:,:,jk) = gdept_0(:,:,jk) * zssht_h(:,:) gdepw_n(:,:,jk) = gdepw_0(:,:,jk) * zssht_h(:,:) gde3w_n(:,:,jk) = gdept_n(:,:,jk) - sshn(:,:) END DO ENDIF ! ! write restart file ! ================== IF( lrst_oce ) CALL dom_vvl_rst( kt, 'WRITE' ) ! IF( ln_timing ) CALL timing_stop('dom_vvl_sf_swp') ! END SUBROUTINE dom_vvl_sf_swp SUBROUTINE dom_vvl_rst( kt, cdrw ) !!--------------------------------------------------------------------- !! *** ROUTINE dom_vvl_rst *** !! !! ** Purpose : Read or write VVL file in restart file !! !! ** Method : use of IOM library !! if the restart does not contain vertical scale factors, !! they are set to the _0 values !! if the restart does not contain vertical scale factors increments (z_tilde), !! they are set to 0. !!---------------------------------------------------------------------- INTEGER , INTENT(in) :: kt ! ocean time-step CHARACTER(len=*), INTENT(in) :: cdrw ! "READ"/"WRITE" flag ! INTEGER :: ji, jj, jk INTEGER :: id1, id2, id3, id4, id5 ! local integers !!---------------------------------------------------------------------- ! IF( TRIM(cdrw) == 'READ' ) THEN ! Read/initialise ! ! =============== IF( ln_rstart ) THEN !* Read the restart file CALL rst_read_open ! open the restart file if necessary ! id1 = iom_varid( numror, 'sshb' , ldstop = .FALSE. ) id2 = iom_varid( numror, 'sshn' , ldstop = .FALSE. ) ! IF( MIN( id1, id2 ) > 0 ) THEN ! all required arrays exist IF(lwp) write(numout,*) 'dom_vvl_rst : both sshb and sshn found in restart files' ! !!gm Question: use jpdom_data above to read data over jpi x jpj (like is dom_hgr_read and dom_zgr_read) !! so that it will work with land processor suppression ! CALL iom_get( numror, jpdom_autoglo, 'sshn' , sshn, ldxios = lrxios ) ! CALL iom_get( numror, jpdom_autoglo, 'sshb' , sshb, ldxios = lrxios ) !!gm CALL iom_get( numror, jpdom_data, 'sshn' , sshn, ldxios = lrxios ) CALL iom_get( numror, jpdom_data, 'sshb' , sshb, ldxios = lrxios ) !!gm end IF( l_1st_euler ) THEN sshb(:,:) = sshn(:,:) ENDIF ELSE IF( id1 > 0 ) THEN IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : sshn not found in restart files' IF(lwp) write(numout,*) ' set sshn = sshb and force l_1st_euler = true' !!gm CALL iom_get( numror, jpdom_autoglo, 'e3t_b', e3t_b(:,:,:), ldxios = lrxios ) CALL iom_get( numror, jpdom_data, 'sshb', sshb, ldxios = lrxios ) sshn(:,:) = sshb(:,:) l_1st_euler = .TRUE. ELSE IF( id2 > 0 ) THEN IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : sshb not found in restart files' IF(lwp) write(numout,*) 'set sshb = sshn and force l_1st_euler = true' CALL iom_get( numror, jpdom_data, 'sshn', sshb, ldxios = lrxios ) sshb(:,:) = sshn(:,:) l_1st_euler = .TRUE. ELSE IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : sshb and sshn not found in restart file' IF(lwp) write(numout,*) 'set sshb = sshn = 0 and force l_1st_euler = true' sshb(:,:) = 0._wp sshn(:,:) = 0._wp l_1st_euler = .TRUE. ENDIF ELSE !* Initialize at "rest" ! IF( ll_wd ) THEN ! MJB ll_wd edits start here - these are essential ! IF( cn_cfg == 'wad' ) THEN ! Wetting and drying test case CALL usr_def_istate( gdept_b, tmask, tsb, ub, vb, sshb ) tsn (:,:,:,:) = tsb (:,:,:,:) ! set now values from to before ones sshn (:,:) = sshb(:,:) un (:,:,:) = ub (:,:,:) vn (:,:,:) = vb (:,:,:) ELSE ! Not the test case sshn(:,:) = -ssh_ref sshb(:,:) = -ssh_ref ! DO jj = 1, jpj DO ji = 1, jpi IF( ht_0(ji,jj)-ssh_ref < rn_wdmin1 ) THEN ! if total depth is less than min depth sshb(ji,jj) = rn_wdmin1 - (ht_0(ji,jj) ) sshn(ji,jj) = rn_wdmin1 - (ht_0(ji,jj) ) ssha(ji,jj) = rn_wdmin1 - (ht_0(ji,jj) ) ENDIF END DO END DO ENDIF ! If test case else ! DO ji = 1, jpi DO jj = 1, jpj IF ( ht_0(ji,jj) /= 0._wp .AND. NINT( ssmask(ji,jj) ) == 1 ) THEN CALL ctl_stop( 'dom_vvl_rst: ht_0 must be positive at potentially wet points' ) ENDIF END DO END DO ! ELSE ! ! Just to read set ssh in fact, called latter once vertical grid is set up: ! CALL usr_def_istate( gdept_0, tmask, tsb, ub, vb, sshb ) sshn(:,:) = 0._wp sshb(:,:) = 0._wp ! END IF ! ENDIF ! ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN ! Create restart file ! ! =================== !!gm DO NOTHING, sshb and sshn are written in restart.F90 ENDIF ! END SUBROUTINE dom_vvl_rst SUBROUTINE dom_vvl_ctl !!--------------------------------------------------------------------- !! *** ROUTINE dom_vvl_ctl *** !! !! ** Purpose : Control the consistency between namelist options !! for vertical coordinate !!---------------------------------------------------------------------- INTEGER :: ioptio, ios !! NAMELIST/nam_vvl/ ln_vvl_zstar, ln_vvl_ztilde, ln_vvl_layer, ln_vvl_ztilde_as_zstar, & & ln_vvl_zstar_at_eqtor , rn_ahe3 , rn_rst_e3t , & & rn_lf_cutoff , rn_zdef_max , ln_vvl_dbg ! not yet implemented: ln_vvl_kepe !!---------------------------------------------------------------------- ! REWIND( numnam_ref ) ! Namelist nam_vvl in reference namelist : READ ( numnam_ref, nam_vvl, IOSTAT = ios, ERR = 901) 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_vvl in reference namelist', lwp ) REWIND( numnam_cfg ) ! Namelist nam_vvl in configuration namelist : Parameters of the run READ ( numnam_cfg, nam_vvl, IOSTAT = ios, ERR = 902 ) 902 IF( ios > 0 ) CALL ctl_nam ( ios , 'nam_vvl in configuration namelist', lwp ) IF(lwm) WRITE ( numond, nam_vvl ) ! IF(lwp) THEN ! Namelist print WRITE(numout,*) WRITE(numout,*) 'dom_vvl_ctl : choice/control of the variable vertical coordinate' WRITE(numout,*) '~~~~~~~~~~~' WRITE(numout,*) ' Namelist nam_vvl : chose a vertical coordinate' WRITE(numout,*) ' zstar ln_vvl_zstar = ', ln_vvl_zstar WRITE(numout,*) ' ztilde ln_vvl_ztilde = ', ln_vvl_ztilde WRITE(numout,*) ' layer ln_vvl_layer = ', ln_vvl_layer WRITE(numout,*) ' ztilde as zstar ln_vvl_ztilde_as_zstar = ', ln_vvl_ztilde_as_zstar WRITE(numout,*) ' ztilde near the equator ln_vvl_zstar_at_eqtor = ', ln_vvl_zstar_at_eqtor WRITE(numout,*) ' !' WRITE(numout,*) ' thickness diffusion coefficient rn_ahe3 = ', rn_ahe3 WRITE(numout,*) ' maximum e3t deformation fractional change rn_zdef_max = ', rn_zdef_max IF( ln_vvl_ztilde_as_zstar ) THEN WRITE(numout,*) ' ztilde running in zstar emulation mode (ln_vvl_ztilde_as_zstar=T) ' WRITE(numout,*) ' ignoring namelist timescale parameters and using:' WRITE(numout,*) ' hard-wired : z-tilde to zstar restoration timescale (days)' WRITE(numout,*) ' rn_rst_e3t = 0.e0' WRITE(numout,*) ' hard-wired : z-tilde cutoff frequency of low-pass filter (days)' WRITE(numout,*) ' rn_lf_cutoff = 1/rn_Dt' ELSE WRITE(numout,*) ' z-tilde to zstar restoration timescale (days) rn_rst_e3t = ', rn_rst_e3t WRITE(numout,*) ' z-tilde cutoff frequency of low-pass filter (days) rn_lf_cutoff = ', rn_lf_cutoff ENDIF WRITE(numout,*) ' debug prints flag ln_vvl_dbg = ', ln_vvl_dbg ENDIF ! !!gm IF ( ln_vvl_ztilde .OR. ln_vvl_ztilde_as_zstar ) CALL ctl_stop( 'z-tilde NOT available in this branch' ) !!gm ! ioptio = 0 ! Parameter control IF( ln_vvl_ztilde_as_zstar ) ln_vvl_ztilde = .true. IF( ln_vvl_zstar ) ioptio = ioptio + 1 IF( ln_vvl_ztilde ) ioptio = ioptio + 1 IF( ln_vvl_layer ) ioptio = ioptio + 1 ! IF( ioptio /= 1 ) CALL ctl_stop( 'Choose ONE vertical coordinate in namelist nam_vvl' ) IF( .NOT. ln_vvl_zstar .AND. ln_isf ) CALL ctl_stop( 'Only vvl_zstar has been tested with ice shelf cavity' ) ! IF(lwp) THEN ! Print the choice WRITE(numout,*) IF( ln_vvl_zstar ) WRITE(numout,*) ' ==>>> zstar vertical coordinate is used' IF( ln_vvl_ztilde ) WRITE(numout,*) ' ==>>> ztilde vertical coordinate is used' IF( ln_vvl_layer ) WRITE(numout,*) ' ==>>> layer vertical coordinate is used' IF( ln_vvl_ztilde_as_zstar ) WRITE(numout,*) ' ==>>> to emulate a zstar coordinate' ENDIF ! #if defined key_agrif IF( (.NOT.Agrif_Root()).AND.(.NOT.ln_vvl_zstar) ) CALL ctl_stop( 'AGRIF is implemented with zstar coordinate only' ) #endif ! END SUBROUTINE dom_vvl_ctl SUBROUTINE ssh2e3_now !!---------------------------------------------------------------------- !! *** ROUTINE ssh2e3_now *** !!---------------------------------------------------------------------- INTEGER :: ji, jj, jk REAL(wp), DIMENSION(jpi,jpj) :: zssht_h, zsshu_h, zsshv_h, zsshf_h !!---------------------------------------------------------------------- ! ! !== ssh at u- and v-points ==! ! DO jj = 1, jpjm1 ! start from 1 due to f-point DO ji = 1, jpim1 zsshu_h(ji,jj) = 0.50_wp * ( sshn(ji ,jj) + sshn(ji+1,jj ) ) * ssumask(ji,jj) zsshv_h(ji,jj) = 0.50_wp * ( sshn(ji ,jj) + sshn(ji ,jj+1) ) * ssvmask(ji,jj) zsshf_h(ji,jj) = 0.25_wp * ( sshn(ji ,jj) + sshn(ji ,jj+1) & & + sshn(ji+1,jj) + sshn(ji+1,jj+1) ) * ssfmask(ji,jj) END DO END DO CALL lbc_lnk_multi( zsshu_h(:,:),'U', 1._wp , zsshv_h(:,:),'V', 1._wp , zsshf_h(:,:),'F', 1._wp ) ! ! !== ht, hu and hv == ! (and their inverse) ! ht_n (:,:) = ht_0(:,:) + sshn (:,:) hu_n (:,:) = hu_0(:,:) + zsshu_h(:,:) hv_n (:,:) = hv_0(:,:) + zsshv_h(:,:) r1_hu_n(:,:) = ssumask(:,:) / ( hu_n(:,:) + 1._wp - ssumask(:,:) ) ! ss mask mask due to ISF r1_hv_n(:,:) = ssvmask(:,:) / ( hv_n(:,:) + 1._wp - ssvmask(:,:) ) ! ! !== ssh / h factor at t-, u- ,v- & f-points ==! ! zssht_h(:,:) = sshn (:,:) * r1_ht_0(:,:) zsshu_h(:,:) = zsshu_h(:,:) * r1_hu_0(:,:) zsshv_h(:,:) = zsshv_h(:,:) * r1_hv_0(:,:) zsshf_h(:,:) = zsshf_h(:,:) * r1_hf_0(:,:) ! ! !== e3t, e3w , e3u, e3uw , e3v, e3vw , and e3f ==! ! DO jk = 1, jpkm1 e3t_n(:,:,jk) = e3t_0(:,:,jk) * ( 1._wp + zssht_h(:,:) * tmask(:,:,jk) ) e3w_n(:,:,jk) = e3w_0(:,:,jk) * ( 1._wp + zssht_h(:,:) * tmask(:,:,jk) ) ! e3u_n(:,:,jk) = e3u_0(:,:,jk) * ( 1._wp + zsshu_h(:,:) * umask(:,:,jk) ) e3uw_n(:,:,jk) = e3uw_0(:,:,jk) * ( 1._wp + zsshu_h(:,:) * wumask(:,:,jk) ) ! e3v_n(:,:,jk) = e3v_0(:,:,jk) * ( 1._wp + zsshv_h(:,:) * vmask(:,:,jk) ) e3vw_n(:,:,jk) = e3vw_0(:,:,jk) * ( 1._wp + zsshv_h(:,:) * wvmask(:,:,jk) ) ! e3f_n(:,:,jk) = e3f_0(:,:,jk) * ( 1._wp + zsshf_h(:,:) * fmask(:,:,jk) ) END DO ! ! !== depth of t- and w-points ==! ! zssht_h(:,:) = 1._wp + zssht_h(:,:) ! = 1 + sshn / ht_0 ! IF( ln_isfcav ) THEN ! ISF cavities : ssh scaling not applied over the iceshelf thickness DO jk = 1, jpkm1 gdept_n(:,:,jk) = ( gdept_0(:,:,jk) - risfdep(:,:) ) * zssht_h(:,:) + risfdep(:,:) gdepw_n(:,:,jk) = ( gdepw_0(:,:,jk) - risfdep(:,:) ) * zssht_h(:,:) + risfdep(:,:) gde3w_n(:,:,jk) = gdept_n(:,:,jk) - sshn(:,:) END DO ELSE ! no ISF cavities DO jk = 1, jpkm1 gdept_n(:,:,jk) = gdept_0(:,:,jk) * zssht_h(:,:) gdepw_n(:,:,jk) = gdepw_0(:,:,jk) * zssht_h(:,:) gde3w_n(:,:,jk) = gdept_n(:,:,jk) - sshn(:,:) END DO ENDIF ! END SUBROUTINE ssh2e3_now SUBROUTINE ssh2e3_before !!---------------------------------------------------------------------- !! *** ROUTINE ssh2e3_before *** !!---------------------------------------------------------------------- INTEGER :: ji, jj, jk REAL(wp), DIMENSION(jpi,jpj) :: zssht_h, zsshu_h, zsshv_h !!---------------------------------------------------------------------- ! ! !== ssh at u- and v-points ==! DO jj = 2, jpjm1 DO ji = 2, jpim1 zsshu_h(ji,jj) = 0.5_wp * ( sshb(ji ,jj) + sshb(ji+1,jj ) ) * ssumask(ji,jj) zsshv_h(ji,jj) = 0.5_wp * ( sshb(ji ,jj) + sshb(ji ,jj+1) ) * ssvmask(ji,jj) END DO END DO CALL lbc_lnk_multi( zsshu_h(:,:),'U', 1._wp , zsshv_h(:,:),'V', 1._wp ) ! ! !== ht, hu and hv == ! (and their inverse) hu_b(:,:) = hu_0(:,:) + zsshu_h(:,:) hv_b(:,:) = hv_0(:,:) + zsshv_h(:,:) r1_hu_b(:,:) = ssumask(:,:) / ( hu_b(:,:) + 1._wp - ssumask(:,:) ) r1_hv_b(:,:) = ssvmask(:,:) / ( hv_b(:,:) + 1._wp - ssvmask(:,:) ) ! ! ! !== ssh / h factor at t-, u- ,v- & f-points ==! zssht_h(:,:) = sshb (:,:) * r1_ht_0(:,:) zsshu_h (:,:) = zsshu_h(:,:) * r1_hu_0(:,:) zsshv_h (:,:) = zsshv_h(:,:) * r1_hv_0(:,:) ! ! !== e3t, e3w , e3u, e3uw , and e3v, e3vw ==! DO jk = 1, jpkm1 e3t_b(:,:,jk) = e3t_0(:,:,jk) * ( 1._wp + zssht_h(:,:) * tmask(:,:,jk) ) e3w_b(:,:,jk) = e3w_0(:,:,jk) * ( 1._wp + zssht_h(:,:) * tmask(:,:,jk) ) ! e3u_b(:,:,jk) = e3u_0(:,:,jk) * ( 1._wp + zsshu_h (:,:) * umask(:,:,jk) ) e3uw_b(:,:,jk) = e3uw_0(:,:,jk) * ( 1._wp + zsshu_h (:,:) * wumask(:,:,jk) ) ! e3v_b(:,:,jk) = e3v_0(:,:,jk) * ( 1._wp + zsshv_h (:,:) * vmask(:,:,jk) ) e3vw_b(:,:,jk) = e3vw_0(:,:,jk) * ( 1._wp + zsshv_h (:,:) * wvmask(:,:,jk) ) END DO ! END SUBROUTINE ssh2e3_before !!====================================================================== END MODULE domvvl