MODULE domzgr !!============================================================================== !! *** MODULE domzgr *** !! Ocean initialization : domain initialization !!============================================================================== !! History : OPA ! 1995-12 (G. Madec) Original code : s vertical coordinate !! ! 1997-07 (G. Madec) lbc_lnk call !! ! 1997-04 (J.-O. Beismann) !! 8.5 ! 2002-09 (A. Bozec, G. Madec) F90: Free form and module !! - ! 2002-09 (A. de Miranda) rigid-lid + islands !! NEMO 1.0 ! 2003-08 (G. Madec) F90: Free form and module !! - ! 2005-10 (A. Beckmann) modifications for hybrid s-ccordinates & new stretching function !! 2.0 ! 2006-04 (R. Benshila, G. Madec) add zgr_zco !! 3.0 ! 2008-06 (G. Madec) insertion of domzgr_zps.h90 & conding style !! 3.2 ! 2009-07 (R. Benshila) Suppression of rigid-lid option !! 3.3 ! 2010-11 (G. Madec) add mbk. arrays associated to the deepest ocean level !! 3.4 ! 2012-12 (R. Bourdalle-Badie and G. Reffray) modify C1D case !!---------------------------------------------------------------------- !!---------------------------------------------------------------------- !! dom_zgr : defined the ocean vertical coordinate system !! zgr_bat : bathymetry fields (levels and meters) !! zgr_bat_zoom : modify the bathymetry field if zoom domain !! zgr_bat_ctl : check the bathymetry files !! zgr_bot_level: deepest ocean level for t-, u, and v-points !! zgr_z : reference z-coordinate !! zgr_zco : z-coordinate !! zgr_zps : z-coordinate with partial steps !! zgr_sco : s-coordinate !! fssig : sigma coordinate non-dimensional function !! dfssig : derivative of the sigma coordinate function !!gm (currently missing!) !!--------------------------------------------------------------------- USE oce ! ocean variables USE dom_oce ! ocean domain USE closea ! closed seas USE c1d ! 1D vertical configuration USE in_out_manager ! I/O manager USE iom ! I/O library USE lbclnk ! ocean lateral boundary conditions (or mpp link) USE lib_mpp ! distributed memory computing library USE wrk_nemo ! Memory allocation USE timing ! Timing IMPLICIT NONE PRIVATE PUBLIC dom_zgr ! called by dom_init.F90 ! !!* Namelist namzgr_sco * REAL(wp) :: rn_sbot_min = 300._wp ! minimum depth of s-bottom surface (>0) (m) REAL(wp) :: rn_sbot_max = 5250._wp ! maximum depth of s-bottom surface (= ocean depth) (>0) (m) REAL(wp) :: rn_theta = 6.00_wp ! surface control parameter (0<=rn_theta<=20) REAL(wp) :: rn_thetb = 0.75_wp ! bottom control parameter (0<=rn_thetb<= 1) REAL(wp) :: rn_rmax = 0.15_wp ! maximum cut-off r-value allowed (0 rn_hmin, dim = 1 ) ! from a depth ENDIF zhmin = gdepw_0(ik+1) ! minimum depth = ik+1 w-levels WHERE( bathy(:,:) <= 0._wp ) ; bathy(:,:) = 0._wp ! min=0 over the lands ELSE WHERE ; bathy(:,:) = MAX( zhmin , bathy(:,:) ) ! min=zhmin over the oceans END WHERE IF(lwp) write(numout,*) 'Minimum ocean depth: ', zhmin, ' minimum number of ocean levels : ', ik ENDIF ! ! CALL wrk_dealloc( jpidta, jpjdta, idta ) CALL wrk_dealloc( jpidta, jpjdta, zdta ) ! IF( nn_timing == 1 ) CALL timing_stop('zgr_bat') ! END SUBROUTINE zgr_bat SUBROUTINE zgr_bat_zoom !!---------------------------------------------------------------------- !! *** ROUTINE zgr_bat_zoom *** !! !! ** Purpose : - Close zoom domain boundary if necessary !! - Suppress Med Sea from ORCA R2 and R05 arctic zoom !! !! ** Method : !! !! ** Action : - update mbathy: level bathymetry (in level index) !!---------------------------------------------------------------------- INTEGER :: ii0, ii1, ij0, ij1 ! temporary integers !!---------------------------------------------------------------------- ! IF(lwp) WRITE(numout,*) IF(lwp) WRITE(numout,*) ' zgr_bat_zoom : modify the level bathymetry for zoom domain' IF(lwp) WRITE(numout,*) ' ~~~~~~~~~~~~' ! ! Zoom domain ! =========== ! ! Forced closed boundary if required IF( lzoom_s ) mbathy( : , mj0(jpjzoom):mj1(jpjzoom) ) = 0 IF( lzoom_w ) mbathy( mi0(jpizoom):mi1(jpizoom) , : ) = 0 IF( lzoom_e ) mbathy( mi0(jpiglo+jpizoom-1):mi1(jpiglo+jpizoom-1) , : ) = 0 IF( lzoom_n ) mbathy( : , mj0(jpjglo+jpjzoom-1):mj1(jpjglo+jpjzoom-1) ) = 0 ! ! Configuration specific domain modifications ! (here, ORCA arctic configuration: suppress Med Sea) IF( cp_cfg == "orca" .AND. lzoom_arct ) THEN SELECT CASE ( jp_cfg ) ! ! ======================= CASE ( 2 ) ! ORCA_R2 configuration ! ! ======================= IF(lwp) WRITE(numout,*) ' ORCA R2 arctic zoom: suppress the Med Sea' ii0 = 141 ; ii1 = 162 ! Sea box i,j indices ij0 = 98 ; ij1 = 110 ! ! ======================= CASE ( 05 ) ! ORCA_R05 configuration ! ! ======================= IF(lwp) WRITE(numout,*) ' ORCA R05 arctic zoom: suppress the Med Sea' ii0 = 563 ; ii1 = 642 ! zero over the Med Sea boxe ij0 = 314 ; ij1 = 370 END SELECT ! mbathy( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0 ! zero over the Med Sea boxe ! ENDIF ! END SUBROUTINE zgr_bat_zoom SUBROUTINE zgr_bat_ctl !!---------------------------------------------------------------------- !! *** ROUTINE zgr_bat_ctl *** !! !! ** Purpose : check the bathymetry in levels !! !! ** Method : The array mbathy is checked to verified its consistency !! with the model options. in particular: !! mbathy must have at least 1 land grid-points (mbathy<=0) !! along closed boundary. !! mbathy must be cyclic IF jperio=1. !! mbathy must be lower or equal to jpk-1. !! isolated ocean grid points are suppressed from mbathy !! since they are only connected to remaining !! ocean through vertical diffusion. !! C A U T I O N : mbathy will be modified during the initializa- !! tion phase to become the number of non-zero w-levels of a water !! column, with a minimum value of 1. !! !! ** Action : - update mbathy: level bathymetry (in level index) !! - update bathy : meter bathymetry (in meters) !!---------------------------------------------------------------------- !! INTEGER :: ji, jj, jl ! dummy loop indices INTEGER :: icompt, ibtest, ikmax ! temporary integers REAL(wp), POINTER, DIMENSION(:,:) :: zbathy !!---------------------------------------------------------------------- ! IF( nn_timing == 1 ) CALL timing_start('zgr_bat_ctl') ! CALL wrk_alloc( jpi, jpj, zbathy ) ! IF(lwp) WRITE(numout,*) IF(lwp) WRITE(numout,*) ' zgr_bat_ctl : check the bathymetry' IF(lwp) WRITE(numout,*) ' ~~~~~~~~~~~' ! ! Suppress isolated ocean grid points IF(lwp) WRITE(numout,*) IF(lwp) WRITE(numout,*)' suppress isolated ocean grid points' IF(lwp) WRITE(numout,*)' -----------------------------------' icompt = 0 DO jl = 1, 2 IF( nperio == 1 .OR. nperio == 4 .OR. nperio == 6 ) THEN mbathy( 1 ,:) = mbathy(jpim1,:) ! local domain is cyclic east-west mbathy(jpi,:) = mbathy( 2 ,:) ENDIF DO jj = 2, jpjm1 DO ji = 2, jpim1 ibtest = MAX( mbathy(ji-1,jj), mbathy(ji+1,jj), & & mbathy(ji,jj-1), mbathy(ji,jj+1) ) IF( ibtest < mbathy(ji,jj) ) THEN IF(lwp) WRITE(numout,*) ' the number of ocean level at ', & & 'grid-point (i,j) = ',ji,jj,' is changed from ', mbathy(ji,jj),' to ', ibtest mbathy(ji,jj) = ibtest icompt = icompt + 1 ENDIF END DO END DO END DO IF( lk_mpp ) CALL mpp_sum( icompt ) IF( icompt == 0 ) THEN IF(lwp) WRITE(numout,*)' no isolated ocean grid points' ELSE IF(lwp) WRITE(numout,*)' ',icompt,' ocean grid points suppressed' ENDIF IF( lk_mpp ) THEN zbathy(:,:) = FLOAT( mbathy(:,:) ) CALL lbc_lnk( zbathy, 'T', 1._wp ) mbathy(:,:) = INT( zbathy(:,:) ) ENDIF ! ! East-west cyclic boundary conditions IF( nperio == 0 ) THEN IF(lwp) WRITE(numout,*) ' mbathy set to 0 along east and west boundary: nperio = ', nperio IF( lk_mpp ) THEN IF( nbondi == -1 .OR. nbondi == 2 ) THEN IF( jperio /= 1 ) mbathy(1,:) = 0 ENDIF IF( nbondi == 1 .OR. nbondi == 2 ) THEN IF( jperio /= 1 ) mbathy(nlci,:) = 0 ENDIF ELSE IF( ln_zco .OR. ln_zps ) THEN mbathy( 1 ,:) = 0 mbathy(jpi,:) = 0 ELSE mbathy( 1 ,:) = jpkm1 mbathy(jpi,:) = jpkm1 ENDIF ENDIF ELSEIF( nperio == 1 .OR. nperio == 4 .OR. nperio == 6 ) THEN IF(lwp) WRITE(numout,*)' east-west cyclic boundary conditions on mbathy: nperio = ', nperio mbathy( 1 ,:) = mbathy(jpim1,:) mbathy(jpi,:) = mbathy( 2 ,:) ELSEIF( nperio == 2 ) THEN IF(lwp) WRITE(numout,*) ' equatorial boundary conditions on mbathy: nperio = ', nperio ELSE IF(lwp) WRITE(numout,*) ' e r r o r' IF(lwp) WRITE(numout,*) ' parameter , nperio = ', nperio ! STOP 'dom_mba' ENDIF ! Boundary condition on mbathy IF( .NOT.lk_mpp ) THEN !!gm !!bug ??? think about it ! ! ... mono- or macro-tasking: T-point, >0, 2D array, no slab zbathy(:,:) = FLOAT( mbathy(:,:) ) CALL lbc_lnk( zbathy, 'T', 1._wp ) mbathy(:,:) = INT( zbathy(:,:) ) ENDIF ! Number of ocean level inferior or equal to jpkm1 ikmax = 0 DO jj = 1, jpj DO ji = 1, jpi ikmax = MAX( ikmax, mbathy(ji,jj) ) END DO END DO !!gm !!! test to do: ikmax = MAX( mbathy(:,:) ) ??? IF( ikmax > jpkm1 ) THEN IF(lwp) WRITE(numout,*) ' maximum number of ocean level = ', ikmax,' > jpk-1' IF(lwp) WRITE(numout,*) ' change jpk to ',ikmax+1,' to use the exact ead bathymetry' ELSE IF( ikmax < jpkm1 ) THEN IF(lwp) WRITE(numout,*) ' maximum number of ocean level = ', ikmax,' < jpk-1' IF(lwp) WRITE(numout,*) ' you can decrease jpk to ', ikmax+1 ENDIF IF( lwp .AND. nprint == 1 ) THEN ! control print WRITE(numout,*) WRITE(numout,*) ' bathymetric field : number of non-zero T-levels ' WRITE(numout,*) ' ------------------' CALL prihin( mbathy, jpi, jpj, 1, jpi, 1, 1, jpj, 1, 3, numout ) WRITE(numout,*) ENDIF ! CALL wrk_dealloc( jpi, jpj, zbathy ) ! IF( nn_timing == 1 ) CALL timing_stop('zgr_bat_ctl') ! END SUBROUTINE zgr_bat_ctl SUBROUTINE zgr_bot_level !!---------------------------------------------------------------------- !! *** ROUTINE zgr_bot_level *** !! !! ** Purpose : defines the vertical index of ocean bottom (mbk. arrays) !! !! ** Method : computes from mbathy with a minimum value of 1 over land !! !! ** Action : mbkt, mbku, mbkv : vertical indices of the deeptest !! ocean level at t-, u- & v-points !! (min value = 1 over land) !!---------------------------------------------------------------------- !! INTEGER :: ji, jj ! dummy loop indices REAL(wp), POINTER, DIMENSION(:,:) :: zmbk !!---------------------------------------------------------------------- ! IF( nn_timing == 1 ) CALL timing_start('zgr_bot_level') ! CALL wrk_alloc( jpi, jpj, zmbk ) ! IF(lwp) WRITE(numout,*) IF(lwp) WRITE(numout,*) ' zgr_bot_level : ocean bottom k-index of T-, U-, V- and W-levels ' IF(lwp) WRITE(numout,*) ' ~~~~~~~~~~~~~' ! mbkt(:,:) = MAX( mbathy(:,:) , 1 ) ! bottom k-index of T-level (=1 over land) ! ! bottom k-index of W-level = mbkt+1 DO jj = 1, jpjm1 ! bottom k-index of u- (v-) level DO ji = 1, jpim1 mbku(ji,jj) = MIN( mbkt(ji+1,jj ) , mbkt(ji,jj) ) mbkv(ji,jj) = MIN( mbkt(ji ,jj+1) , mbkt(ji,jj) ) END DO END DO ! converte into REAL to use lbc_lnk ; impose a min value of 1 as a zero can be set in lbclnk zmbk(:,:) = REAL( mbku(:,:), wp ) ; CALL lbc_lnk(zmbk,'U',1.) ; mbku (:,:) = MAX( INT( zmbk(:,:) ), 1 ) zmbk(:,:) = REAL( mbkv(:,:), wp ) ; CALL lbc_lnk(zmbk,'V',1.) ; mbkv (:,:) = MAX( INT( zmbk(:,:) ), 1 ) ! CALL wrk_dealloc( jpi, jpj, zmbk ) ! IF( nn_timing == 1 ) CALL timing_stop('zgr_bot_level') ! END SUBROUTINE zgr_bot_level SUBROUTINE zgr_zco !!---------------------------------------------------------------------- !! *** ROUTINE zgr_zco *** !! !! ** Purpose : define the z-coordinate system !! !! ** Method : set 3D coord. arrays to reference 1D array !!---------------------------------------------------------------------- INTEGER :: jk !!---------------------------------------------------------------------- ! IF( nn_timing == 1 ) CALL timing_start('zgr_zco') ! DO jk = 1, jpk gdept(:,:,jk) = gdept_0(jk) gdepw(:,:,jk) = gdepw_0(jk) gdep3w(:,:,jk) = gdepw_0(jk) e3t (:,:,jk) = e3t_0(jk) e3u (:,:,jk) = e3t_0(jk) e3v (:,:,jk) = e3t_0(jk) e3f (:,:,jk) = e3t_0(jk) e3w (:,:,jk) = e3w_0(jk) e3uw(:,:,jk) = e3w_0(jk) e3vw(:,:,jk) = e3w_0(jk) END DO ! IF( nn_timing == 1 ) CALL timing_stop('zgr_zco') ! END SUBROUTINE zgr_zco SUBROUTINE zgr_zps !!---------------------------------------------------------------------- !! *** ROUTINE zgr_zps *** !! !! ** Purpose : the depth and vertical scale factor in partial step !! z-coordinate case !! !! ** Method : Partial steps : computes the 3D vertical scale factors !! of T-, U-, V-, W-, UW-, VW and F-points that are associated with !! a partial step representation of bottom topography. !! !! The reference depth of model levels is defined from an analytical !! function the derivative of which gives the reference vertical !! scale factors. !! From depth and scale factors reference, we compute there new value !! with partial steps on 3d arrays ( i, j, k ). !! !! w-level: gdepw(i,j,k) = fsdep(k) !! e3w(i,j,k) = dk(fsdep)(k) = fse3(i,j,k) !! t-level: gdept(i,j,k) = fsdep(k+0.5) !! e3t(i,j,k) = dk(fsdep)(k+0.5) = fse3(i,j,k+0.5) !! !! With the help of the bathymetric file ( bathymetry_depth_ORCA_R2.nc), !! we find the mbathy index of the depth at each grid point. !! This leads us to three cases: !! !! - bathy = 0 => mbathy = 0 !! - 1 < mbathy < jpkm1 !! - bathy > gdepw(jpk) => mbathy = jpkm1 !! !! Then, for each case, we find the new depth at t- and w- levels !! and the new vertical scale factors at t-, u-, v-, w-, uw-, vw- !! and f-points. !! !! This routine is given as an example, it must be modified !! following the user s desiderata. nevertheless, the output as !! well as the way to compute the model levels and scale factors !! must be respected in order to insure second order accuracy !! schemes. !! !! c a u t i o n : gdept_0, gdepw_0 and e3._0 are positives !! - - - - - - - gdept, gdepw and e3. are positives !! !! Reference : Pacanowsky & Gnanadesikan 1997, Mon. Wea. Rev., 126, 3248-3270. !!---------------------------------------------------------------------- !! INTEGER :: ji, jj, jk ! dummy loop indices INTEGER :: ik, it ! temporary integers LOGICAL :: ll_print ! Allow control print for debugging REAL(wp) :: ze3tp , ze3wp ! Last ocean level thickness at T- and W-points REAL(wp) :: zdepwp, zdepth ! Ajusted ocean depth to avoid too small e3t REAL(wp) :: zmax ! Maximum depth REAL(wp) :: zdiff ! temporary scalar REAL(wp) :: zrefdep ! temporary scalar REAL(wp), POINTER, DIMENSION(:,:,:) :: zprt !!--------------------------------------------------------------------- ! IF( nn_timing == 1 ) CALL timing_start('zgr_zps') ! CALL wrk_alloc( jpi, jpj, jpk, zprt ) ! IF(lwp) WRITE(numout,*) IF(lwp) WRITE(numout,*) ' zgr_zps : z-coordinate with partial steps' IF(lwp) WRITE(numout,*) ' ~~~~~~~ ' IF(lwp) WRITE(numout,*) ' mbathy is recomputed : bathy_level file is NOT used' ll_print = .FALSE. ! Local variable for debugging IF(lwp .AND. ll_print) THEN ! control print of the ocean depth WRITE(numout,*) WRITE(numout,*) 'dom_zgr_zps: bathy (in hundred of meters)' CALL prihre( bathy, jpi, jpj, 1,jpi, 1, 1, jpj, 1, 1.e-2, numout ) ENDIF ! bathymetry in level (from bathy_meter) ! =================== zmax = gdepw_0(jpk) + e3t_0(jpk) ! maximum depth (i.e. the last ocean level thickness <= 2*e3t_0(jpkm1) ) bathy(:,:) = MIN( zmax , bathy(:,:) ) ! bounded value of bathy (min already set at the end of zgr_bat) WHERE( bathy(:,:) == 0._wp ) ; mbathy(:,:) = 0 ! land : set mbathy to 0 ELSE WHERE ; mbathy(:,:) = jpkm1 ! ocean : initialize mbathy to the max ocean level END WHERE ! Compute mbathy for ocean points (i.e. the number of ocean levels) ! find the number of ocean levels such that the last level thickness ! is larger than the minimum of e3zps_min and e3zps_rat * e3t_0 (where ! e3t_0 is the reference level thickness DO jk = jpkm1, 1, -1 zdepth = gdepw_0(jk) + MIN( e3zps_min, e3t_0(jk)*e3zps_rat ) WHERE( 0._wp < bathy(:,:) .AND. bathy(:,:) <= zdepth ) mbathy(:,:) = jk-1 END DO ! Scale factors and depth at T- and W-points DO jk = 1, jpk ! intitialization to the reference z-coordinate gdept(:,:,jk) = gdept_0(jk) gdepw(:,:,jk) = gdepw_0(jk) e3t (:,:,jk) = e3t_0 (jk) e3w (:,:,jk) = e3w_0 (jk) END DO ! DO jj = 1, jpj DO ji = 1, jpi ik = mbathy(ji,jj) IF( ik > 0 ) THEN ! ocean point only ! max ocean level case IF( ik == jpkm1 ) THEN zdepwp = bathy(ji,jj) ze3tp = bathy(ji,jj) - gdepw_0(ik) ze3wp = 0.5_wp * e3w_0(ik) * ( 1._wp + ( ze3tp/e3t_0(ik) ) ) e3t(ji,jj,ik ) = ze3tp e3t(ji,jj,ik+1) = ze3tp e3w(ji,jj,ik ) = ze3wp e3w(ji,jj,ik+1) = ze3tp gdepw(ji,jj,ik+1) = zdepwp gdept(ji,jj,ik ) = gdept_0(ik-1) + ze3wp gdept(ji,jj,ik+1) = gdept(ji,jj,ik) + ze3tp ! ELSE ! standard case IF( bathy(ji,jj) <= gdepw_0(ik+1) ) THEN ; gdepw(ji,jj,ik+1) = bathy(ji,jj) ELSE ; gdepw(ji,jj,ik+1) = gdepw_0(ik+1) ENDIF !gm Bug? check the gdepw_0 ! ... on ik gdept(ji,jj,ik) = gdepw_0(ik) + ( gdepw (ji,jj,ik+1) - gdepw_0(ik) ) & & * ((gdept_0( ik ) - gdepw_0(ik) ) & & / ( gdepw_0( ik+1) - gdepw_0(ik) )) e3t (ji,jj,ik) = e3t_0 (ik) * ( gdepw (ji,jj,ik+1) - gdepw_0(ik) ) & & / ( gdepw_0( ik+1) - gdepw_0(ik) ) e3w (ji,jj,ik) = 0.5_wp * ( gdepw(ji,jj,ik+1) + gdepw_0(ik+1) - 2._wp * gdepw_0(ik) ) & & * ( e3w_0(ik) / ( gdepw_0(ik+1) - gdepw_0(ik) ) ) ! ... on ik+1 e3w (ji,jj,ik+1) = e3t (ji,jj,ik) e3t (ji,jj,ik+1) = e3t (ji,jj,ik) gdept(ji,jj,ik+1) = gdept(ji,jj,ik) + e3t(ji,jj,ik) ENDIF ENDIF END DO END DO ! it = 0 DO jj = 1, jpj DO ji = 1, jpi ik = mbathy(ji,jj) IF( ik > 0 ) THEN ! ocean point only e3tp (ji,jj) = e3t(ji,jj,ik ) e3wp (ji,jj) = e3w(ji,jj,ik ) ! test zdiff= gdepw(ji,jj,ik+1) - gdept(ji,jj,ik ) IF( zdiff <= 0._wp .AND. lwp ) THEN it = it + 1 WRITE(numout,*) ' it = ', it, ' ik = ', ik, ' (i,j) = ', ji, jj WRITE(numout,*) ' bathy = ', bathy(ji,jj) WRITE(numout,*) ' gdept = ', gdept(ji,jj,ik), ' gdepw = ', gdepw(ji,jj,ik+1), ' zdiff = ', zdiff WRITE(numout,*) ' e3tp = ', e3t (ji,jj,ik), ' e3wp = ', e3w (ji,jj,ik ) ENDIF ENDIF END DO END DO ! Scale factors and depth at U-, V-, UW and VW-points DO jk = 1, jpk ! initialisation to z-scale factors e3u (:,:,jk) = e3t_0(jk) e3v (:,:,jk) = e3t_0(jk) e3uw(:,:,jk) = e3w_0(jk) e3vw(:,:,jk) = e3w_0(jk) END DO DO jk = 1,jpk ! Computed as the minimum of neighbooring scale factors DO jj = 1, jpjm1 DO ji = 1, fs_jpim1 ! vector opt. e3u (ji,jj,jk) = MIN( e3t(ji,jj,jk), e3t(ji+1,jj,jk) ) e3v (ji,jj,jk) = MIN( e3t(ji,jj,jk), e3t(ji,jj+1,jk) ) e3uw(ji,jj,jk) = MIN( e3w(ji,jj,jk), e3w(ji+1,jj,jk) ) e3vw(ji,jj,jk) = MIN( e3w(ji,jj,jk), e3w(ji,jj+1,jk) ) END DO END DO END DO CALL lbc_lnk( e3u , 'U', 1._wp ) ; CALL lbc_lnk( e3uw, 'U', 1._wp ) ! lateral boundary conditions CALL lbc_lnk( e3v , 'V', 1._wp ) ; CALL lbc_lnk( e3vw, 'V', 1._wp ) ! DO jk = 1, jpk ! set to z-scale factor if zero (i.e. along closed boundaries) WHERE( e3u (:,:,jk) == 0._wp ) e3u (:,:,jk) = e3t_0(jk) WHERE( e3v (:,:,jk) == 0._wp ) e3v (:,:,jk) = e3t_0(jk) WHERE( e3uw(:,:,jk) == 0._wp ) e3uw(:,:,jk) = e3w_0(jk) WHERE( e3vw(:,:,jk) == 0._wp ) e3vw(:,:,jk) = e3w_0(jk) END DO ! Scale factor at F-point DO jk = 1, jpk ! initialisation to z-scale factors e3f(:,:,jk) = e3t_0(jk) END DO DO jk = 1, jpk ! Computed as the minimum of neighbooring V-scale factors DO jj = 1, jpjm1 DO ji = 1, fs_jpim1 ! vector opt. e3f(ji,jj,jk) = MIN( e3v(ji,jj,jk), e3v(ji+1,jj,jk) ) END DO END DO END DO CALL lbc_lnk( e3f, 'F', 1._wp ) ! Lateral boundary conditions ! DO jk = 1, jpk ! set to z-scale factor if zero (i.e. along closed boundaries) WHERE( e3f(:,:,jk) == 0._wp ) e3f(:,:,jk) = e3t_0(jk) END DO !!gm bug ? : must be a do loop with mj0,mj1 ! e3t(:,mj0(1),:) = e3t(:,mj0(2),:) ! we duplicate factor scales for jj = 1 and jj = 2 e3w(:,mj0(1),:) = e3w(:,mj0(2),:) e3u(:,mj0(1),:) = e3u(:,mj0(2),:) e3v(:,mj0(1),:) = e3v(:,mj0(2),:) e3f(:,mj0(1),:) = e3f(:,mj0(2),:) ! Control of the sign IF( MINVAL( e3t (:,:,:) ) <= 0._wp ) CALL ctl_stop( ' zgr_zps : e r r o r e3t <= 0' ) IF( MINVAL( e3w (:,:,:) ) <= 0._wp ) CALL ctl_stop( ' zgr_zps : e r r o r e3w <= 0' ) IF( MINVAL( gdept(:,:,:) ) < 0._wp ) CALL ctl_stop( ' zgr_zps : e r r o r gdepw < 0' ) IF( MINVAL( gdepw(:,:,:) ) < 0._wp ) CALL ctl_stop( ' zgr_zps : e r r o r gdepw < 0' ) ! Compute gdep3w (vertical sum of e3w) gdep3w(:,:,1) = 0.5_wp * e3w(:,:,1) DO jk = 2, jpk gdep3w(:,:,jk) = gdep3w(:,:,jk-1) + e3w(:,:,jk) END DO ! ! ================= ! IF(lwp .AND. ll_print) THEN ! Control print ! ! ! ================= ! DO jj = 1,jpj DO ji = 1, jpi ik = MAX( mbathy(ji,jj), 1 ) zprt(ji,jj,1) = e3t (ji,jj,ik) zprt(ji,jj,2) = e3w (ji,jj,ik) zprt(ji,jj,3) = e3u (ji,jj,ik) zprt(ji,jj,4) = e3v (ji,jj,ik) zprt(ji,jj,5) = e3f (ji,jj,ik) zprt(ji,jj,6) = gdep3w(ji,jj,ik) END DO END DO WRITE(numout,*) WRITE(numout,*) 'domzgr e3t(mbathy)' ; CALL prihre(zprt(:,:,1),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout) WRITE(numout,*) WRITE(numout,*) 'domzgr e3w(mbathy)' ; CALL prihre(zprt(:,:,2),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout) WRITE(numout,*) WRITE(numout,*) 'domzgr e3u(mbathy)' ; CALL prihre(zprt(:,:,3),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout) WRITE(numout,*) WRITE(numout,*) 'domzgr e3v(mbathy)' ; CALL prihre(zprt(:,:,4),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout) WRITE(numout,*) WRITE(numout,*) 'domzgr e3f(mbathy)' ; CALL prihre(zprt(:,:,5),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout) WRITE(numout,*) WRITE(numout,*) 'domzgr gdep3w(mbathy)' ; CALL prihre(zprt(:,:,6),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout) ENDIF ! CALL wrk_dealloc( jpi, jpj, jpk, zprt ) ! IF( nn_timing == 1 ) CALL timing_stop('zgr_zps') ! END SUBROUTINE zgr_zps FUNCTION fssig( pk ) RESULT( pf ) !!---------------------------------------------------------------------- !! *** ROUTINE eos_init *** !! !! ** Purpose : provide the analytical function in s-coordinate !! !! ** Method : the function provide the non-dimensional position of !! T and W (i.e. between 0 and 1) !! T-points at integer values (between 1 and jpk) !! W-points at integer values - 1/2 (between 0.5 and jpk-0.5) !!---------------------------------------------------------------------- REAL(wp), INTENT(in) :: pk ! continuous "k" coordinate REAL(wp) :: pf ! sigma value !!---------------------------------------------------------------------- ! pf = ( TANH( rn_theta * ( -(pk-0.5_wp) / REAL(jpkm1) + rn_thetb ) ) & & - TANH( rn_thetb * rn_theta ) ) & & * ( COSH( rn_theta ) & & + COSH( rn_theta * ( 2._wp * rn_thetb - 1._wp ) ) ) & & / ( 2._wp * SINH( rn_theta ) ) ! END FUNCTION fssig FUNCTION fssig1( pk1, pbb ) RESULT( pf1 ) !!---------------------------------------------------------------------- !! *** ROUTINE eos_init *** !! !! ** Purpose : provide the Song and Haidvogel version of the analytical function in s-coordinate !! !! ** Method : the function provides the non-dimensional position of !! T and W (i.e. between 0 and 1) !! T-points at integer values (between 1 and jpk) !! W-points at integer values - 1/2 (between 0.5 and jpk-0.5) !!---------------------------------------------------------------------- REAL(wp), INTENT(in) :: pk1 ! continuous "k" coordinate REAL(wp), INTENT(in) :: pbb ! Stretching coefficient REAL(wp) :: pf1 ! sigma value !!---------------------------------------------------------------------- ! IF ( rn_theta == 0._wp ) then ! uniform sigma pf1 = - ( pk1 - 0.5_wp ) / REAL( jpkm1 ) ELSE ! stretched sigma pf1 = ( 1._wp - pbb ) * ( SINH( rn_theta*(-(pk1-0.5_wp)/REAL(jpkm1)) ) ) / SINH( rn_theta ) & & + pbb * ( (TANH( rn_theta*( (-(pk1-0.5_wp)/REAL(jpkm1)) + 0.5_wp) ) - TANH( 0.5_wp * rn_theta ) ) & & / ( 2._wp * TANH( 0.5_wp * rn_theta ) ) ) ENDIF ! END FUNCTION fssig1 FUNCTION fssig2 ( pk, kmax ) RESULT( pf2 ) !!---------------------------------------------------------------------- !! *** ROUTINE eos_init *** !! !! ** Purpose : provide the analytical function in s-coordinate !! !! ** Method : the function provide the non-dimensional position of !! T and W (i.e. between 0 and 1) !! T-points at integer values (between 1 and kmax ) !! W-points at integer values - 1/2 (between 0.5 and kmax-0.5) !! !! Reference : ??? !!---------------------------------------------------------------------- REAL(wp), INTENT(in ) :: pk ! continuous "k" coordinate REAL(wp) :: pf2 ! sigma value INTEGER, INTENT (in) :: kmax ! max of sigma)level !!---------------------------------------------------------------------- ! pf2 = ( TANH( rn_theta * ( -(pk-0.5) / REAL(kmax-1,wp) + rn_thetb ) ) & & - TANH( rn_thetb * rn_theta ) ) & & * ( COSH( rn_theta ) & & + COSH( rn_theta * ( 2._wp * rn_thetb - 1._wp ) ) ) & & / ( 2._wp * SINH( rn_theta ) ) ! END FUNCTION fssig2 FUNCTION fssig3( pk1, pbb ,kmax ) RESULT( pf3 ) !!---------------------------------------------------------------------- !! *** ROUTINE eos_init *** !! !! ** Purpose : provide the Song and Haidvogel version of the analytical function in s-coordinate !! !! ** Method : the function provides the non-dimensional position of !! T and W (i.e. between 0 and 1) !! T-points at integer values (between 1 and jpk) !! W-points at integer values - 1/2 (between 0.5 and jpk-0.5) !! !! Reference : ??? !!---------------------------------------------------------------------- REAL(wp), INTENT(in ) :: pk1 ! continuous "k" coordinate REAL(wp), INTENT(in ) :: pbb ! Stretching coefficient REAL(wp) :: pf3 ! sigma value INTEGER, INTENT (in) :: kmax ! max number of s-sigma levels !!---------------------------------------------------------------------- ! IF ( rn_theta == 0 ) then ! uniform sigma pf3 = -(pk1-0.5_wp) / REAL( kmax-1,wp ) ELSE ! stretched sigma pf3 = (1.0-pbb) * (sinh( rn_theta*(-(pk1-0.5_wp)/REAL(kmax-1,wp)) ) ) / sinh(rn_theta) + & & pbb * ( (tanh( rn_theta*( (-(pk1-0.5_wp)/REAL(kmax-1,wp)) + 0.5_wp) ) - tanh(0.5*rn_theta) ) / & & (2._wp*tanh(0.5_wp*rn_theta) ) ) ENDIF END FUNCTION fssig3 SUBROUTINE fszref (zkth, zdzmin, zacr, zhmax,jpup,zhsigm ) INTEGER :: jk ! dummy loop indices REAL(wp) :: zt, zw ! temporary scalars REAL(wp) :: zsur, za0, za1, zkth ! Values set from parameters in REAL(wp) :: zacr, zdzmin, zhmax, zhmax_r ! read from namelist or par_XXX.h90 REAL(wp) :: zhsigm ! depth of sigma layer INTEGER :: jpup, jpkmax ! the last sigma level and number of z-levels !!---------------------------------------------------------------------- ! compute reference depth leveles ! Set variables from parameters ! ------------------------------ ! zkth = rn_kth ; zacr = rn_acr ! zdzmin = rn_dzmin ; zhmax_r = rn_hmax ! za0, za1, zsur are computed from zdzmin , zhmax, zkth, zacr ! jpkmax= jpk - jpup zhmax_r = zhmax - zhsigm za1 = ( zdzmin - zhmax_r / REAL(jpkmax,wp) ) & & / ( TANH((1-zkth)/zacr) - zacr/REAL(jpkmax,wp) * ( LOG( COSH( (jpkmax + 1 - zkth) / zacr) ) & & - LOG( COSH( ( 1 - zkth) / zacr) ) ) ) za0 = zdzmin - za1 * TANH( (1-zkth) / zacr ) zsur = - za0 - za1 * zacr * LOG( COSH( (1-zkth) / zacr ) ) ! Reference z-coordinate (depth - scale factor at T- and W-points) ! ====================== IF( zkth == 0.e0 ) THEN ! uniform vertical grid za1 = zhmax_r / REAL(jpkmax-1,wp) DO jk = 1, jpkmax+1 zw = REAL( jk,wp ) zt = REAL( jk,wp ) + 0.5_wp gdepw_0(jk+jpup-1 ) = ( zw - 1 ) * za1 + zhsigm gdept_0(jk+jpup-1 ) = ( zt - 1 ) * za1 + zhsigm e3w_0 (jk+jpup-1 ) = za1 e3t_0 (jk+jpup-1 ) = za1 END DO ELSE ! Madec & Imbard 1996 function DO jk = 1, jpkmax+1 zw = REAL( jk,wp) zt = REAL( jk,wp ) + 0.5_wp gdepw_0(jk+jpup-1) = ( zsur + za0 * zw + za1 * zacr * LOG ( COSH( (zw-zkth) / zacr ) ) )+zhsigm gdept_0(jk+jpup-1) = ( zsur + za0 * zt + za1 * zacr * LOG ( COSH( (zt-zkth) / zacr ) ) )+zhsigm e3w_0 (jk+jpup-1) = za0 + za1 * TANH( (zw-zkth) / zacr ) e3t_0 (jk+jpup-1) = za0 + za1 * TANH( (zt-zkth) / zacr ) END DO gdepw_0(jpup) = zhsigm ! force first w-level to be exactly at zhsigm ENDIF IF(lwp) WRITE (numout,*) " max and min z-vertical level",jpkmax+1,jpup END SUBROUTINE fszref SUBROUTINE zgr_sco !!---------------------------------------------------------------------- !! *** ROUTINE zgr_sco *** !! !! ** Purpose : define the s-coordinate system !! !! ** Method : s-coordinate !! The depth of model levels is defined as the product of an !! analytical function by the local bathymetry, while the vertical !! scale factors are defined as the product of the first derivative !! of the analytical function by the bathymetry. !! (this solution save memory as depth and scale factors are not !! 3d fields) !! - Read bathymetry (in meters) at t-point and compute the !! bathymetry at u-, v-, and f-points. !! hbatu = mi( hbatt ) !! hbatv = mj( hbatt ) !! hbatf = mi( mj( hbatt ) ) !! - Compute gsigt, gsigw, esigt, esigw from an analytical !! function and its derivative given as function. !! gsigt(k) = fssig (k ) !! gsigw(k) = fssig (k-0.5) !! esigt(k) = fsdsig(k ) !! esigw(k) = fsdsig(k-0.5) !! This routine is given as an example, it must be modified !! following the user s desiderata. nevertheless, the output as !! well as the way to compute the model levels and scale factors !! must be respected in order to insure second order a!!uracy !! schemes. !! !! Reference : Madec, Lott, Delecluse and Crepon, 1996. JPO, 26, 1393-1408. !!---------------------------------------------------------------------- ! INTEGER :: ji, jj, jk, jl ! dummy loop argument INTEGER :: inum ! temporary logical unit INTEGER :: iip1, ijp1, iim1, ijm1, kdep ! temporary integers REAL(wp) :: zcoeft, zcoefw, zrmax, ztaper, maxzenv ! temporary scalars #if defined key_melange REAL(wp) :: rn_hc_bak ! temporary scalars #endif REAL(wp) :: zrfact ! temporary scalars REAL(wp), POINTER, DIMENSION(:,: ) :: ztmpi1, ztmpi2, ztmpj1, ztmpj2 ! #if defined key_fudge REAL(wp), POINTER, DIMENSION(:,: ) :: zenv, zri, zrj, zhbat, fenv #else REAL(wp), POINTER, DIMENSION(:,: ) :: zenv, zri, zrj, zhbat #endif #if defined key_smsh REAL(wp), POINTER, DIMENSION(:,:,:) :: esigtu3, esigtv3, esigtf3, esigwu3, esigwv3 #else REAL(wp), POINTER, DIMENSION(:,:,:) :: gsigw3, gsigt3, gsi3w3 REAL(wp), POINTER, DIMENSION(:,:,:) :: esigt3, esigw3, esigtu3, esigtv3, esigtf3, esigwu3, esigwv3 #endif NAMELIST/namzgr_sco/ rn_sbot_max, rn_sbot_min, rn_theta, rn_thetb, rn_rmax, ln_s_sigma, rn_bb, rn_hc #if defined key_melange NAMELIST/namzgr_sco/ rn_sbot_max, rn_sbot_min, rn_theta, rn_thetb, rn_rmax, ln_s_sigma, rn_bb, rn_hc, nn_sig_lev #endif !!---------------------------------------------------------------------- ! IF( nn_timing == 1 ) CALL timing_start('zgr_sco') ! CALL wrk_alloc( jpi, jpj, ztmpi1, ztmpi2, ztmpj1, ztmpj2 ) #if defined key_fudge CALL wrk_alloc( jpi, jpj, zenv, zri, zrj, zhbat, fenv ) #else CALL wrk_alloc( jpi, jpj, zenv, zri, zrj, zhbat ) #endif #if defined key_smsh CALL wrk_alloc( jpi, jpj, jpk, esigtu3, esigtv3, esigtf3, esigwu3, esigwv3 ) #else CALL wrk_alloc( jpi, jpj, jpk, gsigw3, gsigt3, gsi3w3 ) CALL wrk_alloc( jpi, jpj, jpk, esigt3, esigw3, esigtu3, esigtv3, esigtf3, esigwu3, esigwv3 ) #endif REWIND( numnam ) ! Read Namelist namzgr_sco : sigma-stretching parameters READ ( numnam, namzgr_sco ) IF(lwp) THEN ! control print WRITE(numout,*) WRITE(numout,*) 'dom:zgr_sco : s-coordinate or hybrid z-s-coordinate' WRITE(numout,*) '~~~~~~~~~~~' WRITE(numout,*) ' Namelist namzgr_sco' WRITE(numout,*) ' sigma-stretching coeffs ' WRITE(numout,*) ' maximum depth of s-bottom surface (>0) rn_sbot_max = ' ,rn_sbot_max WRITE(numout,*) ' minimum depth of s-bottom surface (>0) rn_sbot_min = ' ,rn_sbot_min WRITE(numout,*) ' surface control parameter (0<=rn_theta<=20) rn_theta = ', rn_theta WRITE(numout,*) ' bottom control parameter (0<=rn_thetb<= 1) rn_thetb = ', rn_thetb WRITE(numout,*) ' maximum cut-off r-value allowed rn_rmax = ', rn_rmax WRITE(numout,*) ' Hybrid s-sigma-coordinate ln_s_sigma = ', ln_s_sigma WRITE(numout,*) ' stretching parameter (song and haidvogel) rn_bb = ', rn_bb WRITE(numout,*) ' Critical depth rn_hc = ', rn_hc ENDIF #if defined key_melange CALL zgr_zps ! Partial step z-coordinate ! Scale factors and depth at U-, V-, UW and VW-points DO jk = 1, nn_sig_lev ! initialisation to z-scale factors above ln_s_sigma to remove any zps e3u (:,:,jk) = e3t_0(jk) e3v (:,:,jk) = e3t_0(jk) e3uw(:,:,jk) = e3w_0(jk) e3vw(:,:,jk) = e3w_0(jk) END DO #endif gsigw3 = 0._wp ; gsigt3 = 0._wp ; gsi3w3 = 0._wp esigt3 = 0._wp ; esigw3 = 0._wp esigtu3 = 0._wp ; esigtv3 = 0._wp ; esigtf3 = 0._wp esigwu3 = 0._wp ; esigwv3 = 0._wp hift(:,:) = rn_sbot_min ! set the minimum depth for the s-coordinate hifu(:,:) = rn_sbot_min hifv(:,:) = rn_sbot_min hiff(:,:) = rn_sbot_min ! ! set maximum ocean depth bathy(:,:) = MIN( rn_sbot_max, bathy(:,:) ) DO jj = 1, jpj DO ji = 1, jpi IF( bathy(ji,jj) > 0._wp ) bathy(ji,jj) = MAX( rn_sbot_min, bathy(ji,jj) ) END DO END DO ! ! ============================= ! ! Define the envelop bathymetry (hbatt) ! ! ============================= ! use r-value to create hybrid coordinates zenv(:,:) = bathy(:,:) #if defined key_melange DO jj = 1, jpj DO ji = 1, jpi zenv(ji,jj) = MIN( bathy(ji,jj), gdepw_0(nn_sig_lev + 1) ) ENDDO ENDDO #endif #if defined key_fudge CALL iom_open ( 'fudge.nc', inum ) CALL iom_get ( inum, jpdom_data, 'zenv', fenv ) CALL iom_close( inum ) DO jj = 1, jpj DO ji = 1, jpi zenv(ji,jj) = MAX( zenv(ji,jj), fenv(ji,jj) ) ENDDO ENDDO #endif ! ! set first land point adjacent to a wet cell to sbot_min as this needs to be included in smoothing DO jj = 1, jpj DO ji = 1, jpi IF( bathy(ji,jj) == 0._wp ) THEN iip1 = MIN( ji+1, jpi ) ijp1 = MIN( jj+1, jpj ) iim1 = MAX( ji-1, 1 ) ijm1 = MAX( jj-1, 1 ) IF( (bathy(iip1,jj) + bathy(iim1,jj) + bathy(ji,ijp1) + bathy(ji,ijm1) + & & bathy(iip1,ijp1) + bathy(iim1,ijm1) + bathy(iip1,ijp1) + bathy(iim1,ijm1)) > 0._wp ) THEN zenv(ji,jj) = rn_sbot_min ENDIF ENDIF END DO END DO ! apply lateral boundary condition CAUTION: keep the value when the lbc field is zero CALL lbc_lnk( zenv, 'T', 1._wp, 'no0' ) ! ! smooth the bathymetry (if required) scosrf(:,:) = 0._wp ! ocean surface depth (here zero: no under ice-shelf sea) scobot(:,:) = bathy(:,:) ! ocean bottom depth ! jl = 0 zrmax = 1._wp ! ! set scaling factor used in reducing vertical gradients zrfact = ( 1._wp - rn_rmax ) / ( 1._wp + rn_rmax ) ! ! initialise temporary evelope depth arrays ztmpi1(:,:) = zenv(:,:) ztmpi2(:,:) = zenv(:,:) ztmpj1(:,:) = zenv(:,:) ztmpj2(:,:) = zenv(:,:) ! ! initialise temporary r-value arrays zri(:,:) = 1._wp zrj(:,:) = 1._wp ! ! ================ ! DO WHILE( jl <= 10000 .AND. ( zrmax - rn_rmax ) > 1.e-8_wp ) ! Iterative loop ! ! ! ================ ! jl = jl + 1 zrmax = 0._wp ! we set zrmax from previous r-values (zri and zrj) first ! if set after current r-value calculation (as previously) ! we could exit DO WHILE prematurely before checking r-value ! of current zenv DO jj = 1, nlcj DO ji = 1, nlci zrmax = MAX( zrmax, ABS(zri(ji,jj)), ABS(zrj(ji,jj)) ) END DO END DO zri(:,:) = 0._wp zrj(:,:) = 0._wp DO jj = 1, nlcj DO ji = 1, nlci iip1 = MIN( ji+1, nlci ) ! force zri = 0 on last line (ji=ncli+1 to jpi) ijp1 = MIN( jj+1, nlcj ) ! force zrj = 0 on last raw (jj=nclj+1 to jpj) IF( (zenv(ji,jj) > 0._wp) .AND. (zenv(iip1,jj) > 0._wp)) THEN zri(ji,jj) = ( zenv(iip1,jj ) - zenv(ji,jj) ) / ( zenv(iip1,jj ) + zenv(ji,jj) ) END IF IF( (zenv(ji,jj) > 0._wp) .AND. (zenv(ji,ijp1) > 0._wp)) THEN zrj(ji,jj) = ( zenv(ji ,ijp1) - zenv(ji,jj) ) / ( zenv(ji ,ijp1) + zenv(ji,jj) ) END IF IF( zri(ji,jj) > rn_rmax ) ztmpi1(ji ,jj ) = zenv(iip1,jj ) * zrfact IF( zri(ji,jj) < -rn_rmax ) ztmpi2(iip1,jj ) = zenv(ji ,jj ) * zrfact IF( zrj(ji,jj) > rn_rmax ) ztmpj1(ji ,jj ) = zenv(ji ,ijp1) * zrfact IF( zrj(ji,jj) < -rn_rmax ) ztmpj2(ji ,ijp1) = zenv(ji ,jj ) * zrfact END DO END DO IF( lk_mpp ) CALL mpp_max( zrmax ) ! max over the global domain ! IF(lwp)WRITE(numout,*) 'zgr_sco : iter= ',jl, ' rmax= ', zrmax ! DO jj = 1, nlcj DO ji = 1, nlci zenv(ji,jj) = MAX(zenv(ji,jj), ztmpi1(ji,jj), ztmpi2(ji,jj), ztmpj1(ji,jj), ztmpj2(ji,jj) ) END DO END DO ! apply lateral boundary condition CAUTION: keep the value when the lbc field is zero CALL lbc_lnk( zenv, 'T', 1._wp, 'no0' ) ! ! ================ ! END DO ! End loop ! ! ! ================ ! DO jj = 1, jpj DO ji = 1, jpi zenv(ji,jj) = MAX( zenv(ji,jj), rn_sbot_min ) ! set all points to avoid undefined scale value warnings END DO END DO ! ! Envelope bathymetry saved in hbatt hbatt(:,:) = zenv(:,:) IF( MINVAL( gphit(:,:) ) * MAXVAL( gphit(:,:) ) <= 0._wp ) THEN CALL ctl_warn( ' s-coordinates are tapered in vicinity of the Equator' ) DO jj = 1, jpj DO ji = 1, jpi ztaper = EXP( -(gphit(ji,jj)/8._wp)**2._wp ) hbatt(ji,jj) = rn_sbot_max * ztaper + hbatt(ji,jj) * ( 1._wp - ztaper ) END DO END DO ENDIF ! IF(lwp) THEN ! Control print WRITE(numout,*) WRITE(numout,*) ' domzgr: hbatt field; ocean depth in meters' WRITE(numout,*) CALL prihre( hbatt(1,1), jpi, jpj, 1, jpi, 1, 1, jpj, 1, 0._wp, numout ) IF( nprint == 1 ) THEN WRITE(numout,*) ' bathy MAX ', MAXVAL( bathy(:,:) ), ' MIN ', MINVAL( bathy(:,:) ) WRITE(numout,*) ' hbatt MAX ', MAXVAL( hbatt(:,:) ), ' MIN ', MINVAL( hbatt(:,:) ) ENDIF ENDIF ! ! ============================== ! ! hbatu, hbatv, hbatf fields ! ! ============================== IF(lwp) THEN WRITE(numout,*) WRITE(numout,*) ' zgr_sco: minimum depth of the envelop topography set to : ', rn_sbot_min ENDIF hbatu(:,:) = rn_sbot_min hbatv(:,:) = rn_sbot_min hbatf(:,:) = rn_sbot_min DO jj = 1, jpjm1 DO ji = 1, jpim1 ! NO vector opt. hbatu(ji,jj) = 0.50_wp * ( hbatt(ji ,jj) + hbatt(ji+1,jj ) ) hbatv(ji,jj) = 0.50_wp * ( hbatt(ji ,jj) + hbatt(ji ,jj+1) ) hbatf(ji,jj) = 0.25_wp * ( hbatt(ji ,jj) + hbatt(ji ,jj+1) & & + hbatt(ji+1,jj) + hbatt(ji+1,jj+1) ) END DO END DO ! ! Apply lateral boundary condition !!gm ! CAUTION: retain non zero value in the initial file this should be OK for orca cfg, not for EEL zhbat(:,:) = hbatu(:,:) ; CALL lbc_lnk( hbatu, 'U', 1._wp ) DO jj = 1, jpj DO ji = 1, jpi IF( hbatu(ji,jj) == 0._wp ) THEN IF( zhbat(ji,jj) == 0._wp ) hbatu(ji,jj) = rn_sbot_min IF( zhbat(ji,jj) /= 0._wp ) hbatu(ji,jj) = zhbat(ji,jj) ENDIF END DO END DO zhbat(:,:) = hbatv(:,:) ; CALL lbc_lnk( hbatv, 'V', 1._wp ) DO jj = 1, jpj DO ji = 1, jpi IF( hbatv(ji,jj) == 0._wp ) THEN IF( zhbat(ji,jj) == 0._wp ) hbatv(ji,jj) = rn_sbot_min IF( zhbat(ji,jj) /= 0._wp ) hbatv(ji,jj) = zhbat(ji,jj) ENDIF END DO END DO zhbat(:,:) = hbatf(:,:) ; CALL lbc_lnk( hbatf, 'F', 1._wp ) DO jj = 1, jpj DO ji = 1, jpi IF( hbatf(ji,jj) == 0._wp ) THEN IF( zhbat(ji,jj) == 0._wp ) hbatf(ji,jj) = rn_sbot_min IF( zhbat(ji,jj) /= 0._wp ) hbatf(ji,jj) = zhbat(ji,jj) ENDIF END DO END DO !!bug: key_helsinki a verifer hift(:,:) = MIN( hift(:,:), hbatt(:,:) ) hifu(:,:) = MIN( hifu(:,:), hbatu(:,:) ) hifv(:,:) = MIN( hifv(:,:), hbatv(:,:) ) hiff(:,:) = MIN( hiff(:,:), hbatf(:,:) ) IF( nprint == 1 .AND. lwp ) THEN WRITE(numout,*) ' MAX val hif t ', MAXVAL( hift (:,:) ), ' f ', MAXVAL( hiff (:,:) ), & & ' u ', MAXVAL( hifu (:,:) ), ' v ', MAXVAL( hifv (:,:) ) WRITE(numout,*) ' MIN val hif t ', MINVAL( hift (:,:) ), ' f ', MINVAL( hiff (:,:) ), & & ' u ', MINVAL( hifu (:,:) ), ' v ', MINVAL( hifv (:,:) ) WRITE(numout,*) ' MAX val hbat t ', MAXVAL( hbatt(:,:) ), ' f ', MAXVAL( hbatf(:,:) ), & & ' u ', MAXVAL( hbatu(:,:) ), ' v ', MAXVAL( hbatv(:,:) ) WRITE(numout,*) ' MIN val hbat t ', MINVAL( hbatt(:,:) ), ' f ', MINVAL( hbatf(:,:) ), & & ' u ', MINVAL( hbatu(:,:) ), ' v ', MINVAL( hbatv(:,:) ) ENDIF !! helsinki ! ! ======================= ! ! s-ccordinate fields (gdep., e3.) ! ! ======================= ! ! non-dimensional "sigma" for model level depth at w- and t-levels IF( ln_s_sigma ) THEN ! Song and Haidvogel style stretched sigma for depths ! ! below rn_hc, with uniform sigma in shallower waters DO ji = 1, jpi DO jj = 1, jpj IF( hbatt(ji,jj) > rn_hc ) THEN !deep water, stretched sigma DO jk = 1, jpk #if defined key_melange gsigw3(ji,jj,jk) = gdepw_0(jk)/gdepw_0(nn_sig_lev + 1) gsigt3(ji,jj,jk) = gdept_0(jk)/gdepw_0(nn_sig_lev + 1) #else gsigw3(ji,jj,jk) = -fssig1( REAL(jk,wp)-0.5_wp, rn_bb ) gsigt3(ji,jj,jk) = -fssig1( REAL(jk,wp) , rn_bb ) #endif END DO ELSE ! shallow water, uniform sigma DO jk = 1, jpk #if defined key_melange gsigw3(ji,jj,jk) = REAL(jk-1,wp) / REAL(nn_sig_lev,wp) gsigt3(ji,jj,jk) = ( REAL(jk-1,wp) + 0.5_wp ) / REAL(nn_sig_lev,wp) #else gsigw3(ji,jj,jk) = REAL(jk-1,wp) / REAL(jpk-1,wp) gsigt3(ji,jj,jk) = ( REAL(jk-1,wp) + 0.5_wp ) / REAL(jpk-1,wp) #endif END DO ENDIF IF( nprint == 1 .AND. lwp ) WRITE(numout,*) 'gsigw3 1 jpk ', gsigw3(ji,jj,1), gsigw3(ji,jj,jpk) ! DO jk = 1, jpkm1 esigt3(ji,jj,jk ) = gsigw3(ji,jj,jk+1) - gsigw3(ji,jj,jk) esigw3(ji,jj,jk+1) = gsigt3(ji,jj,jk+1) - gsigt3(ji,jj,jk) END DO esigw3(ji,jj,1 ) = 2._wp * ( gsigt3(ji,jj,1 ) - gsigw3(ji,jj,1 ) ) esigt3(ji,jj,jpk) = 2._wp * ( gsigt3(ji,jj,jpk) - gsigw3(ji,jj,jpk) ) ! ! Coefficients for vertical depth as the sum of e3w scale factors gsi3w3(ji,jj,1) = 0.5_wp * esigw3(ji,jj,1) DO jk = 2, jpk gsi3w3(ji,jj,jk) = gsi3w3(ji,jj,jk-1) + esigw3(ji,jj,jk) END DO ! #if defined key_melange DO jk = 1, nn_sig_lev+1 ! DO jk = 1, jpk IF( bathy(ji,jj) < gdepw_0(nn_sig_lev + 1) ) THEN ! should this be bathy or hbatt? #else DO jk = 1, jpk #endif #if defined key_melange zcoeft = ( REAL(jk,wp) - 0.5_wp ) / REAL(nn_sig_lev,wp) zcoefw = ( REAL(jk,wp) - 1.0_wp ) / REAL(nn_sig_lev,wp) ! zcoeft = ( REAL(MIN(jk,nn_sig_lev),wp) - 0.5_wp ) / REAL(nn_sig_lev-1,wp) ! zcoefw = ( REAL(MIN(jk,nn_sig_lev),wp) - 1.0_wp ) / REAL(nn_sig_lev-1,wp) #else zcoeft = ( REAL(jk,wp) - 0.5_wp ) / REAL(jpkm1,wp) zcoefw = ( REAL(jk,wp) - 1.0_wp ) / REAL(jpkm1,wp) #endif #if defined key_melange rn_hc_bak = rn_hc rn_hc = MIN( MAX ( & & (hbatt(ji,jj)-gdepw_0(nn_sig_lev + 1)) / (1._wp - (gdepw_0(nn_sig_lev + 1)/rn_hc)) & & ,0._wp) ,rn_hc) #endif gdept (ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*gsigt3(ji,jj,jk)+rn_hc*zcoeft ) gdepw (ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*gsigw3(ji,jj,jk)+rn_hc*zcoefw ) gdep3w(ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*gsi3w3(ji,jj,jk)+rn_hc*zcoeft ) #if defined key_melange rn_hc = rn_hc_bak #endif IF( gdepw(ji,jj,jk) < 0._wp ) THEN WRITE(*,*) 'zgr_sco : gdepw at point (i,j,k)= ', ji, jj, jk, (gsigw3(ji,jj,jk)*10000._wp-zcoefw*10000._wp) ENDIF #if defined key_melange ENDIF #endif END DO ! END DO ! for all jj's END DO ! for all ji's DO ji = 1, jpim1 DO jj = 1, jpjm1 #if defined key_melange IF( bathy(ji,jj) < gdepw_0(nn_sig_lev + 1) ) THEN ! should this be bathy or hbatt? DO jk = 1, nn_sig_lev+1 ! scale factors should be the same in both zps and sco when H > Hcrit?? ! DO jk = 1, jpk #else DO jk = 1, jpk #endif esigtu3(ji,jj,jk) = ( hbatt(ji,jj)*esigt3(ji,jj,jk)+hbatt(ji+1,jj)*esigt3(ji+1,jj,jk) ) & & / ( hbatt(ji,jj)+hbatt(ji+1,jj) ) esigtv3(ji,jj,jk) = ( hbatt(ji,jj)*esigt3(ji,jj,jk)+hbatt(ji,jj+1)*esigt3(ji,jj+1,jk) ) & & / ( hbatt(ji,jj)+hbatt(ji,jj+1) ) esigtf3(ji,jj,jk) = ( hbatt(ji,jj)*esigt3(ji,jj,jk)+hbatt(ji+1,jj)*esigt3(ji+1,jj,jk) & & + hbatt(ji,jj+1)*esigt3(ji,jj+1,jk)+hbatt(ji+1,jj+1)*esigt3(ji+1,jj+1,jk) ) & & / ( hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1) ) esigwu3(ji,jj,jk) = ( hbatt(ji,jj)*esigw3(ji,jj,jk)+hbatt(ji+1,jj)*esigw3(ji+1,jj,jk) ) & & / ( hbatt(ji,jj)+hbatt(ji+1,jj) ) esigwv3(ji,jj,jk) = ( hbatt(ji,jj)*esigw3(ji,jj,jk)+hbatt(ji,jj+1)*esigw3(ji,jj+1,jk) ) & & / ( hbatt(ji,jj)+hbatt(ji,jj+1) ) ! #if defined key_melange rn_hc_bak = rn_hc rn_hc = MIN( MAX( & & (hbatt(ji,jj)-gdepw_0(nn_sig_lev + 1)) / (1._wp - (gdepw_0(nn_sig_lev + 1)/rn_hc)) & & ,0._wp) ,rn_hc) ! e3t(ji,jj,jk) = ( (hbatt(ji,jj)-rn_hc)*esigt3 (ji,jj,jk) + rn_hc/REAL(nn_sig_lev - 1,wp) ) ! e3w (ji,jj,jk) = ( (hbatt(ji,jj)-rn_hc)*esigw3 (ji,jj,jk) + rn_hc/REAL(nn_sig_lev - 1,wp) ) e3t(ji,jj,jk) = ( (hbatt(ji,jj)-rn_hc)*esigt3 (ji,jj,jk) + rn_hc/REAL(nn_sig_lev ,wp) ) e3w (ji,jj,jk) = ( (hbatt(ji,jj)-rn_hc)*esigw3 (ji,jj,jk) + rn_hc/REAL(nn_sig_lev ,wp) ) #else e3t(ji,jj,jk) = ( (hbatt(ji,jj)-rn_hc)*esigt3 (ji,jj,jk) + rn_hc/REAL(jpkm1,wp) ) e3w (ji,jj,jk) = ( (hbatt(ji,jj)-rn_hc)*esigw3 (ji,jj,jk) + rn_hc/REAL(jpkm1,wp) ) #endif #if defined key_melange rn_hc = MIN( MAX( & & (hbatu(ji,jj)-gdepw_0(nn_sig_lev + 1)) / (1._wp - (gdepw_0(nn_sig_lev + 1)/rn_hc_bak)) & & ,0._wp) ,rn_hc_bak) ! e3u(ji,jj,jk) = ( (hbatu(ji,jj)-rn_hc)*esigtu3(ji,jj,jk) + rn_hc/REAL(nn_sig_lev - 1,wp) ) ! e3uw(ji,jj,jk) = ( (hbatu(ji,jj)-rn_hc)*esigwu3(ji,jj,jk) + rn_hc/REAL(nn_sig_lev - 1,wp) ) e3u(ji,jj,jk) = ( (hbatu(ji,jj)-rn_hc)*esigtu3(ji,jj,jk) + rn_hc/REAL(nn_sig_lev ,wp) ) e3uw(ji,jj,jk) = ( (hbatu(ji,jj)-rn_hc)*esigwu3(ji,jj,jk) + rn_hc/REAL(nn_sig_lev ,wp) ) #else e3u(ji,jj,jk) = ( (hbatu(ji,jj)-rn_hc)*esigtu3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) ) e3uw(ji,jj,jk) = ( (hbatu(ji,jj)-rn_hc)*esigwu3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) ) #endif #if defined key_melange rn_hc = MIN( MAX( & & (hbatv(ji,jj)-gdepw_0(nn_sig_lev + 1)) / (1._wp - (gdepw_0(nn_sig_lev + 1)/rn_hc_bak)) & & ,0._wp) ,rn_hc_bak) ! e3v(ji,jj,jk) = ( (hbatv(ji,jj)-rn_hc)*esigtv3(ji,jj,jk) + rn_hc/REAL(nn_sig_lev - 1,wp) ) ! e3vw(ji,jj,jk) = ( (hbatv(ji,jj)-rn_hc)*esigwv3(ji,jj,jk) + rn_hc/REAL(nn_sig_lev - 1,wp) ) e3v(ji,jj,jk) = ( (hbatv(ji,jj)-rn_hc)*esigtv3(ji,jj,jk) + rn_hc/REAL(nn_sig_lev ,wp) ) e3vw(ji,jj,jk) = ( (hbatv(ji,jj)-rn_hc)*esigwv3(ji,jj,jk) + rn_hc/REAL(nn_sig_lev ,wp) ) #else e3v(ji,jj,jk) = ( (hbatv(ji,jj)-rn_hc)*esigtv3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) ) e3vw(ji,jj,jk) = ( (hbatv(ji,jj)-rn_hc)*esigwv3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) ) #endif #if defined key_melange rn_hc = MIN( MAX( & & (hbatf(ji,jj)-gdepw_0(nn_sig_lev + 1)) / (1._wp - (gdepw_0(nn_sig_lev + 1)/rn_hc_bak)) & & ,0._wp), rn_hc_bak) ! e3f(ji,jj,jk) = ( (hbatf(ji,jj)-rn_hc)*esigtf3(ji,jj,jk) + rn_hc/REAL(nn_sig_lev - 1,wp) ) e3f(ji,jj,jk) = ( (hbatf(ji,jj)-rn_hc)*esigtf3(ji,jj,jk) + rn_hc/REAL(nn_sig_lev ,wp) ) #else e3f(ji,jj,jk) = ( (hbatf(ji,jj)-rn_hc)*esigtf3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp)) #endif ! #if defined key_melange rn_hc = rn_hc_bak #endif END DO #if defined key_melange ENDIF #endif END DO END DO CALL lbc_lnk( e3t , 'T', 1._wp ) CALL lbc_lnk( e3u , 'U', 1._wp ) CALL lbc_lnk( e3v , 'V', 1._wp ) CALL lbc_lnk( e3f , 'F', 1._wp ) CALL lbc_lnk( e3w , 'W', 1._wp ) CALL lbc_lnk( e3uw, 'U', 1._wp ) CALL lbc_lnk( e3vw, 'V', 1._wp ) ! ELSE ! not ln_s_sigma ! DO jk = 1, jpk gsigw(jk) = -fssig( REAL(jk,wp)-0.5_wp ) gsigt(jk) = -fssig( REAL(jk,wp) ) END DO IF( nprint == 1 .AND. lwp ) WRITE(numout,*) 'gsigw 1 jpk ', gsigw(1), gsigw(jpk) ! ! Coefficients for vertical scale factors at w-, t- levels !!gm bug : define it from analytical function, not like juste bellow.... !!gm or betteroffer the 2 possibilities.... DO jk = 1, jpkm1 esigt(jk ) = gsigw(jk+1) - gsigw(jk) esigw(jk+1) = gsigt(jk+1) - gsigt(jk) END DO esigw( 1 ) = 2._wp * ( gsigt(1 ) - gsigw(1 ) ) esigt(jpk) = 2._wp * ( gsigt(jpk) - gsigw(jpk) ) !!gm original form !!org DO jk = 1, jpk !!org esigt(jk)=fsdsig( FLOAT(jk) ) !!org esigw(jk)=fsdsig( FLOAT(jk)-0.5 ) !!org END DO !!gm ! ! Coefficients for vertical depth as the sum of e3w scale factors gsi3w(1) = 0.5_wp * esigw(1) DO jk = 2, jpk gsi3w(jk) = gsi3w(jk-1) + esigw(jk) END DO !!gm: depuw, depvw can be suppressed (modif in ldfslp) and depw=dep3w can be set (save 3 3D arrays) DO jk = 1, jpk zcoeft = ( REAL(jk,wp) - 0.5_wp ) / REAL(jpkm1,wp) zcoefw = ( REAL(jk,wp) - 1.0_wp ) / REAL(jpkm1,wp) gdept (:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*gsigt(jk) + hift(:,:)*zcoeft ) gdepw (:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*gsigw(jk) + hift(:,:)*zcoefw ) gdep3w(:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*gsi3w(jk) + hift(:,:)*zcoeft ) END DO !!gm: e3uw, e3vw can be suppressed (modif in dynzdf, dynzdf_iso, zdfbfr) (save 2 3D arrays) DO jj = 1, jpj DO ji = 1, jpi DO jk = 1, jpk e3t(ji,jj,jk) = ( (hbatt(ji,jj)-hift(ji,jj))*esigt(jk) + hift(ji,jj)/REAL(jpkm1,wp) ) e3u(ji,jj,jk) = ( (hbatu(ji,jj)-hifu(ji,jj))*esigt(jk) + hifu(ji,jj)/REAL(jpkm1,wp) ) e3v(ji,jj,jk) = ( (hbatv(ji,jj)-hifv(ji,jj))*esigt(jk) + hifv(ji,jj)/REAL(jpkm1,wp) ) e3f(ji,jj,jk) = ( (hbatf(ji,jj)-hiff(ji,jj))*esigt(jk) + hiff(ji,jj)/REAL(jpkm1,wp) ) ! e3w (ji,jj,jk) = ( (hbatt(ji,jj)-hift(ji,jj))*esigw(jk) + hift(ji,jj)/REAL(jpkm1,wp) ) e3uw(ji,jj,jk) = ( (hbatu(ji,jj)-hifu(ji,jj))*esigw(jk) + hifu(ji,jj)/REAL(jpkm1,wp) ) e3vw(ji,jj,jk) = ( (hbatv(ji,jj)-hifv(ji,jj))*esigw(jk) + hifv(ji,jj)/REAL(jpkm1,wp) ) END DO END DO END DO ! ENDIF ! ln_s_sigma where (e3t (:,:,:).eq.0._wp) e3t(:,:,:) = 1._wp where (e3u (:,:,:).eq.0._wp) e3u(:,:,:) = 1._wp where (e3v (:,:,:).eq.0._wp) e3v(:,:,:) = 1._wp where (e3f (:,:,:).eq.0._wp) e3f(:,:,:) = 1._wp where (e3w (:,:,:).eq.0._wp) e3w(:,:,:) = 1._wp where (e3uw (:,:,:).eq.0._wp) e3uw(:,:,:) = 1._wp where (e3vw (:,:,:).eq.0._wp) e3vw(:,:,:) = 1._wp fsdept(:,:,:) = gdept (:,:,:) fsdepw(:,:,:) = gdepw (:,:,:) fsde3w(:,:,:) = gdep3w(:,:,:) fse3t (:,:,:) = e3t (:,:,:) fse3u (:,:,:) = e3u (:,:,:) fse3v (:,:,:) = e3v (:,:,:) fse3f (:,:,:) = e3f (:,:,:) fse3w (:,:,:) = e3w (:,:,:) fse3uw(:,:,:) = e3uw (:,:,:) fse3vw(:,:,:) = e3vw (:,:,:) !! ! HYBRID : DO jj = 1, jpj DO ji = 1, jpi #if defined key_melange IF( bathy(ji,jj) < gdepw_0(nn_sig_lev + 1) ) THEN ! should this be hbatt or bathy DO jk = 1, nn_sig_lev ! DO jk = 1, jpkm1 #else DO jk = 1, jpkm1 #endif IF( scobot(ji,jj) >= fsdept(ji,jj,jk) ) mbathy(ji,jj) = MAX( 2, jk ) END DO #if defined key_melange ENDIF #endif IF( scobot(ji,jj) == 0._wp ) mbathy(ji,jj) = 0 END DO END DO IF( nprint == 1 .AND. lwp ) WRITE(numout,*) ' MIN val mbathy h90 ', MINVAL( mbathy(:,:) ), & & ' MAX ', MAXVAL( mbathy(:,:) ) ! ! ============= IF(lwp) THEN ! Control print ! ! ============= WRITE(numout,*) WRITE(numout,*) ' domzgr: vertical coefficients for model level' WRITE(numout, "(9x,' level gsigt gsigw esigt esigw gsi3w')" ) WRITE(numout, "(10x,i4,5f11.4)" ) ( jk, gsigt(jk), gsigw(jk), esigt(jk), esigw(jk), gsi3w(jk), jk=1,jpk ) ENDIF IF( nprint == 1 .AND. lwp ) THEN ! min max values over the local domain WRITE(numout,*) ' MIN val mbathy ', MINVAL( mbathy(:,:) ), ' MAX ', MAXVAL( mbathy(:,:) ) WRITE(numout,*) ' MIN val depth t ', MINVAL( fsdept(:,:,:) ), & & ' w ', MINVAL( fsdepw(:,:,:) ), '3w ' , MINVAL( fsde3w(:,:,:) ) WRITE(numout,*) ' MIN val e3 t ', MINVAL( fse3t (:,:,:) ), ' f ' , MINVAL( fse3f (:,:,:) ), & & ' u ', MINVAL( fse3u (:,:,:) ), ' u ' , MINVAL( fse3v (:,:,:) ), & & ' uw', MINVAL( fse3uw(:,:,:) ), ' vw' , MINVAL( fse3vw(:,:,:) ), & & ' w ', MINVAL( fse3w (:,:,:) ) WRITE(numout,*) ' MAX val depth t ', MAXVAL( fsdept(:,:,:) ), & & ' w ', MAXVAL( fsdepw(:,:,:) ), '3w ' , MAXVAL( fsde3w(:,:,:) ) WRITE(numout,*) ' MAX val e3 t ', MAXVAL( fse3t (:,:,:) ), ' f ' , MAXVAL( fse3f (:,:,:) ), & & ' u ', MAXVAL( fse3u (:,:,:) ), ' u ' , MAXVAL( fse3v (:,:,:) ), & & ' uw', MAXVAL( fse3uw(:,:,:) ), ' vw' , MAXVAL( fse3vw(:,:,:) ), & & ' w ', MAXVAL( fse3w (:,:,:) ) ENDIF ! IF(lwp) THEN ! selected vertical profiles WRITE(numout,*) WRITE(numout,*) ' domzgr: vertical coordinates : point (1,1,k) bathy = ', bathy(1,1), hbatt(1,1) WRITE(numout,*) ' ~~~~~~ --------------------' WRITE(numout,"(9x,' level gdept gdepw gde3w e3t e3w ')") WRITE(numout,"(10x,i4,4f9.2)") ( jk, fsdept(1,1,jk), fsdepw(1,1,jk), & & fse3t (1,1,jk), fse3w (1,1,jk), jk=1,jpk ) DO jj = mj0(20), mj1(20) DO ji = mi0(20), mi1(20) WRITE(numout,*) WRITE(numout,*) ' domzgr: vertical coordinates : point (20,20,k) bathy = ', bathy(ji,jj), hbatt(ji,jj) WRITE(numout,*) ' ~~~~~~ --------------------' WRITE(numout,"(9x,' level gdept gdepw gde3w e3t e3w ')") WRITE(numout,"(10x,i4,4f9.2)") ( jk, fsdept(ji,jj,jk), fsdepw(ji,jj,jk), & & fse3t (ji,jj,jk), fse3w (ji,jj,jk), jk=1,jpk ) END DO END DO DO jj = mj0(74), mj1(74) DO ji = mi0(100), mi1(100) WRITE(numout,*) WRITE(numout,*) ' domzgr: vertical coordinates : point (100,74,k) bathy = ', bathy(ji,jj), hbatt(ji,jj) WRITE(numout,*) ' ~~~~~~ --------------------' WRITE(numout,"(9x,' level gdept gdepw gde3w e3t e3w ')") WRITE(numout,"(10x,i4,4f9.2)") ( jk, fsdept(ji,jj,jk), fsdepw(ji,jj,jk), & & fse3t (ji,jj,jk), fse3w (ji,jj,jk), jk=1,jpk ) END DO END DO ENDIF !!gm bug? no more necessary? if ! defined key_helsinki DO jk = 1, jpk DO jj = 1, jpj DO ji = 1, jpi IF( fse3w(ji,jj,jk) <= 0._wp .OR. fse3t(ji,jj,jk) <= 0._wp ) THEN WRITE(*,*) 'zgr_sco : e3w or e3t =< 0 at point (i,j,k)= ', ji, jj, jk, fse3w(ji,jj,jk), fse3t(ji,jj,jk) ! WRITE(ctmp1,*) 'zgr_sco : e3w or e3t =< 0 at point (i,j,k)= ', ji, jj, jk ! CALL ctl_stop( ctmp1 ) ENDIF IF( fsdepw(ji,jj,jk) < 0._wp .OR. fsdept(ji,jj,jk) < 0._wp ) THEN WRITE(*,*) 'zgr_sco : gdepw or gdept =< 0 at point (i,j,k)= ', ji, jj, jk, fsdepw(ji,jj,jk), fsdept(ji,jj,jk) ! WRITE(ctmp1,*) 'zgr_sco : gdepw or gdept =< 0 at point (i,j,k)= ', ji, jj, jk ! CALL ctl_stop( ctmp1 ) ENDIF END DO END DO END DO !!gm bug #endif ! CALL wrk_dealloc( jpi, jpj, zenv, ztmpi1, ztmpi2, ztmpj1, ztmpj2, zri, zrj, zhbat ) #if defined key_smsh CALL wrk_dealloc( jpi, jpj, jpk, esigtu3, esigtv3, esigtf3, esigwu3, esigwv3 ) #else CALL wrk_dealloc( jpi, jpj, jpk, gsigw3, gsigt3, gsi3w3 ) CALL wrk_dealloc( jpi, jpj, jpk, esigt3, esigw3, esigtu3, esigtv3, esigtf3, esigwu3, esigwv3 ) #endif ! IF( nn_timing == 1 ) CALL timing_stop('zgr_sco') ! END SUBROUTINE zgr_sco SUBROUTINE zgr_hyb !!---------------------------------------------------------------------- !! *** ROUTINE zgr_sco *** !! Combination of zgr_sco in upper layers ( shelf ) and zgr_zps in abyss !! !! ** Purpose : define the s-z coordinate system !! !! ** Method : s-coordinate in upper layers and z-coordinates below !! The depth of model levels is defined as the product of an !! analytical function by the local bathymetry, while the vertical !! scale factors are defined as the product of the first derivative !! of the analytical function by the bathymetry. !! (this solution save memory as depth and scale factors are not !! 3d fields) !! - Read bathymetry (in meters) at t-point and compute the !! bathymetry at u-, v-, and f-points. !! hbatu = mi( hbatt ) !! hbatv = mj( hbatt ) !! hbatf = mi( mj( hbatt ) ) !! - Compute gsigt, gsigw, esigt, esigw from an analytical !! function !! gsigt(k) = fssig (k ) !! gsigw(k) = fssig (k-0.5) !! This routine is given as an example, it must be modified !! following the user s desiderata. nevertheless, the output as !! well as the way to compute the model levels and scale factors !! must be respected in order to insure second order a!!uracy !! schemes. !! !!====================================================================== INTEGER :: ji, jj, jk, jl, ik ! dummy loop argument INTEGER :: iip1, ijp1, iim1, ijm1 ! temporary integers INTEGER :: jpksigm ! temporary integer for maxnumber of s-levels REAL(wp) :: zcoeft, zcoefw, zrmax, ztaper,zrmin,e3t_t,e3w_t ! temporary scalars REAL(wp) :: ze3tp , ze3wp ! Last ocean level thickness at T- and W-points REAL(wp) :: zdepwp, zdepth ! Ajusted ocean depth to avoid too small e3t REAL(wp) :: zmax, zmin ,zsigma ! Maximum and minimum depth and depth of sigma layer REAL(wp) :: zacr , zkth ,za1 ! parameters for z- layer (as ppacr , ppzkth ) ! REAL(wp), POINTER, DIMENSION(:,: ) :: zenv, ztmp, zmsk, zri, zrj, zhbat ,zrpt NAMELIST/namzgr_hyb/ rn_sbot_max, rn_sbot_min, rn_theta, rn_thetb, rn_rmax, & ln_s_sigma, rn_bb, rn_hc,rn_zsigma , nn_sig_lev , rn_kth , rn_acr !!---------------------------------------------------------------------- ! IF( nn_timing == 1 ) CALL timing_start('zgr_hyb') ! CALL wrk_alloc( jpi, jpj, zenv, ztmp, zmsk, zri, zrj, zhbat ) ! CALL wrk_alloc( jpi, jpj, jpk, gsigw3, gsigt3 ) ! REWIND( numnam ) ! Read Namelist namzgr_sco : sigma-stretching parameters READ ( numnam, namzgr_hyb ) IF(lwp) THEN ! control print WRITE(numout,*) WRITE(numout,*) 'dom:zgr_hyb : s-coordinate or hybrid z-s-coordinate' WRITE(numout,*) '~~~~~~~~~~~' WRITE(numout,*) ' Namelist namzgr_hyb' WRITE(numout,*) ' sigma-stretching coeffs ' WRITE(numout,*) ' maximum depth of s-bottom surface (>0) rn_sbot_max = ' ,rn_sbot_max WRITE(numout,*) ' minimum depth of s-bottom surface (>0) rn_sbot_min = ' ,rn_sbot_min WRITE(numout,*) ' surface control parameter (0<=rn_theta<=20) rn_theta = ', rn_theta WRITE(numout,*) ' bottom control parameter (0<=rn_thetb<= 1) rn_thetb = ', rn_thetb WRITE(numout,*) ' maximum cut-off r-value allowed rn_rmax = ', rn_rmax WRITE(numout,*) ' Hybrid s-sigma-coordinate ln_s_sigma = ', ln_s_sigma WRITE(numout,*) ' stretching parameter (song and haidvogel) rn_bb = ', rn_bb WRITE(numout,*) ' Critical depth rn_hc = ', rn_hc WRITE(numout,*) ' Sigma depth rn_zsigma = ', rn_zsigma WRITE(numout,*) ' The same as pp_arc rn_arc = ', rn_acr WRITE(numout,*) ' Number of sigma levels rn_arc = ', nn_sig_lev WRITE(numout,*) ' Number of levels for stretching z-coord rn_kth = ', rn_kth ENDIF zsigma = rn_zsigma jpksigm = nn_sig_lev zmax = rn_sbot_max zacr = rn_acr zkth = rn_kth e3t(:,:,:) = 1._wp e3w(:,:,:) = 1._wp e3u(:,:,:) = 1._wp e3v(:,:,:) = 1._wp e3f(:,:,:) = 1._wp e3uw(:,:,:)= 1._wp e3vw(:,:,:)= 1._wp DO jj = 1, jpj DO ji= 1, jpi IF( bathy(ji,jj) <= 0._wp ) THEN ; bathy(ji,jj) = 0.e0_wp ELSE ; bathy(ji,jj) = MIN( rn_sbot_max, MAX( bathy(ji,jj),rn_sbot_min ) ) ENDIF END DO END DO ! create bathymetry for enveloping DO jj = 1, jpj DO ji = 1, jpi zenv(ji,jj) = MAX( bathy(ji,jj), rn_sbot_min ) zenv(ji,jj) = MIN (zenv(ji,jj), zsigma ) hbatt(ji,jj) = zenv(ji,jj) END DO END DO scosrf(:,:) = 0._wp ! ocean surface depth (here zero: no under ice-shelf sea) scobot(:,:) = bathy(:,:) ! ocean bottom depth jl = 0 zrmax = 1._wp ! ! ================ ! DO WHILE( jl <= 10000 .AND. zrmax > rn_rmax ) ! Iterative loop ! ! ! ================ ! jl = jl + 1 zrmax = 0._wp zmsk(:,:) = 0._wp DO jj = 1, nlcj DO ji = 1, nlci iip1 = MIN( ji+1, nlci ) ! force zri = 0 on last line (ji=ncli+1 to jpi) ijp1 = MIN( jj+1, nlcj ) ! force zrj = 0 on last raw (jj=nclj+1 to jpj) zri(ji,jj) = ABS( zenv(iip1,jj ) - zenv(ji,jj) ) / ( zenv(iip1,jj ) + zenv(ji,jj) ) zrj(ji,jj) = ABS( zenv(ji ,ijp1) - zenv(ji,jj) ) / ( zenv(ji ,ijp1) + zenv(ji,jj) ) zrmax = MAX( zrmax, zri(ji,jj), zrj(ji,jj) ) IF( zri(ji,jj) > rn_rmax ) zmsk(ji ,jj ) = 1._wp IF( zri(ji,jj) > rn_rmax ) zmsk(iip1,jj ) = 1._wp IF( zrj(ji,jj) > rn_rmax ) zmsk(ji ,jj ) = 1._wp IF( zrj(ji,jj) > rn_rmax ) zmsk(ji ,ijp1) = 1._wp END DO END DO IF( lk_mpp ) CALL mpp_max( zrmax ) ! max over the global domain ! lateral boundary condition on zmsk: keep 1 along closed boundary (use of MAX) ztmp(:,:) = zmsk(:,:) ; CALL lbc_lnk( zmsk, 'T', 1._wp ) DO jj = 1, nlcj DO ji = 1, nlci zmsk(ji,jj) = MAX( zmsk(ji,jj), ztmp(ji,jj) ) END DO END DO IF(lwp)WRITE(numout,*) 'zgr_hyb : iter= ',jl, ' rmax= ', zrmax, ' nb of pt= ', INT( SUM(zmsk(:,:) ) ) ! DO jj = 1, nlcj DO ji = 1, nlci iip1 = MIN( ji+1, nlci ) ! last line (ji=nlci) ijp1 = MIN( jj+1, nlcj ) ! last raw (jj=nlcj) iim1 = MAX( ji-1, 1 ) ! first line (ji=nlci) ijm1 = MAX( jj-1, 1 ) ! first raw (jj=nlcj) IF( zmsk(ji,jj) == 1._wp ) THEN ztmp(ji,jj) = ( & & zenv(iim1,ijp1)*zmsk(iim1,ijp1) + zenv(ji,ijp1)*zmsk(ji,ijp1) + zenv(iip1,ijp1)*zmsk(iip1,ijp1) & & + zenv(iim1,jj )*zmsk(iim1,jj ) + zenv(ji,jj )* 2._wp + zenv(iip1,jj )*zmsk(iip1,jj ) & & + zenv(iim1,ijm1)*zmsk(iim1,ijm1) + zenv(ji,ijm1)*zmsk(ji,ijm1) + zenv(iip1,ijm1)*zmsk(iip1,ijm1) & & ) / ( & & zmsk(iim1,ijp1) + zmsk(ji,ijp1) + zmsk(iip1,ijp1) & & + zmsk(iim1,jj ) + 2._wp + zmsk(iip1,jj ) & & + zmsk(iim1,ijm1) + zmsk(ji,ijm1) + zmsk(iip1,ijm1) & & ) ENDIF END DO END DO ! DO jj = 1, nlcj DO ji = 1, nlci IF( zmsk(ji,jj) == 1._wp ) zenv(ji,jj) = MAX( ztmp(ji,jj), hbatt(ji,jj) ) END DO END DO ! ! Apply lateral boundary condition CAUTION: kept the value when the lbc field is zero ztmp(:,:) = zenv(:,:) ; CALL lbc_lnk( zenv, 'T', 1._wp ) DO jj = 1, nlcj DO ji = 1, nlci IF( zenv(ji,jj) == 0._wp ) zenv(ji,jj) = ztmp(ji,jj) END DO END DO ! ! ================ ! END DO ! End loop ! ! ! ================ ! ! ! ! envelop bathymetry saved in hbatt hbatt(:,:) = zenv(:,:) IF( lk_mpp ) CALL mpp_max( nstop ) IF (lwp) write(numout,*)"after envelope", nstop ! define new reference levels ! for s- levels, allows stretching at surface and bottom layer IF( jpksigm > 1 )THEN IF(ln_s_sigma)THEN DO jk = 1, jpksigm gsigw(jk) = -fssig3( REAL(jk,wp)-0.5_wp , rn_bb,jpksigm ) gsigt(jk) = -fssig3( REAL(jk,wp) , rn_bb,jpksigm ) END DO ELSE DO jk = 1, jpksigm gsigw(jk) = -fssig2( REAL(jk,wp)-0.5_wp ,jpksigm ) gsigt(jk) = -fssig2( REAL(jk,wp) ,jpksigm ) END DO ENDIF gsigw(1)=0._wp ! set gsigw exactly to zero IF( lk_mpp ) CALL mpp_max( nstop ) IF (lwp) THEN write(numout,*)"after fssig", nstop,"gsigw,gsigt=" do jk=1,jpksigm write(numout,*)jk,gsigw(jk),gsigt(jk) enddo ENDIF DO jk=1,jpksigm DO jj=1,jpj DO ji=1,jpi zrmin= min( hbatt(ji,jj), zsigma ) IF(hbatt(ji,jj).lt.rn_hc)THEN zcoefw=REAL(jk-1,wp) / REAL(jpksigm-1,wp) zcoeft=(REAL(jk-1,wp)+0.5)/ REAL(jpksigm-1,wp) ELSE zcoefw=gsigw(jk) zcoeft=gsigt(jk) ENDIF gdept(ji,jj,jk)=scosrf(ji,jj)+(zrmin-rn_hc)*zcoeft & +rn_hc* (REAL(jk,wp)- 0.5_wp) / REAL(jpksigm-1,wp) gdepw(ji,jj,jk)=scosrf(ji,jj)+(zrmin-rn_hc)*zcoefw & +rn_hc*( REAL(jk,wp)- 1.0_wp ) / REAL(jpksigm-1,wp) ENDDO ENDDO ENDDO gdepw(:,:,1)= scosrf(:,:) ! redefine gdept_0, gdepw_0 which will be used in diawri.F90 DO jk=1,jpksigm IF(zsigma.lt.rn_hc)THEN zcoefw=REAL(jk-1,wp) / REAL(jpksigm-1,wp) zcoeft=(REAL(jk-1,wp)+0.5)/ REAL(jpksigm-1,wp) ELSE zcoefw=gsigw(jk) zcoeft=gsigt(jk) ENDIF gdept_0(jk)= zcoeft * (zsigma-rn_hc)+rn_hc* (REAL(jk,wp)- 0.5_wp )/REAL(jpksigm-1,wp) gdepw_0(jk)= zcoefw * (zsigma-rn_hc)+rn_hc* (REAL(jk,wp)- 1.0_wp )/REAL(jpksigm-1,wp) ENDDO DO jk=1,jpksigm-1 e3t_0(jk) = gdepw_0(jk+1)-gdepw_0(jk) e3w_0(jk+1)= gdept_0(jk+1)-gdept_0(jk) ENDDO e3w_0(1) = 2._wp * ( gdept_0(1 ) - gdepw_0(1 ) ) ! now for lower z-levels : zmin = e3t_0 (jpksigm -1 ) ! min layer width in z- zone is the same as lowest in s- layer ELSE zsigma = 0._wp hbatt(:,:) = zsigma zmin = 5._wp ENDIF IF(lwp) write(numout,*) ": last vertical level of sigma-coordinates",zmin CALL fszref ( zkth, zmin, zacr, zmax, jpksigm , zsigma ) DO jk=jpksigm,jpkm1 e3t_0(jk) = gdepw_0(jk+1)-gdepw_0(jk) e3w_0(jk+1)= gdept_0(jk+1)-gdept_0(jk) ENDDO e3w_0(1) = 2._wp * ( gdept_0(1 ) - gdepw_0(1 ) ) e3t_0(jpk) = 2._wp * ( gdept_0(jpk) - gdepw_0(jpk) ) IF( lk_mpp ) CALL mpp_max( nstop ) IF (lwp) write(numout,*)"e3t0" ,nstop IF(lwp) THEN ! control print WRITE(numout,*) WRITE(numout,*) ' zhyb Reference z-coordinate depth and scale factors:' WRITE(numout, "(9x,' level gdept gdepw e3t e3w ')" ) WRITE(numout, "(10x, i4, 4f9.2)" ) ( jk, gdept_0(jk), gdepw_0(jk), e3t_0(jk), e3w_0(jk), jk = 1, jpk ) open(333,file='zmesh.dat') WRITE(333,*)'initial zsigma =', zsigma, jpk WRITE(333,*)'initial jpksigm =', jpksigm WRITE(333,*)'rn_bb =', rn_bb,'rn_theta=',rn_theta do jk=1,jpk WRITE(333,'(i4,1x,4(1x,e13.6))') jk, gdept_0(jk), gdepw_0(jk), e3t_0(jk),e3w_0(jk) enddo close(333) ENDIF DO jk=jpksigm,jpk DO jj=1,jpj DO ji=1,jpi IF(jpksigm>1)THEN gdept(ji,jj,jk)=scosrf(ji,jj) + gdept_0(jk) *hbatt(ji,jj)/zsigma ! differ from gdept0+scorf only at land gdepw(ji,jj,jk)=scosrf(ji,jj) + gdepw_0(jk) *hbatt(ji,jj)/zsigma ! as hbatt=zsigma over deep part of basin ELSE gdept(ji,jj,jk)=scosrf(ji,jj) + gdept_0(jk) gdepw(ji,jj,jk)=scosrf(ji,jj) + gdepw_0(jk) ENDIF ENDDO ENDDO ENDDO ! define e3t, e3w for general levels DO jk=1,jpkm1 e3t(:,:,jk) = gdepw(:,:,jk+1)-gdepw(:,:,jk) e3w(:,:,jk+1)= gdept(:,:,jk+1)-gdept(:,:,jk) ENDDO e3w(:,:,1) = 2._wp * ( gdept(:,:,1 ) - gdepw(:,:,1 ) ) e3t(:,:,jpk) = 2._wp * ( gdept(:,:,jpk) - gdepw(:,:,jpk) ) ! and surface : ! ! HYBRID mbathy : IF(lwp) THEN ! control print WRITE(numout,*) WRITE(numout,*) ' zhyb centre of basin s-z-coordinate depth and scale factors:' WRITE(numout, "(9x,' level gdept gdepw e3t e3w ')" ) write(numout,*)"bathy" ,"min e3t" do jk=1,jpk WRITE(numout, "(10x, i4, 4f9.2)" ) jk, gdept(20,20,jk), gdepw(20,20,jk), & & e3t(20,20,jk), e3w(20,20,jk) enddo ENDIF mbathy(:,:)=0 ! WHERE( 0._wp < bathy(:,:)) mbathy(:,:)=jpkm1 DO jj = 1, jpj DO ji = 1, jpi DO jk = 1, jpkm1 IF( bathy(ji,jj) >= gdept(ji,jj,jk) ) mbathy(ji,jj) = MAX( 2, jk ) END DO END DO END DO ! DO jk = jpkm1, jpksigm+1, -1 ! zdepth = gdepw_0(jk) + MIN( e3zps_min, e3t_0(jk)*e3zps_rat ) ! WHERE( 0._wp < bathy(:,:) .AND. bathy(:,:) <= zdepth ) mbathy(:,:) = jk-1 ! END DO ! z-partial steps :goto 20 DO jj = 1, jpj DO ji = 1, jpi ik = mbathy(ji,jj) IF( ik > jpksigm ) THEN ! ocean point only ! max ocean level case IF( ik == jpkm1 ) THEN zdepwp = bathy(ji,jj) ze3tp = bathy(ji,jj) - gdepw_0(ik) ze3wp = 0.5_wp * e3w_0(ik) * ( 1._wp + ( ze3tp/e3t_0(ik) ) ) e3t(ji,jj,ik ) = ze3tp e3t(ji,jj,ik+1) = ze3tp e3w(ji,jj,ik ) = ze3wp e3w(ji,jj,ik+1) = ze3tp gdepw(ji,jj,ik+1) = zdepwp gdept(ji,jj,ik ) = gdept_0(ik-1) + ze3wp gdept(ji,jj,ik+1) = gdept(ji,jj,ik) + ze3tp ! ELSE ! standard case IF( bathy(ji,jj) <= gdepw_0(ik+1) ) THEN ; gdepw(ji,jj,ik+1) = bathy(ji,jj) ELSE ; gdepw(ji,jj,ik+1) = gdepw_0(ik+1) ENDIF !gm Bug? check the gdepw_0 ! ... on ik gdept(ji,jj,ik) = gdepw_0(ik) + ( gdepw (ji,jj,ik+1) - gdepw_0(ik) ) & & * ((gdept_0( ik ) - gdepw_0(ik) ) & & / ( gdepw_0( ik+1) - gdepw_0(ik) )) e3t (ji,jj,ik) = e3t_0 (ik) * ( gdepw (ji,jj,ik+1) - gdepw_0(ik) ) & & / ( gdepw_0( ik+1) - gdepw_0(ik) ) e3w (ji,jj,ik) = 0.5_wp * ( gdepw(ji,jj,ik+1) + gdepw_0(ik+1) - 2._wp * gdepw_0(ik) ) & & * ( e3w_0(ik) / ( gdepw_0(ik+1) - gdepw_0(ik) ) ) ! ... on ik+1 e3w (ji,jj,ik+1) = e3t (ji,jj,ik) e3t (ji,jj,ik+1) = e3t (ji,jj,ik) gdept(ji,jj,ik+1) = gdept(ji,jj,ik) + e3t(ji,jj,ik) ENDIF ENDIF END DO END DO ! jl = 0 DO jj = 1, jpj DO ji = 1, jpi ik = mbathy(ji,jj) IF( ik > jpksigm ) THEN ! ocean point only e3tp (ji,jj) = e3t(ji,jj,ik ) e3wp (ji,jj) = e3w(ji,jj,ik ) ! test zmin= gdepw(ji,jj,ik+1) - gdept(ji,jj,ik ) IF( zmin <= 0._wp .AND. lwp ) THEN jl = jl + 1 WRITE(numout,*) ' it = ', jl, ' ik = ', ik, ' (i,j) = ', ji, jj WRITE(numout,*) ' bathy = ', bathy(ji,jj) WRITE(numout,*) ' gdept = ', gdept(ji,jj,ik), ' gdepw = ', gdepw(ji,jj,ik+1), ' zdiff = ', zmin WRITE(numout,*) ' e3tp = ', e3t (ji,jj,ik), ' e3wp = ', e3w (ji,jj,ik ) ENDIF ENDIF END DO END DO ! Scale factors and depth at U-, V-, UW and VW-points DO jk = 1, jpk ! initialisation to z-scale factors e3u (:,:,jk) = e3t(:,:,jk) e3v (:,:,jk) = e3t(:,:,jk) e3uw(:,:,jk) = e3w(:,:,jk) e3vw(:,:,jk) = e3w(:,:,jk) e3f (:,:,jk) = e3t(:,:,jk) END DO DO jk = 1, jpksigm-1 DO jj = 1, jpjm1 DO ji = 1, fs_jpim1 e3u(ji,jj,jk)=( REAL(MIN(1,mbathy(ji,jj)),wp)* e3t(ji,jj,jk) + & REAL(MIN(1,mbathy(ji+1,jj)),wp)*e3t(ji+1,jj,jk) ) & /MAX( 1, MIN(1,mbathy(ji,jj))+MIN(1,mbathy(ji+1,jj)) ) e3uw(ji,jj,jk)=(REAL(MIN(1,mbathy(ji,jj)),wp)* e3w(ji,jj,jk) + & REAL(MIN(1,mbathy(ji+1,jj)),wp)*e3w(ji+1,jj,jk) ) & /REAL(MAX( 1, MIN(1,mbathy(ji,jj))+MIN(1,mbathy(ji+1,jj))),wp) e3v(ji,jj,jk)=(REAL(MIN(1,mbathy(ji,jj)),wp)* e3t(ji,jj,jk) + & REAL(MIN(1,mbathy(ji,jj+1)),wp)*e3t(ji,jj+1,jk) ) & /REAL (MAX( 1, MIN(1,mbathy(ji,jj))+MIN(1,mbathy(ji,jj+1))),wp) e3vw(ji,jj,jk)=(REAL(MIN(1,mbathy(ji,jj)),wp)* e3w(ji,jj,jk) + & REAL(MIN(1,mbathy(ji,jj+1)),wp)*e3w(ji,jj+1,jk) ) & /REAL(MAX( 1, MIN(1,mbathy(ji,jj))+MIN(1,mbathy(ji,jj+1))),wp) e3f(ji,jj,jk)=(REAL(MIN(1,mbathy(ji,jj)),wp)* e3t(ji,jj,jk) + & REAL(MIN(1,mbathy(ji+1,jj)),wp)*e3t(ji+1,jj,jk) + & REAL(MIN(1,mbathy(ji+1,jj+1)),wp)* e3t(ji+1,jj+1,jk)+ & REAL(MIN(1,mbathy(ji,jj+1)),wp)*e3t(ji,jj+1,jk) ) & /REAL(MAX( 1, MIN(1,mbathy(ji,jj))+MIN(1,mbathy(ji,jj+1)) & + MIN(1,mbathy(ji+1,jj))+MIN(1,mbathy(ji+1,jj+1))),wp) ENDDO END DO END DO DO jk = jpksigm,jpk ! Computed as the minimum of neighbooring scale factors DO jj = 1, jpjm1 DO ji = 1, fs_jpim1 ! vector opt. e3u (ji,jj,jk) = MIN( e3t(ji,jj,jk), e3t(ji+1,jj,jk) ) e3v (ji,jj,jk) = MIN( e3t(ji,jj,jk), e3t(ji,jj+1,jk) ) e3uw(ji,jj,jk) = MIN( e3w(ji,jj,jk), e3w(ji+1,jj,jk) ) e3vw(ji,jj,jk) = MIN( e3w(ji,jj,jk), e3w(ji,jj+1,jk) ) END DO END DO END DO CALL lbc_lnk( e3u , 'U', 1._wp ) ; CALL lbc_lnk( e3uw, 'U', 1._wp ) ! lateral boundary conditions CALL lbc_lnk( e3v , 'V', 1._wp ) ; CALL lbc_lnk( e3vw, 'V', 1._wp ) ! DO jk = 1, jpk ! set to z-scale factor if zero (i.e. along closed boundaries) WHERE( e3u (:,:,jk) == 0._wp ) e3u (:,:,jk) = e3t_0(jk) WHERE( e3v (:,:,jk) == 0._wp ) e3v (:,:,jk) = e3t_0(jk) WHERE( e3uw(:,:,jk) == 0._wp ) e3uw(:,:,jk) = e3w_0(jk) WHERE( e3vw(:,:,jk) == 0._wp ) e3vw(:,:,jk) = e3w_0(jk) END DO DO jk = jpksigm, jpk ! Computed as the minimum of neighbooring V-scale factors DO jj = 1, jpjm1 DO ji = 1, fs_jpim1 ! vector opt. e3f(ji,jj,jk) = MIN( e3v(ji,jj,jk), e3v(ji+1,jj,jk) ) END DO END DO END DO CALL lbc_lnk( e3f, 'F', 1._wp ) ! Lateral boundary conditions ! DO jk = 1, jpk ! set to z-scale factor if zero (i.e. along closed boundaries) WHERE( e3f(:,:,jk) == 0._wp ) e3f(:,:,jk) = e3t_0(jk) END DO !!gm bug ? : must be a do loop with mj0,mj1 ! e3t(:,mj0(1),:) = e3t(:,mj0(2),:) ! we duplicate factor scales for jj = 1 and jj = 2 e3w(:,mj0(1),:) = e3w(:,mj0(2),:) e3u(:,mj0(1),:) = e3u(:,mj0(2),:) e3v(:,mj0(1),:) = e3v(:,mj0(2),:) e3f(:,mj0(1),:) = e3f(:,mj0(2),:) ! Control of the sign Do jk=1,jpk do jj=1,jpj do ji=1,jpi IF( ( e3t (ji,jj,jk) ) <= 0._wp )then write(numout,*)' zgr_hyb : e r r o r e3t <= 0',ji,jj,jk,e3t (ji,jj,jk); endif IF( ( e3w (ji,jj,jk) ) <= 0._wp )then write(numout,*)' zgr_hyb : e r r o r e3t <= 0',ji,jj,jk,e3w (ji,jj,jk); endif IF( ( gdept(ji,jj,jk) ) < 0._wp )THEN write (numout,*)' zgr_hyb : e r r o r gdept < 0',ji,jj,jk ,gdept(ji,jj,jj); endif IF( ( gdepw(ji,jj,jk) ) < 0._wp )then write (numout,*)' zgr_hyb : e r r o r gdepw < 0',ji,jj,jk , gdepw(ji,jj,jj); endif enddo enddo enddo ! Compute gdep3w (vertical sum of e3w) gdep3w(:,:,1) = 0.5_wp * e3w(:,:,1) DO jk = 2, jpk gdep3w(:,:,jk) = gdep3w(:,:,jk-1) + e3w(:,:,jk) END DO IF( lk_mpp ) CALL mpp_max( nstop ) IF (lwp) write(numout,*)"zpartial" ,nstop CALL lbc_lnk( e3f, 'F', 1._wp ) ! Lateral boundary conditions CALL wrk_dealloc( jpi, jpj, zenv, ztmp, zmsk, zri, zrj, zhbat ) ! IF( nn_timing == 1 ) CALL timing_stop('zgr_hyb') END SUBROUTINE zgr_hyb !!====================================================================== END MODULE domzgr