From 388075a8edd315e144c3203d9fd809df8aeac0eb Mon Sep 17 00:00:00 2001 From: Robinson Hordoir Date: Fri, 5 Sep 2025 12:30:45 +0200 Subject: [PATCH] Commit GMD Hordoir et al. 2025 routines --- src/OCE/TRA/dynhpg.F90_dev | 1286 ++++++++++++++++++++++++++++++++++++ src/OCE/TRA/trabbl.F90_dev | 654 ++++++++++++++++++ 2 files changed, 1940 insertions(+) create mode 100755 src/OCE/TRA/dynhpg.F90_dev create mode 100755 src/OCE/TRA/trabbl.F90_dev diff --git a/src/OCE/TRA/dynhpg.F90_dev b/src/OCE/TRA/dynhpg.F90_dev new file mode 100755 index 0000000..d63de79 --- /dev/null +++ b/src/OCE/TRA/dynhpg.F90_dev @@ -0,0 +1,1286 @@ +MODULE dynhpg + !!====================================================================== + !! *** MODULE dynhpg *** + !! Ocean dynamics: hydrostatic pressure gradient trend + !!====================================================================== + !! History : OPA ! 1987-09 (P. Andrich, M.-A. Foujols) hpg_zco: Original code + !! 5.0 ! 1991-11 (G. Madec) + !! 7.0 ! 1996-01 (G. Madec) hpg_sco: Original code for s-coordinates + !! 8.0 ! 1997-05 (G. Madec) split dynber into dynkeg and dynhpg + !! 8.5 ! 2002-07 (G. Madec) F90: Free form and module + !! 8.5 ! 2002-08 (A. Bozec) hpg_zps: Original code + !! NEMO 1.0 ! 2005-10 (A. Beckmann, B.W. An) various s-coordinate options + !! ! Original code for hpg_ctl, hpg_hel hpg_wdj, hpg_djc, hpg_rot + !! - ! 2005-11 (G. Madec) style & small optimisation + !! 3.3 ! 2010-10 (C. Ethe, G. Madec) reorganisation of initialisation phase + !! 3.4 ! 2011-11 (H. Liu) hpg_prj: Original code for s-coordinates + !! ! (A. Coward) suppression of hel, wdj and rot options + !! 3.6 ! 2014-11 (P. Mathiot) hpg_isf: original code for ice shelf cavity + !! 4.2 ! 2020-12 (M. Bell, A. Young) hpg_djc: revised djc scheme + !!---------------------------------------------------------------------- + + !!---------------------------------------------------------------------- + !! dyn_hpg : update the momentum trend with the now horizontal + !! gradient of the hydrostatic pressure + !! dyn_hpg_init : initialisation and control of options + !! hpg_zco : z-coordinate scheme + !! hpg_zps : z-coordinate plus partial steps (interpolation) + !! hpg_sco : s-coordinate (standard jacobian formulation) + !! hpg_isf : s-coordinate (sco formulation) adapted to ice shelf + !! hpg_djc : s-coordinate (Density Jacobian with Cubic polynomial) + !! hpg_prj : s-coordinate (Pressure Jacobian with Cubic polynomial) + !!---------------------------------------------------------------------- + USE oce ! ocean dynamics and tracers + USE isf_oce , ONLY : risfload ! ice shelf (risfload variable) + USE isfload , ONLY : isf_load ! ice shelf (isf_load routine ) + USE sbc_oce ! surface variable (only for the flag with ice shelf) + USE dom_oce ! ocean space and time domain + USE wet_dry ! wetting and drying + USE phycst ! physical constants + USE trd_oce ! trends: ocean variables + USE trddyn ! trend manager: dynamics + USE zpshde ! partial step: hor. derivative (zps_hde routine) + ! + USE in_out_manager ! I/O manager + USE prtctl ! Print control + USE lbclnk ! lateral boundary condition + USE lib_mpp ! MPP library + USE eosbn2 ! compute density + USE timing ! Timing + USE iom +!!Robinson + USE trabbl +!!Robinson + IMPLICIT NONE + PRIVATE + + PUBLIC dyn_hpg ! routine called by step module + PUBLIC dyn_hpg_init ! routine called by opa module + + ! !!* Namelist namdyn_hpg : hydrostatic pressure gradient + LOGICAL, PUBLIC :: ln_hpg_zco !: z-coordinate - full steps + LOGICAL, PUBLIC :: ln_hpg_zps !: z-coordinate - partial steps (interpolation) + LOGICAL, PUBLIC :: ln_hpg_sco !: s-coordinate (standard jacobian formulation) + LOGICAL, PUBLIC :: ln_hpg_djc !: s-coordinate (Density Jacobian with Cubic polynomial) + LOGICAL, PUBLIC :: ln_hpg_prj !: s-coordinate (Pressure Jacobian scheme) + LOGICAL, PUBLIC :: ln_hpg_isf !: s-coordinate similar to sco modify for isf + + ! !! Flag to control the type of hydrostatic pressure gradient + INTEGER, PARAMETER :: np_ERROR =-10 ! error in specification of lateral diffusion + INTEGER, PARAMETER :: np_zco = 0 ! z-coordinate - full steps + INTEGER, PARAMETER :: np_zps = 1 ! z-coordinate - partial steps (interpolation) + INTEGER, PARAMETER :: np_sco = 2 ! s-coordinate (standard jacobian formulation) + INTEGER, PARAMETER :: np_djc = 3 ! s-coordinate (Density Jacobian with Cubic polynomial) + INTEGER, PARAMETER :: np_prj = 4 ! s-coordinate (Pressure Jacobian scheme) + INTEGER, PARAMETER :: np_isf = 5 ! s-coordinate similar to sco modify for isf + ! + INTEGER, PUBLIC :: nhpg !: type of pressure gradient scheme used ! (deduced from ln_hpg_... flags) (PUBLIC for TAM) + ! + LOGICAL :: ln_hpg_djc_vnh, ln_hpg_djc_vnv ! flag to specify hpg_djc boundary condition type + REAL(wp), PUBLIC :: aco_bc_hor, bco_bc_hor, aco_bc_vrt, bco_bc_vrt !: coefficients for hpg_djc hor and vert boundary conditions + + !! * Substitutions +# include "do_loop_substitute.h90" +# include "domzgr_substitute.h90" + + !!---------------------------------------------------------------------- + !! NEMO/OCE 4.0 , NEMO Consortium (2018) + !! $Id: dynhpg.F90 15529 2021-11-23 15:00:19Z techene $ + !! Software governed by the CeCILL license (see ./LICENSE) + !!---------------------------------------------------------------------- +CONTAINS + + SUBROUTINE dyn_hpg( kt, Kmm, puu, pvv, Krhs ) + !!--------------------------------------------------------------------- + !! *** ROUTINE dyn_hpg *** + !! + !! ** Method : Call the hydrostatic pressure gradient routine + !! using the scheme defined in the namelist + !! + !! ** Action : - Update (puu(:,:,:,Krhs),pvv(:,:,:,Krhs)) with the now hydrastatic pressure trend + !! - send trends to trd_dyn for futher diagnostics (l_trddyn=T) + !!---------------------------------------------------------------------- + INTEGER , INTENT( in ) :: kt ! ocean time-step index + INTEGER , INTENT( in ) :: Kmm, Krhs ! ocean time level indices + REAL(dp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) :: puu, pvv ! ocean velocities and RHS of momentum equation + ! + REAL(dp), ALLOCATABLE, DIMENSION(:,:,:) :: ztrdu, ztrdv + !!---------------------------------------------------------------------- + ! + IF( ln_timing ) CALL timing_start('dyn_hpg') + ! + IF( l_trddyn ) THEN ! Temporary saving of puu(:,:,:,Krhs) and pvv(:,:,:,Krhs) trends (l_trddyn) + ALLOCATE( ztrdu(jpi,jpj,jpk) , ztrdv(jpi,jpj,jpk) ) + ztrdu(:,:,:) = puu(:,:,:,Krhs) + ztrdv(:,:,:) = pvv(:,:,:,Krhs) + ENDIF + ! + SELECT CASE ( nhpg ) ! Hydrostatic pressure gradient computation + CASE ( np_zco ) ; CALL hpg_zco ( kt, Kmm, puu, pvv, Krhs ) ! z-coordinate + CASE ( np_zps ) ; CALL hpg_zps ( kt, Kmm, puu, pvv, Krhs ) ! z-coordinate plus partial steps (interpolation) + CASE ( np_sco ) ; CALL hpg_sco ( kt, Kmm, puu, pvv, Krhs ) ! s-coordinate (standard jacobian formulation) + CASE ( np_djc ) ; CALL hpg_djc ( kt, Kmm, puu, pvv, Krhs ) ! s-coordinate (Density Jacobian with Cubic polynomial) + CASE ( np_prj ) ; CALL hpg_prj ( kt, Kmm, puu, pvv, Krhs ) ! s-coordinate (Pressure Jacobian scheme) + CASE ( np_isf ) ; CALL hpg_isf ( kt, Kmm, puu, pvv, Krhs ) ! s-coordinate similar to sco modify for ice shelf + END SELECT + ! + IF( l_trddyn ) THEN ! save the hydrostatic pressure gradient trends for momentum trend diagnostics + ztrdu(:,:,:) = puu(:,:,:,Krhs) - ztrdu(:,:,:) + ztrdv(:,:,:) = pvv(:,:,:,Krhs) - ztrdv(:,:,:) + CALL trd_dyn( ztrdu, ztrdv, jpdyn_hpg, kt, Kmm ) + DEALLOCATE( ztrdu , ztrdv ) + ENDIF + ! + IF(sn_cfctl%l_prtctl) CALL prt_ctl( tab3d_1=puu(:,:,:,Krhs), clinfo1=' hpg - Ua: ', mask1=umask, & + & tab3d_2=pvv(:,:,:,Krhs), clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' ) + ! + IF( ln_timing ) CALL timing_stop('dyn_hpg') + ! + END SUBROUTINE dyn_hpg + + + SUBROUTINE dyn_hpg_init( Kmm ) + !!---------------------------------------------------------------------- + !! *** ROUTINE dyn_hpg_init *** + !! + !! ** Purpose : initializations for the hydrostatic pressure gradient + !! computation and consistency control + !! + !! ** Action : Read the namelist namdyn_hpg and check the consistency + !! with the type of vertical coordinate used (zco, zps, sco) + !!---------------------------------------------------------------------- + INTEGER, INTENT( in ) :: Kmm ! ocean time level index + ! + INTEGER :: ioptio = 0 ! temporary integer + INTEGER :: ios ! Local integer output status for namelist read + !! + INTEGER :: ji, jj, jk, ikt ! dummy loop indices ISF + REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: zts_top, zrhd ! hypothesys on isf density + REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zrhdtop_isf ! density at bottom of ISF + REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: ziceload ! density at bottom of ISF + !! + NAMELIST/namdyn_hpg/ ln_hpg_zco, ln_hpg_zps, ln_hpg_sco, & + & ln_hpg_djc, ln_hpg_prj, ln_hpg_isf, & + & ln_hpg_djc_vnh, ln_hpg_djc_vnv + !!---------------------------------------------------------------------- + ! + READ ( numnam_ref, namdyn_hpg, IOSTAT = ios, ERR = 901) +901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdyn_hpg in reference namelist' ) + ! + READ ( numnam_cfg, namdyn_hpg, IOSTAT = ios, ERR = 902 ) +902 IF( ios > 0 ) CALL ctl_nam ( ios , 'namdyn_hpg in configuration namelist' ) + IF(lwm) WRITE ( numond, namdyn_hpg ) + ! + IF(lwp) THEN ! Control print + WRITE(numout,*) + WRITE(numout,*) 'dyn_hpg_init : hydrostatic pressure gradient initialisation' + WRITE(numout,*) '~~~~~~~~~~~~' + WRITE(numout,*) ' Namelist namdyn_hpg : choice of hpg scheme' + WRITE(numout,*) ' z-coord. - full steps ln_hpg_zco = ', ln_hpg_zco + WRITE(numout,*) ' z-coord. - partial steps (interpolation) ln_hpg_zps = ', ln_hpg_zps + WRITE(numout,*) ' s-coord. (standard jacobian formulation) ln_hpg_sco = ', ln_hpg_sco + WRITE(numout,*) ' s-coord. (standard jacobian formulation) for isf ln_hpg_isf = ', ln_hpg_isf + WRITE(numout,*) ' s-coord. (Density Jacobian: Cubic polynomial) ln_hpg_djc = ', ln_hpg_djc + WRITE(numout,*) ' s-coord. (Pressure Jacobian: Cubic polynomial) ln_hpg_prj = ', ln_hpg_prj + ENDIF + ! + IF( .NOT.ln_linssh .AND. (ln_hpg_zco.OR.ln_hpg_zps) ) & + & CALL ctl_stop( 'dyn_hpg_init : non-linear free surface incompatible with hpg_zco or hpg_zps' ) + ! + IF( (.NOT.ln_hpg_isf .AND. ln_isfcav) .OR. (ln_hpg_isf .AND. .NOT.ln_isfcav) ) & + & CALL ctl_stop( 'dyn_hpg_init : ln_hpg_isf=T requires ln_isfcav=T and vice versa' ) + ! + ! + ! ! Set nhpg from ln_hpg_... flags & consistency check + nhpg = np_ERROR + ioptio = 0 + IF( ln_hpg_zco ) THEN ; nhpg = np_zco ; ioptio = ioptio +1 ; ENDIF + IF( ln_hpg_zps ) THEN ; nhpg = np_zps ; ioptio = ioptio +1 ; ENDIF + IF( ln_hpg_sco ) THEN ; nhpg = np_sco ; ioptio = ioptio +1 ; ENDIF + IF( ln_hpg_djc ) THEN ; nhpg = np_djc ; ioptio = ioptio +1 ; ENDIF + IF( ln_hpg_prj ) THEN ; nhpg = np_prj ; ioptio = ioptio +1 ; ENDIF + IF( ln_hpg_isf ) THEN ; nhpg = np_isf ; ioptio = ioptio +1 ; ENDIF + ! + IF( ioptio /= 1 ) CALL ctl_stop( 'NO or several hydrostatic pressure gradient options used' ) + ! + IF(lwp) THEN + WRITE(numout,*) + SELECT CASE( nhpg ) + CASE( np_zco ) ; WRITE(numout,*) ' ==>>> z-coord. - full steps ' + CASE( np_zps ) ; WRITE(numout,*) ' ==>>> z-coord. - partial steps (interpolation)' + CASE( np_sco ) ; WRITE(numout,*) ' ==>>> s-coord. (standard jacobian formulation)' + CASE( np_djc ) ; WRITE(numout,*) ' ==>>> s-coord. (Density Jacobian: Cubic polynomial)' + CASE( np_prj ) ; WRITE(numout,*) ' ==>>> s-coord. (Pressure Jacobian: Cubic polynomial)' + CASE( np_isf ) ; WRITE(numout,*) ' ==>>> s-coord. (standard jacobian formulation) for isf' + END SELECT + WRITE(numout,*) + ENDIF + ! + IF ( ln_hpg_djc ) THEN + IF (ln_hpg_djc_vnh) THEN ! Von Neumann boundary condition + IF(lwp) WRITE(numout,*) ' horizontal bc: von Neumann ' + aco_bc_hor = 6.0_wp/5.0_wp + bco_bc_hor = 7.0_wp/15.0_wp + ELSE ! Linear extrapolation + IF(lwp) WRITE(numout,*) ' horizontal bc: linear extrapolation' + aco_bc_hor = 3.0_wp/2.0_wp + bco_bc_hor = 1.0_wp/2.0_wp + END IF + IF (ln_hpg_djc_vnv) THEN ! Von Neumann boundary condition + IF(lwp) WRITE(numout,*) ' vertical bc: von Neumann ' + aco_bc_vrt = 6.0_wp/5.0_wp + bco_bc_vrt = 7.0_wp/15.0_wp + ELSE ! Linear extrapolation + IF(lwp) WRITE(numout,*) ' vertical bc: linear extrapolation' + aco_bc_vrt = 3.0_wp/2.0_wp + bco_bc_vrt = 1.0_wp/2.0_wp + END IF + END IF + ! + END SUBROUTINE dyn_hpg_init + + + SUBROUTINE hpg_zco( kt, Kmm, puu, pvv, Krhs ) + !!--------------------------------------------------------------------- + !! *** ROUTINE hpg_zco *** + !! + !! ** Method : z-coordinate case, levels are horizontal surfaces. + !! The now hydrostatic pressure gradient at a given level, jk, + !! is computed by taking the vertical integral of the in-situ + !! density gradient along the model level from the suface to that + !! level: zhpi = grav ..... + !! zhpj = grav ..... + !! add it to the general momentum trend (puu(:,:,:,Krhs),pvv(:,:,:,Krhs)). + !! puu(:,:,:,Krhs) = puu(:,:,:,Krhs) - 1/e1u * zhpi + !! pvv(:,:,:,Krhs) = pvv(:,:,:,Krhs) - 1/e2v * zhpj + !! + !! ** Action : - Update (puu(:,:,:,Krhs),pvv(:,:,:,Krhs)) with the now hydrastatic pressure trend + !!---------------------------------------------------------------------- + INTEGER , INTENT( in ) :: kt ! ocean time-step index + INTEGER , INTENT( in ) :: Kmm, Krhs ! ocean time level indices + 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) :: zcoef0, zcoef1 ! temporary scalars + REAL(wp), DIMENSION(A2D(nn_hls)) :: zhpi, zhpj + !!---------------------------------------------------------------------- + ! + IF( .NOT. l_istiled .OR. ntile == 1 ) THEN ! Do only on the first tile + IF( kt == nit000 ) THEN + IF(lwp) WRITE(numout,*) + IF(lwp) WRITE(numout,*) 'dyn:hpg_zco : hydrostatic pressure gradient trend' + IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ z-coordinate case ' + ENDIF + ENDIF + ! + zcoef0 = - grav * 0.5_wp ! Local constant initialization + ! + DO_2D( 0, 0, 0, 0 ) ! Surface value + zcoef1 = zcoef0 * e3w(ji,jj,1,Kmm) + ! ! hydrostatic pressure gradient + zhpi(ji,jj) = zcoef1 * ( rhd(ji+1,jj,1) - rhd(ji,jj,1) ) * r1_e1u(ji,jj) + zhpj(ji,jj) = zcoef1 * ( rhd(ji,jj+1,1) - rhd(ji,jj,1) ) * r1_e2v(ji,jj) + ! ! add to the general momentum trend + puu(ji,jj,1,Krhs) = puu(ji,jj,1,Krhs) + zhpi(ji,jj) + pvv(ji,jj,1,Krhs) = pvv(ji,jj,1,Krhs) + zhpj(ji,jj) + END_2D + ! + DO_3D( 0, 0, 0, 0, 2, jpkm1 ) ! interior value (2= 1 ) THEN ! on i-direction (level 2 or more) + puu (ji,jj,iku,Krhs) = puu(ji,jj,iku,Krhs) - zhpi(ji,jj,iku) ! subtract old value + zhpi(ji,jj,iku) = zhpi(ji,jj,iku-1) & ! compute the new one + & + zcoef2 * ( rhd(ji+1,jj,iku-1) - rhd(ji,jj,iku-1) + zgru(ji,jj) ) * r1_e1u(ji,jj) + puu (ji,jj,iku,Krhs) = puu(ji,jj,iku,Krhs) + zhpi(ji,jj,iku) ! add the new one to the general momentum trend + ENDIF + IF( ikv > 1 ) THEN ! on j-direction (level 2 or more) + pvv (ji,jj,ikv,Krhs) = pvv(ji,jj,ikv,Krhs) - zhpj(ji,jj,ikv) ! subtract old value + zhpj(ji,jj,ikv) = zhpj(ji,jj,ikv-1) & ! compute the new one + & + zcoef3 * ( rhd(ji,jj+1,ikv-1) - rhd(ji,jj,ikv-1) + zgrv(ji,jj) ) * r1_e2v(ji,jj) + pvv (ji,jj,ikv,Krhs) = pvv(ji,jj,ikv,Krhs) + zhpj(ji,jj,ikv) ! add the new one to the general momentum trend + ENDIF + END_2D + ! + END SUBROUTINE hpg_zps + + + SUBROUTINE hpg_sco( kt, Kmm, puu, pvv, Krhs ) + !!--------------------------------------------------------------------- + !! *** ROUTINE hpg_sco *** + !! + !! ** Method : s-coordinate case. Jacobian scheme. + !! The now hydrostatic pressure gradient at a given level, jk, + !! is computed by taking the vertical integral of the in-situ + !! density gradient along the model level from the suface to that + !! level. s-coordinates (ln_sco): a corrective term is added + !! to the horizontal pressure gradient : + !! zhpi = grav ..... + 1/e1u mi(rhd) di[ grav dep3w ] + !! zhpj = grav ..... + 1/e2v mj(rhd) dj[ grav dep3w ] + !! add it to the general momentum trend (puu(:,:,:,Krhs),pvv(:,:,:,Krhs)). + !! puu(:,:,:,Krhs) = puu(:,:,:,Krhs) - 1/e1u * zhpi + !! pvv(:,:,:,Krhs) = pvv(:,:,:,Krhs) - 1/e2v * zhpj + !! + !! ** Action : - Update (puu(:,:,:,Krhs),pvv(:,:,:,Krhs)) with the now hydrastatic pressure trend + !!---------------------------------------------------------------------- + INTEGER , INTENT( in ) :: kt ! ocean time-step index + INTEGER , INTENT( in ) :: Kmm, Krhs ! ocean time level indices + REAL(dp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) :: puu, pvv ! ocean velocities and RHS of momentum equation + !! + INTEGER :: ji, jj, jk, jii, jjj ! dummy loop indices + REAL(wp) :: zcoef0, zuap, zvap, ztmp ! local scalars + REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: zhpi, zhpj + !!---------------------------------------------------------------------- + ! + ! + IF( .NOT. l_istiled .OR. ntile == 1 ) THEN ! Do only on the first tile + IF( kt == nit000 ) THEN + IF(lwp) WRITE(numout,*) + IF(lwp) WRITE(numout,*) 'dyn:hpg_sco : hydrostatic pressure gradient trend' + IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ s-coordinate case, OCE original scheme used' + ENDIF + ENDIF + ! + zcoef0 = - grav * 0.5_wp + ! + ! + DO_2D( 0, 0, 0, 0 ) ! Surface value + ! ! hydrostatic pressure gradient along s-surfaces + zhpi(ji,jj,1) = zcoef0 * r1_e1u(ji,jj) & + & * ( e3w(ji+1,jj ,1,Kmm) * rhd(ji+1,jj ,1) & + & - e3w(ji ,jj ,1,Kmm) * rhd(ji ,jj ,1) ) + zhpj(ji,jj,1) = zcoef0 * r1_e2v(ji,jj) & + & * ( e3w(ji ,jj+1,1,Kmm) * rhd(ji ,jj+1,1) & + & - e3w(ji ,jj ,1,Kmm) * rhd(ji ,jj ,1) ) + ! ! s-coordinate pressure gradient correction + zuap = -zcoef0 * ( rhd (ji+1,jj,1) + rhd (ji,jj,1) ) & + & * ( gde3w(ji+1,jj,1) - gde3w(ji,jj,1) ) * r1_e1u(ji,jj) + zvap = -zcoef0 * ( rhd (ji,jj+1,1) + rhd (ji,jj,1) ) & + & * ( gde3w(ji,jj+1,1) - gde3w(ji,jj,1) ) * r1_e2v(ji,jj) + ! + ! ! add to the general momentum trend + puu(ji,jj,1,Krhs) = puu(ji,jj,1,Krhs) + zhpi(ji,jj,1) + zuap + pvv(ji,jj,1,Krhs) = pvv(ji,jj,1,Krhs) + zhpj(ji,jj,1) + zvap + END_2D + ! + DO_3D( 0, 0, 0, 0, 2, jpkm1 ) ! interior value (2=1 ) THEN + iktb = mbkt(ji,jj) + zdrho_k(ji,jj,iktb) = aco_bc_vrt * ( rhd(ji,jj,iktb) - rhd(ji,jj,iktb-1) ) - bco_bc_vrt * zdrho_k(ji,jj,iktb-1) + zdz_k (ji,jj,iktb) = aco_bc_vrt * (-gde3w(ji,jj,iktb) + gde3w(ji,jj,iktb-1) ) - bco_bc_vrt * zdz_k (ji,jj,iktb-1) + END IF + END_2D + + !-------------------------------------------------------------- + ! 4. Compute side face integrals + !------------------------------------------------------------- + +!! ssh replaces e3w_n ; gde3w is a depth; the formulae involve heights +!! rho_k stores grav * FX / rho_0 + + !-------------------------------------------------------------- + ! 4. a) Upper half of top-most grid box, compute and store + !------------------------------------------------------------- +! *** AY note: ssh(ji,jj,Kmm) + gde3w(ji,jj,1) = e3w(ji,jj,1,Kmm) + DO_2D( 0, 1, 0, 1) + z_rho_k(ji,jj,1) = grav * gdept(ji,jj,1,Kmm) & + & * ( rhd(ji,jj,1) & + & -0.5_wp * ( rhd(ji,jj,2) - rhd(ji,jj,1) ) & + & * gdept(ji,jj,1,Kmm) / e3w(ji,jj,2,Kmm) & + & ) + END_2D + + !-------------------------------------------------------------- + ! 4. b) Interior faces, compute and store + !------------------------------------------------------------- + + DO_3D( 0, 1, 0, 1, 2, jpkm1 ) + z_rho_k(ji,jj,jk) = zcoef0 * ( rhd (ji,jj,jk) + rhd (ji,jj,jk-1) ) & + & * ( - gde3w(ji,jj,jk) + gde3w(ji,jj,jk-1) ) & + & + z_grav_10 * ( & + & ( zdrho_k (ji,jj,jk) - zdrho_k (ji,jj,jk-1) ) & + & * ( - gde3w(ji,jj,jk) + gde3w(ji,jj,jk-1) - z1_12 * ( zdz_k (ji,jj,jk) + zdz_k (ji,jj,jk-1) ) ) & + & - ( zdz_k (ji,jj,jk) - zdz_k (ji,jj,jk-1) ) & + & * ( rhd (ji,jj,jk) - rhd (ji,jj,jk-1) - z1_12 * ( zdrho_k(ji,jj,jk) + zdrho_k(ji,jj,jk-1) ) ) & + & ) + END_3D + + !---------------------------------------------------------------------------------------- + ! 5. compute and store elementary horizontal differences in provisional arrays + !---------------------------------------------------------------------------------------- + zdrhox(:,:,:) = 0._wp + zdzx (:,:,:) = 0._wp + zdrhoy(:,:,:) = 0._wp + zdzy (:,:,:) = 0._wp + + DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpkm1 ) + zdrhox(ji,jj,jk) = rhd (ji+1,jj ,jk) - rhd (ji ,jj ,jk) + zdzx (ji,jj,jk) = gde3w(ji ,jj ,jk) - gde3w(ji+1,jj ,jk) + zdrhoy(ji,jj,jk) = rhd (ji ,jj+1,jk) - rhd (ji ,jj ,jk) + zdzy (ji,jj,jk) = gde3w(ji ,jj ,jk) - gde3w(ji ,jj+1,jk) + END_3D + + IF( nn_hls == 1 ) CALL lbc_lnk( 'dynhpg', zdrhox, 'U', -1._wp, zdzx, 'U', -1._wp, zdrhoy, 'V', -1._wp, zdzy, 'V', -1._wp ) + + !------------------------------------------------------------------------- + ! 6. compute harmonic averages using eq. 5.18 + !------------------------------------------------------------------------- + + DO_3D( 0, 1, 0, 1, 1, jpkm1 ) + cffu = MAX( 2._wp * zdrhox(ji-1,jj,jk) * zdrhox(ji,jj,jk), 0._wp ) + z1_cff = zdrhox(ji-1,jj,jk) + zdrhox(ji,jj,jk) + zdrho_i(ji,jj,jk) = cffu / SIGN( MAX( ABS(z1_cff), zep ), z1_cff ) + + cffx = MAX( 2._wp * zdzx(ji-1,jj,jk) * zdzx(ji,jj,jk), 0._wp ) + z1_cff = zdzx(ji-1,jj,jk) + zdzx(ji,jj,jk) + zdz_i(ji,jj,jk) = cffx / SIGN( MAX( ABS(z1_cff), zep ), z1_cff ) + + cffv = MAX( 2._wp * zdrhoy(ji,jj-1,jk) * zdrhoy(ji,jj,jk), 0._wp ) + z1_cff = zdrhoy(ji,jj-1,jk) + zdrhoy(ji,jj,jk) + zdrho_j(ji,jj,jk) = cffv / SIGN( MAX( ABS(z1_cff), zep ), z1_cff ) + + cffy = MAX( 2._wp * zdzy(ji,jj-1,jk) * zdzy(ji,jj,jk), 0._wp ) + z1_cff = zdzy(ji,jj-1,jk) + zdzy(ji,jj,jk) + zdz_j(ji,jj,jk) = cffy / SIGN( MAX( ABS(z1_cff), zep ), z1_cff ) + END_3D + +!!! Note that zdzx, zdzy, zdzz, zdrhox, zdrhoy and zdrhoz should NOT be used beyond this point + + !---------------------------------------------------------------------------------- + ! 6B. apply boundary conditions at side boundaries using 5.36-5.37 + !---------------------------------------------------------------------------------- + + DO jk = 1, jpkm1 + zz_drho_i(:,:) = zdrho_i(:,:,jk) + zz_dz_i (:,:) = zdz_i (:,:,jk) + zz_drho_j(:,:) = zdrho_j(:,:,jk) + zz_dz_j (:,:) = zdz_j (:,:,jk) + ! Walls coming from left: should check from 2 to jpi-1 (and jpj=2-jpj) + DO_2D( 0, 0, 0, 1 ) + IF ( umask(ji,jj,jk) > 0.5_wp .AND. umask(ji-1,jj,jk) < 0.5_wp .AND. umask(ji+1,jj,jk) > 0.5_wp) THEN + zz_drho_i(ji,jj) = aco_bc_hor * ( rhd (ji+1,jj,jk) - rhd (ji,jj,jk) ) - bco_bc_hor * zdrho_i(ji+1,jj,jk) + zz_dz_i (ji,jj) = aco_bc_hor * (-gde3w(ji+1,jj,jk) + gde3w(ji,jj,jk) ) - bco_bc_hor * zdz_i (ji+1,jj,jk) + END IF + END_2D + ! Walls coming from right: should check from 3 to jpi (and jpj=2-jpj) + DO_2D( -1, 1, 0, 1 ) + IF ( umask(ji,jj,jk) < 0.5_wp .AND. umask(ji-1,jj,jk) > 0.5_wp .AND. umask(ji-2,jj,jk) > 0.5_wp) THEN + zz_drho_i(ji,jj) = aco_bc_hor * ( rhd (ji,jj,jk) - rhd (ji-1,jj,jk) ) - bco_bc_hor * zdrho_i(ji-1,jj,jk) + zz_dz_i (ji,jj) = aco_bc_hor * (-gde3w(ji,jj,jk) + gde3w(ji-1,jj,jk) ) - bco_bc_hor * zdz_i (ji-1,jj,jk) + END IF + END_2D + ! Walls coming from left: should check from 2 to jpj-1 (and jpi=2-jpi) + DO_2D( 0, 1, 0, 0 ) + IF ( vmask(ji,jj,jk) > 0.5_wp .AND. vmask(ji,jj-1,jk) < 0.5_wp .AND. vmask(ji,jj+1,jk) > 0.5_wp) THEN + zz_drho_j(ji,jj) = aco_bc_hor * ( rhd (ji,jj+1,jk) - rhd (ji,jj,jk) ) - bco_bc_hor * zdrho_j(ji,jj+1,jk) + zz_dz_j (ji,jj) = aco_bc_hor * (-gde3w(ji,jj+1,jk) + gde3w(ji,jj,jk) ) - bco_bc_hor * zdz_j (ji,jj+1,jk) + END IF + END_2D + ! Walls coming from right: should check from 3 to jpj (and jpi=2-jpi) + DO_2D( 0, 1, -1, 1 ) + IF ( vmask(ji,jj,jk) < 0.5_wp .AND. vmask(ji,jj-1,jk) > 0.5_wp .AND. vmask(ji,jj-2,jk) > 0.5_wp) THEN + zz_drho_j(ji,jj) = aco_bc_hor * ( rhd (ji,jj,jk) - rhd (ji,jj-1,jk) ) - bco_bc_hor * zdrho_j(ji,jj-1,jk) + zz_dz_j (ji,jj) = aco_bc_hor * (-gde3w(ji,jj,jk) + gde3w(ji,jj-1,jk) ) - bco_bc_hor * zdz_j (ji,jj-1,jk) + END IF + END_2D + zdrho_i(:,:,jk) = zz_drho_i(:,:) + zdz_i (:,:,jk) = zz_dz_i (:,:) + zdrho_j(:,:,jk) = zz_drho_j(:,:) + zdz_j (:,:,jk) = zz_dz_j (:,:) + END DO + + !-------------------------------------------------------------- + ! 7. Calculate integrals on side faces + !------------------------------------------------------------- + + DO_3D( 0, 0, 0, 0, 1, jpkm1 ) +! two -ve signs cancel in next two lines (within zcoef0 and because gde3w is a depth not a height) + z_rho_i(ji,jj,jk) = zcoef0 * ( rhd (ji+1,jj,jk) + rhd (ji,jj,jk) ) & + & * ( gde3w(ji+1,jj,jk) - gde3w(ji,jj,jk) ) + IF ( umask(ji-1, jj, jk) > 0.5 .OR. umask(ji+1, jj, jk) > 0.5 ) THEN + z_rho_i(ji,jj,jk) = z_rho_i(ji,jj,jk) - z_grav_10 * ( & + & ( zdrho_i (ji+1,jj,jk) - zdrho_i (ji,jj,jk) ) & + & * ( - gde3w(ji+1,jj,jk) + gde3w(ji,jj,jk) - z1_12 * ( zdz_i (ji+1,jj,jk) + zdz_i (ji,jj,jk) ) ) & + & - ( zdz_i (ji+1,jj,jk) - zdz_i (ji,jj,jk) ) & + & * ( rhd (ji+1,jj,jk) - rhd (ji,jj,jk) - z1_12 * ( zdrho_i(ji+1,jj,jk) + zdrho_i(ji,jj,jk) ) ) & + & ) + END IF + + z_rho_j(ji,jj,jk) = zcoef0 * ( rhd (ji,jj+1,jk) + rhd (ji,jj,jk) ) & + & * ( gde3w(ji,jj+1,jk) - gde3w(ji,jj,jk) ) + IF ( vmask(ji, jj-1, jk) > 0.5 .OR. vmask(ji, jj+1, jk) > 0.5 ) THEN + z_rho_j(ji,jj,jk) = z_rho_j(ji,jj,jk) - z_grav_10 * ( & + & ( zdrho_j (ji,jj+1,jk) - zdrho_j (ji,jj,jk) ) & + & * ( - gde3w(ji,jj+1,jk) + gde3w(ji,jj,jk) - z1_12 * ( zdz_j (ji,jj+1,jk) + zdz_j (ji,jj,jk) ) ) & + & - ( zdz_j (ji,jj+1,jk) - zdz_j (ji,jj,jk) ) & + & * ( rhd (ji,jj+1,jk) - rhd (ji,jj,jk) - z1_12 * ( zdrho_j(ji,jj+1,jk) + zdrho_j(ji,jj,jk) ) ) & + & ) + END IF + END_3D + + !-------------------------------------------------------------- + ! 8. Integrate in the vertical + !------------------------------------------------------------- + ! + ! --------------- + ! Surface value + ! --------------- + DO_2D( 0, 0, 0, 0 ) + zhpi(ji,jj,1) = ( z_rho_k(ji,jj,1) - z_rho_k(ji+1,jj ,1) - z_rho_i(ji,jj,1) ) * r1_e1u(ji,jj) + zhpj(ji,jj,1) = ( z_rho_k(ji,jj,1) - z_rho_k(ji ,jj+1,1) - z_rho_j(ji,jj,1) ) * r1_e2v(ji,jj) + ! add to the general momentum trend + puu(ji,jj,1,Krhs) = puu(ji,jj,1,Krhs) + zhpi(ji,jj,1) + pvv(ji,jj,1,Krhs) = pvv(ji,jj,1,Krhs) + zhpj(ji,jj,1) + END_2D + + ! ---------------- + ! interior value (2== -zdept(ji,jj,jk) ) THEN + jis = ji + 1; jid = ji + ELSE + jis = ji; jid = ji +1 + ENDIF + + ! integrate the pressure on the shallow side + jk1 = jk + DO WHILE ( -zdept(jis,jj,jk1) > zuijk ) + IF( jk1 == mbku(ji,jj) ) THEN + zuijk = -zdept(jis,jj,jk1) + EXIT + ENDIF + zdeps = MIN(zdept(jis,jj,jk1+1), -zuijk) + zpwes = zpwes + & + integ_spline(zdept(jis,jj,jk1), zdeps, & + asp(jis,jj,jk1), bsp(jis,jj,jk1), & + csp(jis,jj,jk1), dsp(jis,jj,jk1)) + jk1 = jk1 + 1 + END DO + + ! integrate the pressure on the deep side + jk1 = jk + DO WHILE ( -zdept(jid,jj,jk1) < zuijk ) + IF( jk1 == 1 ) THEN + zdeps = zdept(jid,jj,1) + MIN(zuijk, ssh(jid,jj,Kmm)*znad) + zrhdt1 = zrhh(jid,jj,1) - interp3(zdept(jid,jj,1), asp(jid,jj,1), & + bsp(jid,jj,1) , csp(jid,jj,1), & + dsp(jid,jj,1)) * zdeps + zpwed = zpwed + 0.5_wp * (zrhh(jid,jj,1) + zrhdt1) * zdeps + EXIT + ENDIF + zdeps = MAX(zdept(jid,jj,jk1-1), -zuijk) + zpwed = zpwed + & + integ_spline(zdeps, zdept(jid,jj,jk1), & + asp(jid,jj,jk1-1), bsp(jid,jj,jk1-1), & + csp(jid,jj,jk1-1), dsp(jid,jj,jk1-1) ) + jk1 = jk1 - 1 + END DO + + ! update the momentum trends in u direction + zdpdx1 = zcoef0 * r1_e1u(ji,jj) * ( zhpi(ji+1,jj,jk) - zhpi(ji,jj,jk) ) + IF( .NOT.ln_linssh ) THEN + zdpdx2 = zcoef0 * r1_e1u(ji,jj) * & + & ( REAL(jis-jid, wp) * (zpwes + zpwed) + (ssh(ji+1,jj,Kmm)-ssh(ji,jj,Kmm)) ) + ELSE + zdpdx2 = zcoef0 * r1_e1u(ji,jj) * REAL(jis-jid, wp) * (zpwes + zpwed) + ENDIF + puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) + (zdpdx1 + zdpdx2 - zpgu(ji,jj)) * umask(ji,jj,jk) + ENDIF + + !!!!! for v equation + IF( jk <= mbkv(ji,jj) ) THEN + IF( -zdept(ji,jj+1,jk) >= -zdept(ji,jj,jk) ) THEN + jjs = jj + 1; jjd = jj + ELSE + jjs = jj ; jjd = jj + 1 + ENDIF + + ! integrate the pressure on the shallow side + jk1 = jk + DO WHILE ( -zdept(ji,jjs,jk1) > zvijk ) + IF( jk1 == mbkv(ji,jj) ) THEN + zvijk = -zdept(ji,jjs,jk1) + EXIT + ENDIF + zdeps = MIN(zdept(ji,jjs,jk1+1), -zvijk) + zpnss = zpnss + & + integ_spline(zdept(ji,jjs,jk1), zdeps, & + asp(ji,jjs,jk1), bsp(ji,jjs,jk1), & + csp(ji,jjs,jk1), dsp(ji,jjs,jk1) ) + jk1 = jk1 + 1 + END DO + + ! integrate the pressure on the deep side + jk1 = jk + DO WHILE ( -zdept(ji,jjd,jk1) < zvijk ) + IF( jk1 == 1 ) THEN + zdeps = zdept(ji,jjd,1) + MIN(zvijk, ssh(ji,jjd,Kmm)*znad) + zrhdt1 = zrhh(ji,jjd,1) - interp3(zdept(ji,jjd,1), asp(ji,jjd,1), & + bsp(ji,jjd,1) , csp(ji,jjd,1), & + dsp(ji,jjd,1) ) * zdeps + zpnsd = zpnsd + 0.5_wp * (zrhh(ji,jjd,1) + zrhdt1) * zdeps + EXIT + ENDIF + zdeps = MAX(zdept(ji,jjd,jk1-1), -zvijk) + zpnsd = zpnsd + & + integ_spline(zdeps, zdept(ji,jjd,jk1), & + asp(ji,jjd,jk1-1), bsp(ji,jjd,jk1-1), & + csp(ji,jjd,jk1-1), dsp(ji,jjd,jk1-1) ) + jk1 = jk1 - 1 + END DO + + ! update the momentum trends in v direction + zdpdy1 = zcoef0 * r1_e2v(ji,jj) * ( zhpi(ji,jj+1,jk) - zhpi(ji,jj,jk) ) + IF( .NOT.ln_linssh ) THEN + zdpdy2 = zcoef0 * r1_e2v(ji,jj) * & + ( REAL(jjs-jjd, wp) * (zpnss + zpnsd) + (ssh(ji,jj+1,Kmm)-ssh(ji,jj,Kmm)) ) + ELSE + zdpdy2 = zcoef0 * r1_e2v(ji,jj) * REAL(jjs-jjd, wp) * (zpnss + zpnsd ) + ENDIF + + pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) + (zdpdy1 + zdpdy2 - zpgv(ji,jj)) * vmask(ji,jj,jk) + ENDIF + ! + END_3D + ! + ! + END SUBROUTINE hpg_prj + + + SUBROUTINE cspline( fsp, xsp, asp, bsp, csp, dsp, polynomial_type ) + !!---------------------------------------------------------------------- + !! *** ROUTINE cspline *** + !! + !! ** Purpose : constrained cubic spline interpolation + !! + !! ** Method : f(x) = asp + bsp*x + csp*x^2 + dsp*x^3 + !! + !! Reference: CJC Kruger, Constrained Cubic Spline Interpoltation + !!---------------------------------------------------------------------- + REAL(wp), DIMENSION(A2D(nn_hls),jpk), INTENT(in ) :: fsp, xsp ! value and coordinate + REAL(wp), DIMENSION(A2D(nn_hls),jpk), INTENT( out) :: asp, bsp, csp, dsp ! coefficients of the interpoated function + INTEGER , INTENT(in ) :: polynomial_type ! 1: cubic spline ; 2: Linear + ! + INTEGER :: ji, jj, jk ! dummy loop indices + REAL(wp) :: zdf1, zdf2, zddf1, zddf2, ztmp1, ztmp2, zdxtmp + REAL(wp) :: zdxtmp1, zdxtmp2, zalpha + REAL(wp), DIMENSION(jpk) :: zdf + !!---------------------------------------------------------------------- + ! + IF (polynomial_type == 1) THEN ! Constrained Cubic Spline + DO_2D( 1, 1, 1, 1 ) + !!Fritsch&Butland's method, 1984 (preferred, but more computation) + ! DO jk = 2, jpkm1-1 + ! zdxtmp1 = xsp(ji,jj,jk) - xsp(ji,jj,jk-1) + ! zdxtmp2 = xsp(ji,jj,jk+1) - xsp(ji,jj,jk) + ! zdf1 = ( fsp(ji,jj,jk) - fsp(ji,jj,jk-1) ) / zdxtmp1 + ! zdf2 = ( fsp(ji,jj,jk+1) - fsp(ji,jj,jk) ) / zdxtmp2 + ! + ! zalpha = ( zdxtmp1 + 2._wp * zdxtmp2 ) / ( zdxtmp1 + zdxtmp2 ) / 3._wp + ! + ! IF(zdf1 * zdf2 <= 0._wp) THEN + ! zdf(jk) = 0._wp + ! ELSE + ! zdf(jk) = zdf1 * zdf2 / ( ( 1._wp - zalpha ) * zdf1 + zalpha * zdf2 ) + ! ENDIF + ! END DO + + !!Simply geometric average + DO jk = 2, jpk-2 + zdf1 = (fsp(ji,jj,jk ) - fsp(ji,jj,jk-1)) / (xsp(ji,jj,jk ) - xsp(ji,jj,jk-1)) + zdf2 = (fsp(ji,jj,jk+1) - fsp(ji,jj,jk )) / (xsp(ji,jj,jk+1) - xsp(ji,jj,jk )) + + IF(zdf1 * zdf2 <= 0._wp) THEN + zdf(jk) = 0._wp + ELSE + zdf(jk) = 2._wp * zdf1 * zdf2 / (zdf1 + zdf2) + ENDIF + END DO + + zdf(1) = 1.5_wp * ( fsp(ji,jj,2) - fsp(ji,jj,1) ) / & + & ( xsp(ji,jj,2) - xsp(ji,jj,1) ) - 0.5_wp * zdf(2) + zdf(jpkm1) = 1.5_wp * ( fsp(ji,jj,jpkm1) - fsp(ji,jj,jpkm1-1) ) / & + & ( xsp(ji,jj,jpkm1) - xsp(ji,jj,jpkm1-1) ) - 0.5_wp * zdf(jpk - 2) + + DO jk = 1, jpk-2 + zdxtmp = xsp(ji,jj,jk+1) - xsp(ji,jj,jk) + ztmp1 = (zdf(jk+1) + 2._wp * zdf(jk)) / zdxtmp + ztmp2 = 6._wp * (fsp(ji,jj,jk+1) - fsp(ji,jj,jk)) / zdxtmp / zdxtmp + zddf1 = -2._wp * ztmp1 + ztmp2 + ztmp1 = (2._wp * zdf(jk+1) + zdf(jk)) / zdxtmp + zddf2 = 2._wp * ztmp1 - ztmp2 + + dsp(ji,jj,jk) = (zddf2 - zddf1) / 6._wp / zdxtmp + csp(ji,jj,jk) = ( xsp(ji,jj,jk+1) * zddf1 - xsp(ji,jj,jk)*zddf2 ) / 2._wp / zdxtmp + bsp(ji,jj,jk) = ( fsp(ji,jj,jk+1) - fsp(ji,jj,jk) ) / zdxtmp - & + & csp(ji,jj,jk) * ( xsp(ji,jj,jk+1) + xsp(ji,jj,jk) ) - & + & dsp(ji,jj,jk) * ((xsp(ji,jj,jk+1) + xsp(ji,jj,jk))**2 - & + & xsp(ji,jj,jk+1) * xsp(ji,jj,jk)) + asp(ji,jj,jk) = fsp(ji,jj,jk) - xsp(ji,jj,jk) * (bsp(ji,jj,jk) + & + & (xsp(ji,jj,jk) * (csp(ji,jj,jk) + & + & dsp(ji,jj,jk) * xsp(ji,jj,jk)))) + END DO + END_2D + + ELSEIF ( polynomial_type == 2 ) THEN ! Linear + DO_3D( 1, 1, 1, 1, 1, jpk-2 ) + zdxtmp =xsp(ji,jj,jk+1) - xsp(ji,jj,jk) + ztmp1 = fsp(ji,jj,jk+1) - fsp(ji,jj,jk) + + dsp(ji,jj,jk) = 0._wp + csp(ji,jj,jk) = 0._wp + bsp(ji,jj,jk) = ztmp1 / zdxtmp + asp(ji,jj,jk) = fsp(ji,jj,jk) - bsp(ji,jj,jk) * xsp(ji,jj,jk) + END_3D + ! + ELSE + CALL ctl_stop( 'invalid polynomial type in cspline' ) + ENDIF + ! + END SUBROUTINE cspline + + + FUNCTION interp1(x, xl, xr, fl, fr) RESULT(f) + !!---------------------------------------------------------------------- + !! *** ROUTINE interp1 *** + !! + !! ** Purpose : 1-d linear interpolation + !! + !! ** Method : interpolation is straight forward + !! extrapolation is also permitted (no value limit) + !!---------------------------------------------------------------------- + REAL(wp), INTENT(in) :: fl, fr + REAL(dp), INTENT(in) :: x, xl, xr + REAL(wp) :: f ! result of the interpolation (extrapolation) + REAL(wp) :: zdeltx + !!---------------------------------------------------------------------- + ! + zdeltx = xr - xl + IF( abs(zdeltx) <= 10._wp * EPSILON(x) ) THEN + f = 0.5_wp * (fl + fr) + ELSE + f = ( (x - xl ) * fr - ( x - xr ) * fl ) / zdeltx + ENDIF + ! + END FUNCTION interp1 + + + FUNCTION interp2( x, a, b, c, d ) RESULT(f) + !!---------------------------------------------------------------------- + !! *** ROUTINE interp1 *** + !! + !! ** Purpose : 1-d constrained cubic spline interpolation + !! + !! ** Method : cubic spline interpolation + !! + !!---------------------------------------------------------------------- + REAL(wp), INTENT(in) :: x, a, b, c, d + REAL(wp) :: f ! value from the interpolation + !!---------------------------------------------------------------------- + ! + f = a + x* ( b + x * ( c + d * x ) ) + ! + END FUNCTION interp2 + + + FUNCTION interp3( x, a, b, c, d ) RESULT(f) + !!---------------------------------------------------------------------- + !! *** ROUTINE interp1 *** + !! + !! ** Purpose : Calculate the first order of derivative of + !! a cubic spline function y=a+b*x+c*x^2+d*x^3 + !! + !! ** Method : f=dy/dx=b+2*c*x+3*d*x^2 + !! + !!---------------------------------------------------------------------- + REAL(wp), INTENT(in) :: x, a, b, c, d + REAL(wp) :: f ! value from the interpolation + !!---------------------------------------------------------------------- + ! + f = b + x * ( 2._wp * c + 3._wp * d * x) + ! + END FUNCTION interp3 + + + FUNCTION integ_spline( xl, xr, a, b, c, d ) RESULT(f) + !!---------------------------------------------------------------------- + !! *** ROUTINE interp1 *** + !! + !! ** Purpose : 1-d constrained cubic spline integration + !! + !! ** Method : integrate polynomial a+bx+cx^2+dx^3 from xl to xr + !! + !!---------------------------------------------------------------------- + REAL(wp), INTENT(in) :: xl, xr, a, b, c, d + REAL(wp) :: za1, za2, za3 + REAL(wp) :: f ! integration result + !!---------------------------------------------------------------------- + ! + za1 = 0.5_wp * b + za2 = c / 3.0_wp + za3 = 0.25_wp * d + ! + f = xr * ( a + xr * ( za1 + xr * ( za2 + za3 * xr ) ) ) - & + & xl * ( a + xl * ( za1 + xl * ( za2 + za3 * xl ) ) ) + ! + END FUNCTION integ_spline + + !!====================================================================== +END MODULE dynhpg diff --git a/src/OCE/TRA/trabbl.F90_dev b/src/OCE/TRA/trabbl.F90_dev new file mode 100755 index 0000000..d58b1d1 --- /dev/null +++ b/src/OCE/TRA/trabbl.F90_dev @@ -0,0 +1,654 @@ +MODULE trabbl + !!============================================================================== + !! *** MODULE trabbl *** + !! Ocean physics : advective and/or diffusive bottom boundary layer scheme + !!============================================================================== + !! History : OPA ! 1996-06 (L. Mortier) Original code + !! 8.0 ! 1997-11 (G. Madec) Optimization + !! NEMO 1.0 ! 2002-08 (G. Madec) free form + modules + !! - ! 2004-01 (A. de Miranda, G. Madec, J.M. Molines ) add advective bbl + !! 3.3 ! 2009-11 (G. Madec) merge trabbl and trabbl_adv + style + optimization + !! - ! 2010-04 (G. Madec) Campin & Goosse advective bbl + !! - ! 2010-06 (C. Ethe, G. Madec) merge TRA-TRC + !! - ! 2010-11 (G. Madec) add mbk. arrays associated to the deepest ocean level + !! - ! 2013-04 (F. Roquet, G. Madec) use of eosbn2 instead of local hard coded alpha and beta + !! 4.0 ! 2017-04 (G. Madec) ln_trabbl namelist variable instead of a CPP key + !!---------------------------------------------------------------------- + + !!---------------------------------------------------------------------- + !! tra_bbl_alloc : allocate trabbl arrays + !! tra_bbl : update the tracer trends due to the bottom boundary layer (advective and/or diffusive) + !! tra_bbl_dif : generic routine to compute bbl diffusive trend + !! tra_bbl_adv : generic routine to compute bbl advective trend + !! bbl : computation of bbl diffu. flux coef. & transport in bottom boundary layer + !! tra_bbl_init : initialization, namelist read, parameters control + !!---------------------------------------------------------------------- + USE oce ! ocean dynamics and active tracers + USE dom_oce ! ocean space and time domain + USE phycst ! physical constant + USE eosbn2 ! equation of state + USE trd_oce ! trends: ocean variables + USE trdtra ! trends: active tracers + ! + USE iom ! IOM library + USE in_out_manager ! I/O manager + USE lbclnk ! ocean lateral boundary conditions + USE prtctl ! Print control + USE timing ! Timing + USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) + !!GMDDEV + USE zdfdrg + !!GMDDEV + IMPLICIT NONE + PRIVATE + + PUBLIC tra_bbl ! routine called by step.F90 + PUBLIC tra_bbl_init ! routine called by nemogcm.F90 + PUBLIC tra_bbl_dif ! routine called by trcbbl.F90 + PUBLIC tra_bbl_adv ! - - - + PUBLIC bbl ! routine called by trcbbl.F90 and dtadyn.F90 + + ! !!* Namelist nambbl * + LOGICAL , PUBLIC :: ln_trabbl !: bottom boundary layer flag + INTEGER , PUBLIC :: nn_bbl_ldf !: =1 : diffusive bbl or not (=0) + INTEGER , PUBLIC :: nn_bbl_adv !: =1/2 : advective bbl or not (=0) + ! ! =1 : advective bbl using the bottom ocean velocity + ! ! =2 : - - using utr_bbl proportional to grad(rho) + REAL(wp), PUBLIC :: rn_ahtbbl !: along slope bbl diffusive coefficient [m2/s] + REAL(wp), PUBLIC :: rn_gambbl !: lateral coeff. for bottom boundary layer scheme [s] + + REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:), PUBLIC :: utr_bbl , vtr_bbl ! u- (v-) transport in the bottom boundary layer + REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:), PUBLIC :: ahu_bbl , ahv_bbl ! masked diffusive bbl coeff. at u & v-pts + + INTEGER , ALLOCATABLE, SAVE, DIMENSION(:,:), PUBLIC :: mbku_d , mbkv_d ! vertical index of the "lower" bottom ocean U/V-level (PUBLIC for TAM) + INTEGER , ALLOCATABLE, SAVE, DIMENSION(:,:), PUBLIC :: mgrhu , mgrhv ! = +/-1, sign of grad(H) in u-(v-)direction (PUBLIC for TAM) + !!GMDDEV + REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:), PUBLIC :: grhu,grhv,dprgu,dprgv,deltax,deltay + !!GMDDEV + REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: ahu_bbl_0, ahv_bbl_0 ! diffusive bbl flux coefficients at u and v-points + REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:), PUBLIC :: e3u_bbl_0, e3v_bbl_0 ! thichness of the bbl (e3) at u and v-points (PUBLIC for TAM) + + !! * Substitutions +# include "do_loop_substitute.h90" +# include "domzgr_substitute.h90" + !!---------------------------------------------------------------------- + !! NEMO/OCE 4.0 , NEMO Consortium (2018) + !! $Id: trabbl.F90 15053 2021-06-24 15:39:38Z clem $ + !! Software governed by the CeCILL license (see ./LICENSE) + !!---------------------------------------------------------------------- +CONTAINS + + INTEGER FUNCTION tra_bbl_alloc() + !!---------------------------------------------------------------------- + !! *** FUNCTION tra_bbl_alloc *** + !!---------------------------------------------------------------------- + ALLOCATE( utr_bbl (jpi,jpj) , ahu_bbl (jpi,jpj) , mbku_d(jpi,jpj) , mgrhu(jpi,jpj) , & + & vtr_bbl (jpi,jpj) , ahv_bbl (jpi,jpj) , mbkv_d(jpi,jpj) , mgrhv(jpi,jpj) , & + & ahu_bbl_0(jpi,jpj) , ahv_bbl_0(jpi,jpj) , & + & e3u_bbl_0(jpi,jpj) , e3v_bbl_0(jpi,jpj) , & +!!GMDDEV + & grhu(jpi,jpj), grhv(jpi,jpj), dprgu(jpi,jpj), dprgv(jpi,jpj),deltax(jpi,jpj),deltay(jpi,jpj),STAT=tra_bbl_alloc ) +!!GMDDEV + ! + CALL mpp_sum ( 'trabbl', tra_bbl_alloc ) + IF( tra_bbl_alloc > 0 ) CALL ctl_warn('tra_bbl_alloc: allocation of arrays failed.') + END FUNCTION tra_bbl_alloc + + + SUBROUTINE tra_bbl( kt, Kbb, Kmm, pts, Krhs ) + !!---------------------------------------------------------------------- + !! *** ROUTINE bbl *** + !! + !! ** Purpose : Compute the before tracer (t & s) trend associated + !! with the bottom boundary layer and add it to the general + !! trend of tracer equations. + !! + !! ** Method : Depending on namtra_bbl namelist parameters the bbl + !! diffusive and/or advective contribution to the tracer trend + !! is added to the general tracer trend + !!---------------------------------------------------------------------- + INTEGER, INTENT(in ) :: kt ! ocean time-step + INTEGER, INTENT(in ) :: Kbb, Kmm, Krhs ! time level indices + REAL(dp), DIMENSION(jpi,jpj,jpk,jpts,jpt), INTENT(inout) :: pts ! active tracers and RHS of tracer equation + ! + INTEGER :: ji, jj, jk ! Dummy loop indices + REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: ztrdt, ztrds + !!---------------------------------------------------------------------- + ! + IF( ln_timing ) CALL timing_start( 'tra_bbl') + ! + IF( l_trdtra ) THEN !* Save the T-S input trends + ALLOCATE( ztrdt(jpi,jpj,jpk), ztrds(jpi,jpj,jpk) ) + ztrdt(:,:,:) = pts(:,:,:,jp_tem,Krhs) + ztrds(:,:,:) = pts(:,:,:,jp_sal,Krhs) + ENDIF + + IF( nn_bbl_ldf == 1 ) THEN !* Diffusive bbl + ! + CALL tra_bbl_dif( pts(:,:,:,:,Kbb), pts(:,:,:,:,Krhs), jpts, Kmm ) + IF( sn_cfctl%l_prtctl ) & + CALL prt_ctl( tab3d_1=pts(:,:,:,jp_tem,Krhs), clinfo1=' bbl_ldf - Ta: ', mask1=tmask, & + & tab3d_2=pts(:,:,:,jp_sal,Krhs), clinfo2= ' Sa: ', mask2=tmask, clinfo3='tra' ) + CALL iom_put( "ahu_bbl", ahu_bbl ) ! bbl diffusive flux i-coef + CALL iom_put( "ahv_bbl", ahv_bbl ) ! bbl diffusive flux j-coef + ! + ENDIF + ! + IF( nn_bbl_adv /= 0 ) THEN !* Advective bbl + ! + CALL tra_bbl_adv( pts(:,:,:,:,Kbb), pts(:,:,:,:,Krhs), jpts, Kmm ) + IF(sn_cfctl%l_prtctl) & + CALL prt_ctl( tab3d_1=pts(:,:,:,jp_tem,Krhs), clinfo1=' bbl_adv - Ta: ', mask1=tmask, & + & tab3d_2=pts(:,:,:,jp_sal,Krhs), clinfo2= ' Sa: ', mask2=tmask, clinfo3='tra' ) + CALL iom_put( "uoce_bbl", utr_bbl ) ! bbl i-transport + CALL iom_put( "voce_bbl", vtr_bbl ) ! bbl j-transport + ! + ENDIF + + IF( l_trdtra ) THEN ! send the trends for further diagnostics + ztrdt(:,:,:) = pts(:,:,:,jp_tem,Krhs) - ztrdt(:,:,:) + ztrds(:,:,:) = pts(:,:,:,jp_sal,Krhs) - ztrds(:,:,:) + CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_tem, jptra_bbl, ztrdt ) + CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_sal, jptra_bbl, ztrds ) + DEALLOCATE( ztrdt, ztrds ) + ENDIF + ! + IF( ln_timing ) CALL timing_stop( 'tra_bbl') + ! + END SUBROUTINE tra_bbl + + + SUBROUTINE tra_bbl_dif( pt, pt_rhs, kjpt, Kmm ) + !!---------------------------------------------------------------------- + !! *** ROUTINE tra_bbl_dif *** + !! + !! ** Purpose : Computes the bottom boundary horizontal and vertical + !! advection terms. + !! + !! ** Method : * diffusive bbl only (nn_bbl_ldf=1) : + !! When the product grad( rho) * grad(h) < 0 (where grad is an + !! along bottom slope gradient) an additional lateral 2nd order + !! diffusion along the bottom slope is added to the general + !! tracer trend, otherwise the additional trend is set to 0. + !! A typical value of ahbt is 2000 m2/s (equivalent to + !! a downslope velocity of 20 cm/s if the condition for slope + !! convection is satified) + !! + !! ** Action : pt_rhs increased by the bbl diffusive trend + !! + !! References : Beckmann, A., and R. Doscher, 1997, J. Phys.Oceanogr., 581-591. + !! Campin, J.-M., and H. Goosse, 1999, Tellus, 412-430. + !!---------------------------------------------------------------------- + INTEGER , INTENT(in ) :: kjpt ! number of tracers + REAL(dp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in ) :: pt ! before and now tracer fields + REAL(dp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) :: pt_rhs ! tracer trend + INTEGER , INTENT(in ) :: Kmm ! time level indices + ! + INTEGER :: ji, jj, jn ! dummy loop indices + INTEGER :: ik ! local integers + REAL(wp) :: zbtr ! local scalars + REAL(wp), DIMENSION(A2D(nn_hls)) :: zptb ! workspace + !!---------------------------------------------------------------------- + ! + DO jn = 1, kjpt ! tracer loop + ! ! =========== + DO_2D( 1, 1, 1, 1 ) + ik = mbkt(ji,jj) ! bottom T-level index + zptb(ji,jj) = pt(ji,jj,ik,jn) ! bottom before T and S + END_2D + ! + DO_2D( 0, 0, 0, 0 ) ! Compute the trend + ik = mbkt(ji,jj) ! bottom T-level index + pt_rhs(ji,jj,ik,jn) = pt_rhs(ji,jj,ik,jn) & + & + ( ahu_bbl(ji ,jj ) * ( zptb(ji+1,jj ) - zptb(ji ,jj ) ) & + & - ahu_bbl(ji-1,jj ) * ( zptb(ji ,jj ) - zptb(ji-1,jj ) ) & + & + ahv_bbl(ji ,jj ) * ( zptb(ji ,jj+1) - zptb(ji ,jj ) ) & + & - ahv_bbl(ji ,jj-1) * ( zptb(ji ,jj ) - zptb(ji ,jj-1) ) ) & + & * r1_e1e2t(ji,jj) / e3t(ji,jj,ik,Kmm) + END_2D + ! ! =========== + END DO ! end tracer + ! ! =========== + END SUBROUTINE tra_bbl_dif + + + ! NOTE: [tiling] tiling changes the results, but only the order of floating point operations is different + SUBROUTINE tra_bbl_adv( pt, pt_rhs, kjpt, Kmm ) + !!---------------------------------------------------------------------- + !! *** ROUTINE trc_bbl *** + !! + !! ** Purpose : Compute the before passive tracer trend associated + !! with the bottom boundary layer and add it to the general trend + !! of tracer equations. + !! ** Method : advective bbl (nn_bbl_adv = 1 or 2) : + !! nn_bbl_adv = 1 use of the ocean near bottom velocity as bbl velocity + !! nn_bbl_adv = 2 follow Campin and Goosse (1999) implentation i.e. + !! transport proportional to the along-slope density gradient + !! + !! References : Beckmann, A., and R. Doscher, 1997, J. Phys.Oceanogr., 581-591. + !! Campin, J.-M., and H. Goosse, 1999, Tellus, 412-430. + !!---------------------------------------------------------------------- + INTEGER , INTENT(in ) :: kjpt ! number of tracers + REAL(dp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in ) :: pt ! before and now tracer fields + REAL(dp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) :: pt_rhs ! tracer trend + INTEGER , INTENT(in ) :: Kmm ! time level indices + ! + INTEGER :: ji, jj, jk, jn ! dummy loop indices + INTEGER :: iis , iid , ijs , ijd ! local integers + INTEGER :: ikus, ikud, ikvs, ikvd ! - - + REAL(wp) :: zbtr, ztra ! local scalars + REAL(wp) :: zu_bbl, zv_bbl ! - - + !!---------------------------------------------------------------------- + ! ! =========== + DO jn = 1, kjpt ! tracer loop + ! ! =========== + DO_2D_OVR( 1, 0, 1, 0 ) ! CAUTION start from i=1 to update i=2 when cyclic east-west + IF( utr_bbl(ji,jj) /= 0.e0 ) THEN ! non-zero i-direction bbl advection + ! down-slope i/k-indices (deep) & up-slope i/k indices (shelf) + iid = ji + MAX( 0, mgrhu(ji,jj) ) ; iis = ji + 1 - MAX( 0, mgrhu(ji,jj) ) + ikud = mbku_d(ji,jj) ; ikus = mbku(ji,jj) + zu_bbl = ABS( utr_bbl(ji,jj) ) + ! + ! ! up -slope T-point (shelf bottom point) + zbtr = r1_e1e2t(iis,jj) / e3t(iis,jj,ikus,Kmm) + ztra = zu_bbl * ( pt(iid,jj,ikus,jn) - pt(iis,jj,ikus,jn) ) * zbtr + pt_rhs(iis,jj,ikus,jn) = pt_rhs(iis,jj,ikus,jn) + ztra + ! + DO jk = ikus, ikud-1 ! down-slope upper to down T-point (deep column) + zbtr = r1_e1e2t(iid,jj) / e3t(iid,jj,jk,Kmm) + ztra = zu_bbl * ( pt(iid,jj,jk+1,jn) - pt(iid,jj,jk,jn) ) * zbtr + pt_rhs(iid,jj,jk,jn) = pt_rhs(iid,jj,jk,jn) + ztra + END DO + ! + zbtr = r1_e1e2t(iid,jj) / e3t(iid,jj,ikud,Kmm) + ztra = zu_bbl * ( pt(iis,jj,ikus,jn) - pt(iid,jj,ikud,jn) ) * zbtr + pt_rhs(iid,jj,ikud,jn) = pt_rhs(iid,jj,ikud,jn) + ztra + ENDIF + ! + IF( vtr_bbl(ji,jj) /= 0.e0 ) THEN ! non-zero j-direction bbl advection + ! down-slope j/k-indices (deep) & up-slope j/k indices (shelf) + ijd = jj + MAX( 0, mgrhv(ji,jj) ) ; ijs = jj + 1 - MAX( 0, mgrhv(ji,jj) ) + ikvd = mbkv_d(ji,jj) ; ikvs = mbkv(ji,jj) + zv_bbl = ABS( vtr_bbl(ji,jj) ) + ! + ! up -slope T-point (shelf bottom point) + zbtr = r1_e1e2t(ji,ijs) / e3t(ji,ijs,ikvs,Kmm) + ztra = zv_bbl * ( pt(ji,ijd,ikvs,jn) - pt(ji,ijs,ikvs,jn) ) * zbtr + pt_rhs(ji,ijs,ikvs,jn) = pt_rhs(ji,ijs,ikvs,jn) + ztra + ! + DO jk = ikvs, ikvd-1 ! down-slope upper to down T-point (deep column) + zbtr = r1_e1e2t(ji,ijd) / e3t(ji,ijd,jk,Kmm) + ztra = zv_bbl * ( pt(ji,ijd,jk+1,jn) - pt(ji,ijd,jk,jn) ) * zbtr + pt_rhs(ji,ijd,jk,jn) = pt_rhs(ji,ijd,jk,jn) + ztra + END DO + ! ! down-slope T-point (deep bottom point) + zbtr = r1_e1e2t(ji,ijd) / e3t(ji,ijd,ikvd,Kmm) + ztra = zv_bbl * ( pt(ji,ijs,ikvs,jn) - pt(ji,ijd,ikvd,jn) ) * zbtr + pt_rhs(ji,ijd,ikvd,jn) = pt_rhs(ji,ijd,ikvd,jn) + ztra + ENDIF + END_2D + ! ! =========== + END DO ! end tracer + ! ! =========== + END SUBROUTINE tra_bbl_adv + + + SUBROUTINE bbl( kt, kit000, Kbb, Kmm ) + !!---------------------------------------------------------------------- + !! *** ROUTINE bbl *** + !! + !! ** Purpose : Computes the bottom boundary horizontal and vertical + !! advection terms. + !! + !! ** Method : * diffusive bbl (nn_bbl_ldf=1) : + !! When the product grad( rho) * grad(h) < 0 (where grad is an + !! along bottom slope gradient) an additional lateral 2nd order + !! diffusion along the bottom slope is added to the general + !! tracer trend, otherwise the additional trend is set to 0. + !! A typical value of ahbt is 2000 m2/s (equivalent to + !! a downslope velocity of 20 cm/s if the condition for slope + !! convection is satified) + !! * advective bbl (nn_bbl_adv=1 or 2) : + !! nn_bbl_adv = 1 use of the ocean velocity as bbl velocity + !! nn_bbl_adv = 2 follow Campin and Goosse (1999) implentation + !! i.e. transport proportional to the along-slope density gradient + !! + !! NB: the along slope density gradient is evaluated using the + !! local density (i.e. referenced at a common local depth). + !! + !! References : Beckmann, A., and R. Doscher, 1997, J. Phys.Oceanogr., 581-591. + !! Campin, J.-M., and H. Goosse, 1999, Tellus, 412-430. + !!---------------------------------------------------------------------- + INTEGER , INTENT(in ) :: kt ! ocean time-step index + INTEGER , INTENT(in ) :: kit000 ! first time step index + INTEGER , INTENT(in ) :: Kbb, Kmm ! ocean time level index + ! + INTEGER :: ji, jj ! dummy loop indices + INTEGER :: ik ! local integers + INTEGER :: iis, iid ! - - + INTEGER :: ijs, ijd ! - - + REAL(wp) :: za, zb, zgdrho ! local scalars + REAL(wp) :: zsign, zsigna, zgbbl ! - - + !!GMDDEV + REAL(wp) :: zbblvel, zhx, zhy, zx_flag, zy_flag + !!GMDDEV + REAL(wp), DIMENSION(A2D(nn_hls),jpts) :: zts, zab ! 3D workspace + REAL(wp), DIMENSION(A2D(nn_hls)) :: zub, zvb, zdep ! 2D workspace + !!GMDDEV + REAL(wp), DIMENSION(A2D(nn_hls)) :: zpvelu, zpvelv, zvgx, zvgy,zcdragx,zcdragy,zu_bbl,zv_bbl + !!GMDDEV + !!---------------------------------------------------------------------- + ! + IF( .NOT. l_istiled .OR. ntile == 1 ) THEN ! Do only on the first tile + IF( kt == kit000 ) THEN + IF(lwp) WRITE(numout,*) + IF(lwp) WRITE(numout,*) 'trabbl:bbl : Compute bbl velocities and diffusive coefficients' + IF(lwp) WRITE(numout,*) '~~~~~~~~~~' + ENDIF + ENDIF + ! !* bottom variables (T, S, alpha, beta, depth, velocity) + DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) + ik = mbkt(ji,jj) ! bottom T-level index + zts (ji,jj,jp_tem) = ts(ji,jj,ik,jp_tem,Kbb) ! bottom before T and S + zts (ji,jj,jp_sal) = ts(ji,jj,ik,jp_sal,Kbb) + ! + zdep(ji,jj) = gdept(ji,jj,ik,Kmm) ! bottom T-level reference depth + zub (ji,jj) = uu(ji,jj,mbku(ji,jj),Kmm) ! bottom velocity + zvb (ji,jj) = vv(ji,jj,mbkv(ji,jj),Kmm) + END_2D + ! + CALL eos_rab( zts, zdep, zab, Kmm ) + ! + ! !-------------------! + IF( nn_bbl_ldf == 1 ) THEN ! diffusive bbl ! + ! !-------------------! + DO_2D_OVR( 1, 0, 1, 0 ) ! (criteria for non zero flux: grad(rho).grad(h) < 0 ) + ! ! i-direction + za = zab(ji+1,jj,jp_tem) + zab(ji,jj,jp_tem) ! 2*(alpha,beta) at u-point + zb = zab(ji+1,jj,jp_sal) + zab(ji,jj,jp_sal) + ! ! 2*masked bottom density gradient + zgdrho = ( za * ( zts(ji+1,jj,jp_tem) - zts(ji,jj,jp_tem) ) & + & - zb * ( zts(ji+1,jj,jp_sal) - zts(ji,jj,jp_sal) ) ) * umask(ji,jj,1) + ! + zsign = SIGN( 0.5_wp, -zgdrho * REAL( mgrhu(ji,jj) ) ) ! sign of ( i-gradient * i-slope ) + ahu_bbl(ji,jj) = ( 0.5 - zsign ) * ahu_bbl_0(ji,jj) ! masked diffusive flux coeff. + ! + ! ! j-direction + za = zab(ji,jj+1,jp_tem) + zab(ji,jj,jp_tem) ! 2*(alpha,beta) at v-point + zb = zab(ji,jj+1,jp_sal) + zab(ji,jj,jp_sal) + ! ! 2*masked bottom density gradient + zgdrho = ( za * ( zts(ji,jj+1,jp_tem) - zts(ji,jj,jp_tem) ) & + & - zb * ( zts(ji,jj+1,jp_sal) - zts(ji,jj,jp_sal) ) ) * vmask(ji,jj,1) + ! + zsign = SIGN( 0.5_wp, -zgdrho * REAL( mgrhv(ji,jj) ) ) ! sign of ( j-gradient * j-slope ) + ahv_bbl(ji,jj) = ( 0.5 - zsign ) * ahv_bbl_0(ji,jj) + END_2D + ! + ENDIF + ! + ! !-------------------! + IF( nn_bbl_adv /= 0 ) THEN ! advective bbl ! + ! !-------------------! + SELECT CASE ( nn_bbl_adv ) !* bbl transport type + ! + CASE( 1 ) != use of upper velocity + DO_2D_OVR( 1, 0, 1, 0 ) ! criteria: grad(rho).grad(h)<0 and grad(rho).grad(h)<0 + ! ! i-direction + za = zab(ji+1,jj,jp_tem) + zab(ji,jj,jp_tem) ! 2*(alpha,beta) at u-point + zb = zab(ji+1,jj,jp_sal) + zab(ji,jj,jp_sal) + ! ! 2*masked bottom density gradient + zgdrho = ( za * ( zts(ji+1,jj,jp_tem) - zts(ji,jj,jp_tem) ) & + - zb * ( zts(ji+1,jj,jp_sal) - zts(ji,jj,jp_sal) ) ) * umask(ji,jj,1) + ! + zsign = SIGN( 0.5_wp, - zgdrho * REAL( mgrhu(ji,jj) ) ) ! sign of i-gradient * i-slope + zsigna= SIGN( 0.5_wp, zub(ji,jj) * REAL( mgrhu(ji,jj) ) ) ! sign of u * i-slope + ! + ! ! bbl velocity + utr_bbl(ji,jj) = ( 0.5 + zsigna ) * ( 0.5 - zsign ) * e2u(ji,jj) * e3u_bbl_0(ji,jj) * zub(ji,jj) + ! + ! ! j-direction + za = zab(ji,jj+1,jp_tem) + zab(ji,jj,jp_tem) ! 2*(alpha,beta) at v-point + zb = zab(ji,jj+1,jp_sal) + zab(ji,jj,jp_sal) + ! ! 2*masked bottom density gradient + zgdrho = ( za * ( zts(ji,jj+1,jp_tem) - zts(ji,jj,jp_tem) ) & + & - zb * ( zts(ji,jj+1,jp_sal) - zts(ji,jj,jp_sal) ) ) * vmask(ji,jj,1) + zsign = SIGN( 0.5_wp, - zgdrho * REAL( mgrhv(ji,jj) ) ) ! sign of j-gradient * j-slope + zsigna= SIGN( 0.5_wp, zvb(ji,jj) * REAL( mgrhv(ji,jj) ) ) ! sign of u * i-slope + ! + ! ! bbl transport + vtr_bbl(ji,jj) = ( 0.5 + zsigna ) * ( 0.5 - zsign ) * e1v(ji,jj) * e3v_bbl_0(ji,jj) * zvb(ji,jj) + END_2D + ! + CASE( 2 ) != bbl velocity = F( delta rho ) + !!GMDDEV + dprgu(:,:)=0.e0 + dprgv(:,:)=0.e0 + deltax(:,:)=0.e0 + deltay(:,:)=0.e0 + ahu_bbl_0(:,:)=0.e0 + ahv_bbl_0(:,:)=0.e0 + zpvelu(:,:)=0.e0 + zpvelv(:,:)=0.e0 + zvgx(:,:)=0.e0 + zvgy(:,:)=0.e0 + zcdragx(:,:)=0.e0 + zcdragy(:,:)=0.e0 + zu_bbl(:,:)=0.e0 + zv_bbl(:,:)=0.e0 + utr_bbl(:,:)=0.e0 + vtr_bbl(:,:)=0.e0 + !!GMDDEV + DO_2D_OVR( 1, 0, 1, 0 ) ! criteria: rho_up > rho_down + ! ! i-direction + ! down-slope T-point i/k-index (deep) & up-slope T-point i/k-index (shelf) + iid = ji + MAX( 0, mgrhu(ji,jj) ) + iis = ji + 1 - MAX( 0, mgrhu(ji,jj) ) + ! + ! + za = zab(ji+1,jj,jp_tem) + zab(ji,jj,jp_tem) ! 2*(alpha,beta) at u-point + zb = zab(ji+1,jj,jp_sal) + zab(ji,jj,jp_sal) + ! ! masked bottom density gradient + zgdrho = 0.5 * ( za * ( zts(iid,jj,jp_tem) - zts(iis,jj,jp_tem) ) & + & - zb * ( zts(iid,jj,jp_sal) - zts(iis,jj,jp_sal) ) ) * umask(ji,jj,1) + zgdrho = MAX( 0.e0, zgdrho ) ! only if shelf is denser than deep + !!GMDDEV + dprgu(ji,jj)=grav*zgdrho*grhu(ji,jj) + zvgy(ji,jj)=-dprgu(ji,jj)/ff_t(ji,jj) + zhx=MAX(e3u_bbl_0(ji,jj),1.e-1) + zcdragx(ji,jj)=MIN(0.1,( vkarmn / LOG( 0.5*zhx / 3.4e-2 ) )**2) + deltax(ji,jj)=ABS(zcdragx(ji,jj)*zvgy(ji,jj)/ff_t(ji,jj))/zhx + zbblvel=ABS(((deltax(ji,jj))/(1.+(deltax(ji,jj))**2.))*zvgy(ji,jj)) + zpvelv(ji,jj)=(1./(1.+(deltax(ji,jj))**2.))*zvgy(ji,jj) + ahu_bbl_0(ji,jj) = zbblvel * e1u(ji,jj) * e2_e1u(ji,jj) * e3u_bbl_0(ji,jj) * umask(ji,jj,1) + !!GMDDEV + ! + ! + ! ! j-direction + ! down-slope T-point j/k-index (deep) & of the up -slope T-point j/k-index (shelf) + ijd = jj + MAX( 0, mgrhv(ji,jj) ) + ijs = jj + 1 - MAX( 0, mgrhv(ji,jj) ) + ! + ! + za = zab(ji,jj+1,jp_tem) + zab(ji,jj,jp_tem) ! 2*(alpha,beta) at v-point + zb = zab(ji,jj+1,jp_sal) + zab(ji,jj,jp_sal) + ! ! masked bottom density gradient + zgdrho = 0.5 * ( za * ( zts(ji,ijd,jp_tem) - zts(ji,ijs,jp_tem) ) & + & - zb * ( zts(ji,ijd,jp_sal) - zts(ji,ijs,jp_sal) ) ) * vmask(ji,jj,1) + zgdrho = MAX( 0.e0, zgdrho ) ! only if shelf is denser than deep + !!GMDDEV + dprgv(ji,jj)=grav*zgdrho*grhv(ji,jj) + zvgx(ji,jj)=dprgv(ji,jj)/ff_t(ji,jj) + zhy=MAX(e3v_bbl_0(ji,jj),1.e-1) + zcdragy(ji,jj)=MIN(0.1,( vkarmn / LOG( 0.5*zhy / 3.4e-2 ) )**2) + deltay(ji,jj)=ABS(zcdragy(ji,jj)*zvgx(ji,jj)/ff_t(ji,jj))/zhy + zbblvel=ABS(((deltay(ji,jj))/(1.+(deltay(ji,jj))**2.))*zvgx(ji,jj)) + zpvelu(ji,jj)=(1./(1.+(deltay(ji,jj))**2.))*zvgx(ji,jj) + ahv_bbl_0(ji,jj) = zbblvel * e2v(ji,jj) * e1_e2v(ji,jj) * e3v_bbl_0(ji,jj) * vmask(ji,jj,1) + !!GMDDEV + ! + END_2D + !!GMDDEV + + CALL lbc_lnk( 'trabbl', dprgu,'U',1.0_wp) + CALL lbc_lnk( 'trabbl', dprgv,'V',1.0_wp) + CALL lbc_lnk( 'trabbl', zvgx,'V',1.0_wp) + CALL lbc_lnk( 'trabbl', zvgy,'U',1.0_wp) + CALL lbc_lnk( 'trabbl', zpvelu,'V',1.0_wp) + CALL lbc_lnk( 'trabbl', zpvelv,'U',1.0_wp) + + DO_2D_OVR( 1, 0, 1, 0 ) + zu_bbl(ji,jj)=0.25*(zpvelu(ji,jj)+zpvelu(ji+1,jj)+zpvelu(ji,jj-1)+zpvelu(ji+1,jj-1)) + zv_bbl(ji,jj)=0.25*(zpvelv(ji,jj)+zpvelv(ji-1,jj)+zpvelv(ji,jj+1)+zpvelv(ji-1,jj+1)) + + !!ahu_bbl_0(ji,jj) = ahu_bbl_0(ji,jj) + zx_flag*ABS(zub(ji,jj)) * e1u(ji,jj) * e2_e1u(ji,jj) * e3u_bbl_0(ji,jj) * umask(ji,jj,1) + !!ahv_bbl_0(ji,jj) = ahv_bbl_0(ji,jj) + zy_flag*ABS(zvb(ji,jj)) * e2v(ji,jj) * e1_e2v(ji,jj) * e3v_bbl_0(ji,jj) * vmask(ji,jj,1) + + ahu_bbl_0(ji,jj) = ahu_bbl_0(ji,jj) + ABS(zu_bbl(ji,jj)) * e1u(ji,jj) * e2_e1u(ji,jj) * e3u_bbl_0(ji,jj) * umask(ji,jj,1) + ahv_bbl_0(ji,jj) = ahv_bbl_0(ji,jj) + ABS(zv_bbl(ji,jj)) * e2v(ji,jj) * e1_e2v(ji,jj) * e3v_bbl_0(ji,jj) * vmask(ji,jj,1) + + END_2D + + CALL lbc_lnk( 'trabbl', ahu_bbl_0,'U',1.0_wp) + CALL lbc_lnk( 'trabbl', ahv_bbl_0,'V',1.0_wp) + !!CALL iom_put( "deltax", deltax ) ! Friction layer ratio to bottom e3 along u + !!CALL iom_put( "deltay", deltay ) ! Friction layer ratio to bottom e3 along v + !!CALL iom_put( "vgx", zvgx ) ! Plume vg along u + !!CALL iom_put( "vgy", zvgy ) ! Plume vg along v + !!CALL iom_put( "dprgu", dprgu ) ! Plume dprgu + !!CALL iom_put( "dprgv", dprgv ) ! Plume dprgv + + DO_2D_OVR( 1, 0, 1, 0 ) ! criteria: grad(rho).grad(h)<0 and grad(rho).grad(h)<0 + ! ! i-direction + za = zab(ji+1,jj,jp_tem) + zab(ji,jj,jp_tem) ! 2*(alpha,beta) at u-point + zb = zab(ji+1,jj,jp_sal) + zab(ji,jj,jp_sal) + ! ! 2*masked bottom density gradient + zgdrho = ( za * ( zts(ji+1,jj,jp_tem) - zts(ji,jj,jp_tem) ) & + - zb * ( zts(ji+1,jj,jp_sal) - zts(ji,jj,jp_sal) ) ) * umask(ji,jj,1) + ! + zsign = SIGN( 0.5_wp, - zgdrho * REAL( mgrhu(ji,jj) ) ) ! sign of i-gradient * i-slope + zsigna= SIGN( 0.5_wp, zu_bbl(ji,jj) * REAL( mgrhu(ji,jj) ) ) ! sign of u * i-slope + ! + ! ! bbl velocity + utr_bbl(ji,jj) = 0.e0 !!( 0.5 + zsigna ) * ( 0.5 - zsign ) * e2u(ji,jj) * e3u_bbl_0(ji,jj) * zu_bbl(ji,jj) + ! + ! ! j-direction + za = zab(ji,jj+1,jp_tem) + zab(ji,jj,jp_tem) ! 2*(alpha,beta) at v-point + zb = zab(ji,jj+1,jp_sal) + zab(ji,jj,jp_sal) + ! ! 2*masked bottom density gradient + zgdrho = ( za * ( zts(ji,jj+1,jp_tem) - zts(ji,jj,jp_tem) ) & + & - zb * ( zts(ji,jj+1,jp_sal) - zts(ji,jj,jp_sal) ) ) * vmask(ji,jj,1) + zsign = SIGN( 0.5_wp, - zgdrho * REAL( mgrhv(ji,jj) ) ) ! sign of j-gradient * j-slope + zsigna= SIGN( 0.5_wp, zv_bbl(ji,jj) * REAL( mgrhv(ji,jj) ) ) ! sign of u * i-slope + ! + ! ! bbl transport + vtr_bbl(ji,jj) = 0.e0 !!( 0.5 + zsigna ) * ( 0.5 - zsign ) * e1v(ji,jj) * e3v_bbl_0(ji,jj) * zv_bbl(ji,jj) + END_2D + ! + + !!GMDDEV + END SELECT + ! + ENDIF + ! + END SUBROUTINE bbl + + + SUBROUTINE tra_bbl_init + !!---------------------------------------------------------------------- + !! *** ROUTINE tra_bbl_init *** + !! + !! ** Purpose : Initialization for the bottom boundary layer scheme. + !! + !! ** Method : Read the nambbl namelist and check the parameters + !! called by nemo_init at the first timestep (kit000) + !!---------------------------------------------------------------------- + INTEGER :: ji, jj ! dummy loop indices + INTEGER :: ii0, ii1, ij0, ij1, ios ! local integer + REAL(wp), DIMENSION(jpi,jpj) :: zmbku, zmbkv ! workspace + !! + NAMELIST/nambbl/ ln_trabbl, nn_bbl_ldf, nn_bbl_adv, rn_ahtbbl, rn_gambbl + !!---------------------------------------------------------------------- + ! + READ ( numnam_ref, nambbl, IOSTAT = ios, ERR = 901) +901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nambbl in reference namelist' ) + ! + READ ( numnam_cfg, nambbl, IOSTAT = ios, ERR = 902 ) +902 IF( ios > 0 ) CALL ctl_nam ( ios , 'nambbl in configuration namelist' ) + IF(lwm) WRITE ( numond, nambbl ) + ! + IF(lwp) THEN !* Parameter control and print + WRITE(numout,*) + WRITE(numout,*) 'tra_bbl_init : bottom boundary layer initialisation' + WRITE(numout,*) '~~~~~~~~~~~~' + WRITE(numout,*) ' Namelist nambbl : set bbl parameters' + WRITE(numout,*) ' bottom boundary layer flag ln_trabbl = ', ln_trabbl + ENDIF + IF( .NOT.ln_trabbl ) RETURN + ! + IF(lwp) THEN + WRITE(numout,*) ' diffusive bbl (=1) or not (=0) nn_bbl_ldf = ', nn_bbl_ldf + WRITE(numout,*) ' advective bbl (=1/2) or not (=0) nn_bbl_adv = ', nn_bbl_adv + WRITE(numout,*) ' diffusive bbl coefficient rn_ahtbbl = ', rn_ahtbbl, ' m2/s' + WRITE(numout,*) ' advective bbl coefficient rn_gambbl = ', rn_gambbl, ' s' + ENDIF + ! + ! ! allocate trabbl arrays + IF( tra_bbl_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'tra_bbl_init : unable to allocate arrays' ) + ! + IF(lwp) THEN + IF( nn_bbl_adv == 1 ) WRITE(numout,*) ' * Advective BBL using upper velocity' + IF( nn_bbl_adv == 2 ) WRITE(numout,*) ' * Advective BBL using velocity = F( delta rho)' + ENDIF + ! + ! !* vertical index of "deep" bottom u- and v-points + DO_2D( 1, 0, 1, 0 ) ! (the "shelf" bottom k-indices are mbku and mbkv) + mbku_d(ji,jj) = MAX( mbkt(ji+1,jj ) , mbkt(ji,jj) ) ! >= 1 as mbkt=1 over land + mbkv_d(ji,jj) = MAX( mbkt(ji ,jj+1) , mbkt(ji,jj) ) + END_2D + ! converte into REAL to use lbc_lnk ; impose a min value of 1 as a zero can be set in lbclnk + zmbku(:,:) = REAL( mbku_d(:,:), wp ) ; zmbkv(:,:) = REAL( mbkv_d(:,:), wp ) + CALL lbc_lnk( 'trabbl', zmbku,'U',1.0_wp, zmbkv,'V',1.0_wp) + mbku_d(:,:) = MAX( INT( zmbku(:,:) ), 1 ) ; mbkv_d(:,:) = MAX( NINT( zmbkv(:,:) ), 1 ) + ! + ! !* sign of grad(H) at u- and v-points; zero if grad(H) = 0 + !!GMDDEV + grhu(:,:) = 0. ; grhv(:,:) = 0. + !!GMDDEV + mgrhu(:,:) = 0 ; mgrhv(:,:) = 0 + DO_2D( 1, 0, 1, 0 ) + IF( gdept_0(ji+1,jj,mbkt(ji+1,jj)) - gdept_0(ji,jj,mbkt(ji,jj)) /= 0._wp ) THEN + mgrhu(ji,jj) = INT( SIGN( 1.0_wp, gdept_0(ji+1,jj,mbkt(ji+1,jj)) - gdept_0(ji,jj,mbkt(ji,jj)) ) ) + !!GMDDEV + !!grhu(ji,jj) = (gdept_0(ji+1,jj,mbkt(ji+1,jj)) - gdept_0(ji,jj,mbkt(ji,jj)))/e1u(ji,jj) + grhu(ji,jj) = (gdept_0(ji+1,jj,mbkt(ji+1,jj)) - e3t_0(ji+1,jj,mbkt(ji+1,jj)) - gdept_0(ji,jj,mbkt(ji,jj)) + e3t_0(ji,jj,mbkt(ji,jj)))/e1u(ji,jj) + !!GMDDEV + ENDIF + ! + IF( gdept_0(ji,jj+1,mbkt(ji,jj+1)) - gdept_0(ji,jj,mbkt(ji,jj)) /= 0._wp ) THEN + mgrhv(ji,jj) = INT( SIGN( 1.0_wp, gdept_0(ji,jj+1,mbkt(ji,jj+1)) - gdept_0(ji,jj,mbkt(ji,jj)) ) ) + !!GMDDEV + !!grhv(ji,jj) = (gdept_0(ji,jj+1,mbkt(ji,jj+1)) - gdept_0(ji,jj,mbkt(ji,jj)))/e2v(ji,jj) + grhv(ji,jj) = (gdept_0(ji,jj+1,mbkt(ji,jj+1)) - e3t_0(ji,jj+1,mbkt(ji,jj+1)) - gdept_0(ji,jj,mbkt(ji,jj)) + e3t_0(ji,jj,mbkt(ji,jj)))/e2v(ji,jj) + !!GMDDEV + ENDIF + END_2D + ! + ! + !!GMDDEV + CALL lbc_lnk( 'trabbl', grhu,'U',1.0_wp) + CALL lbc_lnk( 'trabbl', grhv,'V',1.0_wp) + !!GMDDEV + DO_2D( 1, 0, 1, 0 ) !* bbl thickness at u- (v-) point; minimum of top & bottom e3u_0 (e3v_0) + e3u_bbl_0(ji,jj) = MIN( e3u_0(ji,jj,mbkt(ji+1,jj )), e3u_0(ji,jj,mbkt(ji,jj)) ) + e3v_bbl_0(ji,jj) = MIN( e3v_0(ji,jj,mbkt(ji ,jj+1)), e3v_0(ji,jj,mbkt(ji,jj)) ) + END_2D + CALL lbc_lnk( 'trabbl', e3u_bbl_0, 'U', 1.0_wp , e3v_bbl_0, 'V', 1.0_wp ) ! lateral boundary conditions + ! + ! !* masked diffusive flux coefficients + ahu_bbl_0(:,:) = rn_ahtbbl * e2_e1u(:,:) * e3u_bbl_0(:,:) * umask(:,:,1) + ahv_bbl_0(:,:) = rn_ahtbbl * e1_e2v(:,:) * e3v_bbl_0(:,:) * vmask(:,:,1) + ! + END SUBROUTINE tra_bbl_init + + !!====================================================================== +END MODULE trabbl