PROGRAM SCOORD_GEN !!---------------------------------------------------------------------- !! *** PROGRAM scoord_gen *** !! !! ** 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 z_gsigt, z_gsigw, z_esigt, z_esigw from an analytical !! function and its derivative given as function. !! z_gsigt(k) = fssig (k ) !! z_gsigw(k) = fssig (k-0.5) !! z_esigt(k) = fsdsig(k ) !! z_esigw(k) = fsdsig(k-0.5) !! Three options for stretching are give, and they can be modified !! following the users requirements. 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. !! !! The three methods for stretching available are: !! !! s_sh94 (Song and Haidvogel 1994) !! a sinh/tanh function that allows sigma and stretched sigma !! !! s_sf12 (Siddorn and Furner 2012?) !! allows the maintenance of fixed surface and or !! bottom cell resolutions (cf. geopotential coordinates) !! within an analytically derived stretched S-coordinate framework. !! !! s_tanh (Madec et al 1996) !! a cosh/tanh function that gives stretched coordinates !! !!---------------------------------------------------------------------- ! USE utils IMPLICIT NONE !!---------------------------------------------------------------------- ! ! OPEN( UNIT=numnam, FILE='namelist', FORM='FORMATTED', STATUS='OLD' ) READ( numnam, namzgr_sco ) CLOSE( numnam ) jpk = INT(rn_jpk) WRITE(*,*) WRITE(*,*) 'domzgr_sco : s-coordinate or hybrid z-s-coordinate' WRITE(*,*) '~~~~~~~~~~~' WRITE(*,*) ' Namelist namzgr_sco' WRITE(*,*) ' stretching coeffs ' WRITE(*,*) ' Number of levels rn_jpk = ',rn_jpk WRITE(*,*) ' maximum depth of s-bottom surface (>0) rn_sbot_max = ',rn_sbot_max WRITE(*,*) ' minimum depth of s-bottom surface (>0) rn_sbot_min = ',rn_sbot_min WRITE(*,*) ' Critical depth rn_hc = ',rn_hc WRITE(*,*) ' maximum cut-off r-value allowed rn_rmax = ',rn_rmax WRITE(*,*) ' Song and Haidvogel 1994 stretching ln_s_sh94 = ',ln_s_sh94 WRITE(*,*) ' Song and Haidvogel 1994 stretching coefficients' WRITE(*,*) ' surface control parameter (0<=rn_theta<=20) rn_theta = ',rn_theta WRITE(*,*) ' bottom control parameter (0<=rn_thetb<= 1) rn_thetb = ',rn_thetb WRITE(*,*) ' stretching parameter (song and haidvogel) rn_bb = ',rn_bb WRITE(*,*) ' Siddorn and Furner 2012 stretching ln_s_sf12 = ',ln_s_sf12 WRITE(*,*) ' switching to sigma (T) or Z (F) at H1 surface; <1 bottom) rn_alpha = ',rn_alpha WRITE(*,*) ' e-fold length scale for transition region rn_efold = ',rn_efold WRITE(*,*) ' Surface cell depth (Zs) (m) rn_zs = ',rn_zs WRITE(*,*) ' Bathymetry multiplier for Zb rn_zb_a = ',rn_zb_a WRITE(*,*) ' Offset for Zb rn_zb_b = ',rn_zb_b WRITE(*,*) ' Bottom cell (Zb) (m) = H*rn_zb_a + rn_zb_b' ! Read in bathy, jpi and jpj from bathy.nc CALL read_bathy() !Allocate all other allocatable variables ios = dom_oce_alloc() IF( ios .NE. 0) THEN WRITE(0,*) 'Unable to allocate all arrays' STOP 1 ENDIF 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. ) 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(:,:) ! ! 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. ) 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. ) THEN zenv(ji,jj) = rn_sbot_min ENDIF ENDIF END DO END DO ! ! smooth the bathymetry (if required) scosrf(:,:) = 0. ! ocean surface depth (here zero: no under ice-shelf sea) scobot(:,:) = bathy(:,:) ! ocean bottom depth ! jl = 0 zrmax = 1. ! ! ! set scaling factor used in reducing vertical gradients zrfact = ( 1. - rn_rmax ) / ( 1. + rn_rmax ) ! ! initialise temporary evelope depth arrays ztmpi1(:,:) = zenv(:,:) ztmpi2(:,:) = zenv(:,:) ztmpj1(:,:) = zenv(:,:) ztmpj2(:,:) = zenv(:,:) ! ! initialise temporary r-value arrays zri(:,:) = 1. zrj(:,:) = 1. ! ! ================ ! DO WHILE( jl <= 10000 .AND. ( zrmax - rn_rmax ) > 1.e-8 ) ! Iterative loop ! ! ! ================ ! jl = jl + 1 zrmax = 0. ! 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 DO jj = 1, jpi !jpi or jpim1? DO ji = 1, jpj zrmax = MAX( zrmax, ABS(zri(ji,jj)), ABS(zrj(ji,jj)) ) END DO END DO zri(:,:) = 0. zrj(:,:) = 0. ! DO jj = 1, nlci ! DO ji = 1, nlcj ! 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.) .AND. (zenv(iip1,jj) > 0.)) THEN ! zri(ji,jj) = ( zenv(iip1,jj ) - zenv(ji,jj) ) / ( zenv(iip1,jj ) + zenv(ji,jj) ) ! END IF ! IF( (zenv(ji,jj) > 0.) .AND. (zenv(ji,ijp1) > 0.)) 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 DO jj = 1, jpi-1 DO ji = 1, jpj-1 iip1 = ji+1 ijp1 = jj+1 IF( (zenv(ji,jj) > 0.) .AND. (zenv(iip1,jj) > 0.)) THEN zri(ji,jj) = ( zenv(iip1,jj ) - zenv(ji,jj) ) / ( zenv(iip1,jj ) + zenv(ji,jj) ) END IF IF( (zenv(ji,jj) > 0.) .AND. (zenv(ji,ijp1) > 0.)) 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 ! WRITE(*,*) 'zgr_sco : iter= ',jl, ' rmax= ', zrmax ! DO jj = 1, jpi DO ji = 1, jpj zenv(ji,jj) = MAX(zenv(ji,jj), ztmpi1(ji,jj), ztmpi2(ji,jj), ztmpj1(ji,jj), ztmpj2(ji,jj) ) END DO END DO ! ! ================ ! 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 ! TODO - get this section to work hbatt(:,:) = zenv(:,:) ! IF( MINVAL( gphit(:,:) ) * MAXVAL( gphit(:,:) ) <= 0. ) 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.)**2. ) ! hbatt(ji,jj) = rn_sbot_max * ztaper + hbatt(ji,jj) * ( 1. - ztaper ) ! END DO ! END DO ! ENDIF ! ! ! ============================== ! ! hbatu, hbatv, hbatf fields ! ! ============================== WRITE(*,*) WRITE(*,*) ' zgr_sco: minimum depth of the envelop topography set to : ', rn_sbot_min hbatu(:,:) = rn_sbot_min hbatv(:,:) = rn_sbot_min hbatf(:,:) = rn_sbot_min DO jj = 1, jpj-1 DO ji = 1, jpi-1 ! NO vector opt. hbatu(ji,jj) = 0.50 * ( hbatt(ji ,jj) + hbatt(ji+1,jj ) ) hbatv(ji,jj) = 0.50 * ( hbatt(ji ,jj) + hbatt(ji ,jj+1) ) hbatf(ji,jj) = 0.25 * ( hbatt(ji ,jj) + hbatt(ji ,jj+1) & & + hbatt(ji+1,jj) + hbatt(ji+1,jj+1) ) END DO END DO ! zhbat(:,:) = hbatu(:,:) DO jj = 1, jpj DO ji = 1, jpi IF( hbatu(ji,jj) == 0. ) THEN IF( zhbat(ji,jj) == 0. ) hbatu(ji,jj) = rn_sbot_min IF( zhbat(ji,jj) /= 0. ) hbatu(ji,jj) = zhbat(ji,jj) ENDIF END DO END DO zhbat(:,:) = hbatv(:,:) DO jj = 1, jpj DO ji = 1, jpi IF( hbatv(ji,jj) == 0. ) THEN IF( zhbat(ji,jj) == 0. ) hbatv(ji,jj) = rn_sbot_min IF( zhbat(ji,jj) /= 0. ) hbatv(ji,jj) = zhbat(ji,jj) ENDIF END DO END DO zhbat(:,:) = hbatf(:,:) DO jj = 1, jpj DO ji = 1, jpi IF( hbatf(ji,jj) == 0. ) THEN IF( zhbat(ji,jj) == 0. ) hbatf(ji,jj) = rn_sbot_min IF( zhbat(ji,jj) /= 0. ) 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(:,:) ) WRITE(*,*) ' MAX val hif t ', MAXVAL( hift (:,:) ), ' f ', MAXVAL( hiff (:,:) ), & & ' u ', MAXVAL( hifu (:,:) ), ' v ', MAXVAL( hifv (:,:) ) WRITE(*,*) ' MIN val hif t ', MINVAL( hift (:,:) ), ' f ', MINVAL( hiff (:,:) ), & & ' u ', MINVAL( hifu (:,:) ), ' v ', MINVAL( hifv (:,:) ) WRITE(*,*) ' MAX val hbat t ', MAXVAL( hbatt(:,:) ), ' f ', MAXVAL( hbatf(:,:) ), & & ' u ', MAXVAL( hbatu(:,:) ), ' v ', MAXVAL( hbatv(:,:) ) WRITE(*,*) ' MIN val hbat t ', MINVAL( hbatt(:,:) ), ' f ', MINVAL( hbatf(:,:) ), & & ' u ', MINVAL( hbatu(:,:) ), ' v ', MINVAL( hbatv(:,:) ) !! helsinki ! ! ======================= ! ! s-coordinate fields (gdep., e3.) ! ! ======================= ! ! non-dimensional "sigma" for model level depth at w- and t-levels !======================================================================== ! Song and Haidvogel 1994 (ln_s_sh94=T) ! Siddorn and Furner 2012 (ln_sf12=T) ! or tanh function (both false) !======================================================================== IF ( ln_s_sh94 ) THEN CALL s_sh94() ELSE IF ( ln_s_sf12 ) THEN CALL s_sf12() ELSE CALL s_tanh() ENDIF ! where (e3t_0 (:,:,:).eq.0.0) e3t_0(:,:,:) = 1.0 where (e3u_0 (:,:,:).eq.0.0) e3u_0(:,:,:) = 1.0 where (e3v_0 (:,:,:).eq.0.0) e3v_0(:,:,:) = 1.0 where (e3f_0 (:,:,:).eq.0.0) e3f_0(:,:,:) = 1.0 where (e3w_0 (:,:,:).eq.0.0) e3w_0(:,:,:) = 1.0 where (e3uw_0 (:,:,:).eq.0.0) e3uw_0(:,:,:) = 1.0 where (e3vw_0 (:,:,:).eq.0.0) e3vw_0(:,:,:) = 1.0 #if defined key_agrif ! Ensure meaningful vertical scale factors in ghost lines/columns ! TODO: How can we deal with this in offline tool? IF( .NOT. Agrif_Root() ) THEN ! IF((nbondi == -1).OR.(nbondi == 2)) THEN e3u_0(1,:,:) = e3u_0(2,:,:) ENDIF ! IF((nbondi == 1).OR.(nbondi == 2)) THEN e3u_0(nlci-1,:,:) = e3u_0(nlci-2,:,:) ENDIF ! IF((nbondj == -1).OR.(nbondj == 2)) THEN e3v_0(:,1,:) = e3v_0(:,2,:) ENDIF ! IF((nbondj == 1).OR.(nbondj == 2)) THEN e3v_0(:,nlcj-1,:) = e3v_0(:,nlcj-2,:) ENDIF ! ENDIF #endif !! ! HYBRID : DO jj = 1, jpj DO ji = 1, jpi DO jk = 1, jpk-1 IF( scobot(ji,jj) >= gdept_0(ji,jj,jk) ) mbathy(ji,jj) = MAX( 2, jk ) END DO IF( scobot(ji,jj) == 0. ) mbathy(ji,jj) = 0 END DO END DO WRITE(*,*) ' MIN val mbathy h90 ', MINVAL( mbathy(:,:) ), ' MAX ', MAXVAL( mbathy(:,:) ) ! min max values over the domain WRITE(*,*) ' MIN val mbathy ', MINVAL( mbathy(:,:) ), ' MAX ', MAXVAL( mbathy(:,:) ) WRITE(*,*) ' MIN val depth t ', MINVAL( gdept_0(:,:,:) ), & & ' w ', MINVAL( gdepw_0(:,:,:) ), '3w ' , MINVAL( gdep3w_0(:,:,:) ) WRITE(*,*) ' MIN val e3 t ', MINVAL( e3t_0 (:,:,:) ), ' f ' , MINVAL( e3f_0 (:,:,:) ), & & ' u ', MINVAL( e3u_0 (:,:,:) ), ' u ' , MINVAL( e3v_0 (:,:,:) ), & & ' uw', MINVAL( e3uw_0 (:,:,:) ), ' vw' , MINVAL( e3vw_0 (:,:,:) ), & & ' w ', MINVAL( e3w_0 (:,:,:) ) WRITE(*,*) ' MAX val depth t ', MAXVAL( gdept_0(:,:,:) ), & & ' w ', MAXVAL( gdepw_0(:,:,:) ), '3w ' , MAXVAL( gdep3w_0(:,:,:) ) WRITE(*,*) ' MAX val e3 t ', MAXVAL( e3t_0 (:,:,:) ), ' f ' , MAXVAL( e3f_0 (:,:,:) ), & & ' u ', MAXVAL( e3u_0 (:,:,:) ), ' u ' , MAXVAL( e3v_0 (:,:,:) ), & & ' uw', MAXVAL( e3uw_0 (:,:,:) ), ' vw' , MAXVAL( e3vw_0 (:,:,:) ), & & ' w ', MAXVAL( e3w_0 (:,:,:) ) ! selected vertical profiles WRITE(*,*) WRITE(*,*) ' domzgr: vertical coordinates : point (1,1,k) bathy = ', bathy(1,1), hbatt(1,1) WRITE(*,*) ' ~~~~~~ --------------------' WRITE(*,"(9x,' level gdept_0 gdepw_0 e3t_0 e3w_0')") WRITE(*,"(10x,i4,4f9.2)") ( jk, gdept_0(1,1,jk), gdepw_0(1,1,jk), & & e3t_0 (1,1,jk) , e3w_0 (1,1,jk) , jk=1,jpk ) DO jj = 20, 20 DO ji = 20, 20 WRITE(*,*) WRITE(*,*) ' domzgr: vertical coordinates : point (20,20,k) bathy = ', bathy(ji,jj), hbatt(ji,jj) WRITE(*,*) ' ~~~~~~ --------------------' WRITE(*,"(9x,' level gdept_0 gdepw_0 e3t_0 e3w_0')") WRITE(*,"(10x,i4,4f9.2)") ( jk, gdept_0(ji,jj,jk), gdepw_0(ji,jj,jk), & & e3t_0 (ji,jj,jk) , e3w_0 (ji,jj,jk) , jk=1,jpk ) END DO END DO DO jj = 74,74 DO ji = 100, 100 WRITE(*,*) WRITE(*,*) ' domzgr: vertical coordinates : point (100,74,k) bathy = ', bathy(ji,jj), hbatt(ji,jj) WRITE(*,*) ' ~~~~~~ --------------------' WRITE(*,"(9x,' level gdept_0 gdepw_0 e3t_0 e3w_0')") WRITE(*,"(10x,i4,4f9.2)") ( jk, gdept_0(ji,jj,jk), gdepw_0(ji,jj,jk), & & e3t_0 (ji,jj,jk) , e3w_0 (ji,jj,jk) , jk=1,jpk ) END DO END DO !================================================================================ ! check the coordinate makes sense !================================================================================ DO ji = 1, jpi DO jj = 1, jpj IF( hbatt(ji,jj) > 0.) THEN DO jk = 1, mbathy(ji,jj) ! check coordinate is monotonically increasing IF (e3w_0(ji,jj,jk) <= 0. .OR. e3t_0(ji,jj,jk) <= 0. ) THEN WRITE(*,*) 'ERROR zgr_sco : e3w or e3t =< 0 at point (i,j,k)= ', ji, jj, jk WRITE(*,*) 'e3w',e3w_0(ji,jj,:) WRITE(*,*) 'e3t',e3t_0(ji,jj,:) STOP 1 ENDIF ! and check it has never gone negative IF( gdepw_0(ji,jj,jk) < 0. .OR. gdept_0(ji,jj,jk) < 0. ) THEN WRITE(*,*) 'ERROR zgr_sco : gdepw or gdept =< 0 at point (i,j,k)= ', ji, jj, jk WRITE(*,*) 'gdepw',gdepw_0(ji,jj,:) WRITE(*,*) 'gdept',gdept_0(ji,jj,:) STOP 1 ENDIF ! and check it never exceeds the total depth IF( gdepw_0(ji,jj,jk) > hbatt(ji,jj) ) THEN WRITE(*,*) 'ERROR zgr_sco : gdepw > hbatt at point (i,j,k)= ', ji, jj, jk WRITE(*,*) 'gdepw',gdepw_0(ji,jj,:) STOP 1 ENDIF END DO DO jk = 1, mbathy(ji,jj)-1 ! and check it never exceeds the total depth IF( gdept_0(ji,jj,jk) > hbatt(ji,jj) ) THEN WRITE(*,*) 'ERROR zgr_sco : gdept > hbatt at point (i,j,k)= ', ji, jj, jk WRITE(*,*) 'gdept',gdept_0(ji,jj,:) STOP 1 ENDIF END DO ENDIF END DO END DO ! ! Write output file CALL write_coord_file() ! ! END PROGRAM SCOORD_GEN !!====================================================================== SUBROUTINE s_sh94() !!---------------------------------------------------------------------- !! *** ROUTINE s_sh94 *** !! !! ** Purpose : stretch the s-coordinate system !! !! ** Method : s-coordinate stretch using the Song and Haidvogel 1994 !! mixed S/sigma coordinate !! !! Reference : Song and Haidvogel 1994. !!---------------------------------------------------------------------- ! USE utils REAL(wp) :: zcoeft, zcoefw ! temporary scalars ! REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: z_gsigw3, z_gsigt3, z_gsi3w3 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: z_esigt3, z_esigw3, z_esigtu3 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 ALLOCATE( z_gsigw3(jpi,jpj,jpk), z_gsigt3(jpi,jpj,jpk), z_gsi3w3(jpi,jpj,jpk) ) ALLOCATE( z_esigt3(jpi,jpj,jpk), z_esigw3(jpi,jpj,jpk), z_esigtu3(jpi,jpj,jpk)) ALLOCATE( z_esigtv3(jpi,jpj,jpk), z_esigtf3(jpi,jpj,jpk), z_esigwu3(jpi,jpj,jpk)) ALLOCATE( z_esigwv3(jpi,jpj,jpk) ) z_gsigw3 = 0. ; z_gsigt3 = 0. ; z_gsi3w3 = 0. z_esigt3 = 0. ; z_esigw3 = 0. z_esigtu3 = 0. ; z_esigtv3 = 0. ; z_esigtf3 = 0. z_esigwu3 = 0. ; z_esigwv3 = 0. DO ji = 1, jpi DO jj = 1, jpj IF( hbatt(ji,jj) > rn_hc ) THEN !deep water, stretched sigma DO jk = 1, jpk z_gsigw3(ji,jj,jk) = -fssig1( REAL(jk,wp)-0.5, rn_bb ) z_gsigt3(ji,jj,jk) = -fssig1( REAL(jk,wp) , rn_bb ) END DO ELSE ! shallow water, uniform sigma DO jk = 1, jpk z_gsigw3(ji,jj,jk) = REAL(jk-1,wp) / REAL(jpk-1,wp) z_gsigt3(ji,jj,jk) = ( REAL(jk-1,wp) + 0.5 ) / REAL(jpk-1,wp) END DO ENDIF ! DO jk = 1, jpk-1 z_esigt3(ji,jj,jk ) = z_gsigw3(ji,jj,jk+1) - z_gsigw3(ji,jj,jk) z_esigw3(ji,jj,jk+1) = z_gsigt3(ji,jj,jk+1) - z_gsigt3(ji,jj,jk) END DO z_esigw3(ji,jj,1 ) = 2. * ( z_gsigt3(ji,jj,1 ) - z_gsigw3(ji,jj,1 ) ) z_esigt3(ji,jj,jpk) = 2. * ( z_gsigt3(ji,jj,jpk) - z_gsigw3(ji,jj,jpk) ) ! ! Coefficients for vertical depth as the sum of e3w scale factors z_gsi3w3(ji,jj,1) = 0.5 * z_esigw3(ji,jj,1) DO jk = 2, jpk z_gsi3w3(ji,jj,jk) = z_gsi3w3(ji,jj,jk-1) + z_esigw3(ji,jj,jk) END DO ! DO jk = 1, jpk zcoeft = ( REAL(jk,wp) - 0.5 ) / REAL(jpk-1,wp) zcoefw = ( REAL(jk,wp) - 1.0 ) / REAL(jpk-1,wp) gdept_0(ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*z_gsigt3(ji,jj,jk)+rn_hc*zcoeft ) gdepw_0(ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*z_gsigw3(ji,jj,jk)+rn_hc*zcoefw ) gdep3w_0(ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*z_gsi3w3(ji,jj,jk)+rn_hc*zcoeft ) END DO ! END DO ! for all jj's END DO ! for all ji's DO ji = 1, jpim1 DO jj = 1, jpjm1 DO jk = 1, jpk z_esigtu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk) ) & & / ( hbatt(ji,jj)+hbatt(ji+1,jj) ) z_esigtv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk) ) & & / ( hbatt(ji,jj)+hbatt(ji,jj+1) ) z_esigtf3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk) & & + hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk)+hbatt(ji+1,jj+1)*z_esigt3(ji+1,jj+1,jk) ) & & / ( hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1) ) z_esigwu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigw3(ji+1,jj,jk) ) & & / ( hbatt(ji,jj)+hbatt(ji+1,jj) ) z_esigwv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigw3(ji,jj+1,jk) ) & & / ( hbatt(ji,jj)+hbatt(ji,jj+1) ) ! e3t_0(ji,jj,jk) = ( (hbatt(ji,jj)-rn_hc)*z_esigt3 (ji,jj,jk) + rn_hc/REAL(jpk-1,wp) ) e3u_0(ji,jj,jk) = ( (hbatu(ji,jj)-rn_hc)*z_esigtu3(ji,jj,jk) + rn_hc/REAL(jpk-1,wp) ) e3v_0(ji,jj,jk) = ( (hbatv(ji,jj)-rn_hc)*z_esigtv3(ji,jj,jk) + rn_hc/REAL(jpk-1,wp) ) e3f_0(ji,jj,jk) = ( (hbatf(ji,jj)-rn_hc)*z_esigtf3(ji,jj,jk) + rn_hc/REAL(jpk-1,wp) ) ! e3w_0 (ji,jj,jk) = ( (hbatt(ji,jj)-rn_hc)*z_esigw3 (ji,jj,jk) + rn_hc/REAL(jpk-1,wp) ) e3uw_0(ji,jj,jk) = ( (hbatu(ji,jj)-rn_hc)*z_esigwu3(ji,jj,jk) + rn_hc/REAL(jpk-1,wp) ) e3vw_0(ji,jj,jk) = ( (hbatv(ji,jj)-rn_hc)*z_esigwv3(ji,jj,jk) + rn_hc/REAL(jpk-1,wp) ) END DO END DO END DO DEALLOCATE( z_gsigw3, z_gsigt3, z_gsi3w3) DEALLOCATE( z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 ) END SUBROUTINE s_sh94 SUBROUTINE s_sf12 !!---------------------------------------------------------------------- !! *** ROUTINE s_sf12 *** !! !! ** Purpose : stretch the s-coordinate system !! !! ** Method : s-coordinate stretch using the Siddorn and Furner 2012? !! mixed S/sigma/Z coordinate !! !! This method allows the maintenance of fixed surface and or !! bottom cell resolutions (cf. geopotential coordinates) !! within an analytically derived stretched S-coordinate framework. !! !! !! Reference : Siddorn and Furner 2012 (submitted Ocean modelling). !!---------------------------------------------------------------------- ! USE utils REAL(wp) :: zsmth ! smoothing around critical depth REAL(wp) :: zzs, zzb ! Surface and bottom cell thickness in sigma space ! REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: z_gsigw3, z_gsigt3, z_gsi3w3 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: z_esigt3, z_esigw3, z_esigtu3 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 ! ALLOCATE( z_gsigw3(jpi,jpj,jpk), z_gsigt3(jpi,jpj,jpk), z_gsi3w3(jpi,jpj,jpk) ) ALLOCATE( z_esigt3(jpi,jpj,jpk), z_esigw3(jpi,jpj,jpk), z_esigtu3(jpi,jpj,jpk)) ALLOCATE( z_esigtv3(jpi,jpj,jpk), z_esigtf3(jpi,jpj,jpk), z_esigwu3(jpi,jpj,jpk)) ALLOCATE( z_esigwv3(jpi,jpj,jpk) ) z_gsigw3 = 0. ; z_gsigt3 = 0. ; z_gsi3w3 = 0. z_esigt3 = 0. ; z_esigw3 = 0. z_esigtu3 = 0. ; z_esigtv3 = 0. ; z_esigtf3 = 0. z_esigwu3 = 0. ; z_esigwv3 = 0. DO ji = 1, jpi DO jj = 1, jpj IF (hbatt(ji,jj)>rn_hc) THEN !deep water, stretched sigma zzb = hbatt(ji,jj)*rn_zb_a + rn_zb_b ! this forces a linear bottom cell depth relationship with H,. ! could be changed by users but care must be taken to do so carefully zzb = 1.0-(zzb/hbatt(ji,jj)) zzs = rn_zs / hbatt(ji,jj) IF (rn_efold /= 0.0) THEN zsmth = tanh( (hbatt(ji,jj)- rn_hc ) / rn_efold ) ELSE zsmth = 1.0 ENDIF DO jk = 1, jpk z_gsigw3(ji,jj,jk) = REAL(jk-1,wp) /REAL(jpk-1,wp) z_gsigt3(ji,jj,jk) = (REAL(jk-1,wp)+0.5)/REAL(jpk-1,wp) ENDDO z_gsigw3(ji,jj,:) = fgamma( z_gsigw3(ji,jj,:), zzb, zzs, zsmth ) z_gsigt3(ji,jj,:) = fgamma( z_gsigt3(ji,jj,:), zzb, zzs, zsmth ) ELSE IF(ln_sigcrit) THEN ! shallow water, uniform sigma DO jk = 1, jpk z_gsigw3(ji,jj,jk) = REAL(jk-1,wp) /REAL(jpk-1,wp) z_gsigt3(ji,jj,jk) = (REAL(jk-1,wp)+0.5)/REAL(jpk-1,wp) END DO ELSE ! shallow water, z coordinates DO jk = 1, jpk z_gsigw3(ji,jj,jk) = REAL(jk-1,wp) /REAL(jpk-1,wp)*(rn_hc/hbatt(ji,jj)) z_gsigt3(ji,jj,jk) = (REAL(jk-1,wp)+0.5)/REAL(jpk-1,wp)*(rn_hc/hbatt(ji,jj)) END DO ENDIF DO jk = 1, jpk-1 z_esigt3(ji,jj,jk) = z_gsigw3(ji,jj,jk+1) - z_gsigw3(ji,jj,jk) z_esigw3(ji,jj,jk+1) = z_gsigt3(ji,jj,jk+1) - z_gsigt3(ji,jj,jk) END DO z_esigw3(ji,jj,1 ) = 2.0 * (z_gsigt3(ji,jj,1 ) - z_gsigw3(ji,jj,1 )) z_esigt3(ji,jj,jpk) = 2.0 * (z_gsigt3(ji,jj,jpk) - z_gsigw3(ji,jj,jpk)) ! Coefficients for vertical depth as the sum of e3w scale factors z_gsi3w3(ji,jj,1) = 0.5 * z_esigw3(ji,jj,1) DO jk = 2, jpk z_gsi3w3(ji,jj,jk) = z_gsi3w3(ji,jj,jk-1) + z_esigw3(ji,jj,jk) END DO DO jk = 1, jpk gdept_0(ji,jj,jk) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsigt3(ji,jj,jk) gdepw_0(ji,jj,jk) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsigw3(ji,jj,jk) gdep3w_0(ji,jj,jk) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsi3w3(ji,jj,jk) END DO ENDDO ! for all jj's ENDDO ! for all ji's DO ji=1,jpi-1 DO jj=1,jpj-1 DO jk = 1, jpk z_esigtu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk) ) / & ( hbatt(ji,jj)+hbatt(ji+1,jj) ) z_esigtv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk) ) / & ( hbatt(ji,jj)+hbatt(ji,jj+1) ) z_esigtf3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk) + & hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk)+hbatt(ji+1,jj+1)*z_esigt3(ji+1,jj+1,jk) ) / & ( hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1) ) z_esigwu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigw3(ji+1,jj,jk) ) / & ( hbatt(ji,jj)+hbatt(ji+1,jj) ) z_esigwv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigw3(ji,jj+1,jk) ) / & ( hbatt(ji,jj)+hbatt(ji,jj+1) ) e3t_0(ji,jj,jk)=(scosrf(ji,jj)+hbatt(ji,jj))*z_esigt3(ji,jj,jk) e3u_0(ji,jj,jk)=(scosrf(ji,jj)+hbatu(ji,jj))*z_esigtu3(ji,jj,jk) e3v_0(ji,jj,jk)=(scosrf(ji,jj)+hbatv(ji,jj))*z_esigtv3(ji,jj,jk) e3f_0(ji,jj,jk)=(scosrf(ji,jj)+hbatf(ji,jj))*z_esigtf3(ji,jj,jk) ! e3w_0(ji,jj,jk)=hbatt(ji,jj)*z_esigw3(ji,jj,jk) e3uw_0(ji,jj,jk)=hbatu(ji,jj)*z_esigwu3(ji,jj,jk) e3vw_0(ji,jj,jk)=hbatv(ji,jj)*z_esigwv3(ji,jj,jk) END DO ENDDO ENDDO ! ! ! ============= DEALLOCATE( z_gsigw3, z_gsigt3, z_gsi3w3) DEALLOCATE( z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 ) END SUBROUTINE s_sf12 SUBROUTINE s_tanh() !!---------------------------------------------------------------------- !! *** ROUTINE s_tanh*** !! !! ** Purpose : stretch the s-coordinate system !! !! ** Method : s-coordinate stretch !! !! Reference : Madec, Lott, Delecluse and Crepon, 1996. JPO, 26, 1393-1408. !!---------------------------------------------------------------------- USE utils REAL(wp) :: zcoeft, zcoefw ! temporary scalars REAL(wp), ALLOCATABLE, DIMENSION(:) :: z_gsigw, z_gsigt, z_gsi3w REAL(wp), ALLOCATABLE, DIMENSION(:) :: z_esigt, z_esigw ALLOCATE( z_gsigw(jpk), z_gsigt(jpk), z_gsi3w(jpk) ) ALLOCATE( z_esigt(jpk), z_esigw(jpk) ) z_gsigw = 0. ; z_gsigt = 0. ; z_gsi3w = 0. z_esigt = 0. ; z_esigw = 0. DO jk = 1, jpk z_gsigw(jk) = -fssig( REAL(jk,wp)-0.5 ) z_gsigt(jk) = -fssig( REAL(jk,wp) ) END DO WRITE(*,*) 'z_gsigw 1 jpk ', z_gsigw(1), z_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, jpk-1 z_esigt(jk ) = z_gsigw(jk+1) - z_gsigw(jk) z_esigw(jk+1) = z_gsigt(jk+1) - z_gsigt(jk) END DO z_esigw( 1 ) = 2. * ( z_gsigt(1 ) - z_gsigw(1 ) ) z_esigt(jpk) = 2. * ( z_gsigt(jpk) - z_gsigw(jpk) ) ! ! Coefficients for vertical depth as the sum of e3w scale factors z_gsi3w(1) = 0.5 * z_esigw(1) DO jk = 2, jpk z_gsi3w(jk) = z_gsi3w(jk-1) + z_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 ) / REAL(jpk-1,wp) zcoefw = ( REAL(jk,wp) - 1.0 ) / REAL(jpk-1,wp) gdept_0(:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*z_gsigt(jk) + hift(:,:)*zcoeft ) gdepw_0(:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*z_gsigw(jk) + hift(:,:)*zcoefw ) gdep3w_0(:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*z_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_0(ji,jj,jk) = ( (hbatt(ji,jj)-hift(ji,jj))*z_esigt(jk) + hift(ji,jj)/REAL(jpk-1,wp) ) e3u_0(ji,jj,jk) = ( (hbatu(ji,jj)-hifu(ji,jj))*z_esigt(jk) + hifu(ji,jj)/REAL(jpk-1,wp) ) e3v_0(ji,jj,jk) = ( (hbatv(ji,jj)-hifv(ji,jj))*z_esigt(jk) + hifv(ji,jj)/REAL(jpk-1,wp) ) e3f_0(ji,jj,jk) = ( (hbatf(ji,jj)-hiff(ji,jj))*z_esigt(jk) + hiff(ji,jj)/REAL(jpk-1,wp) ) ! e3w_0(ji,jj,jk) = ( (hbatt(ji,jj)-hift(ji,jj))*z_esigw(jk) + hift(ji,jj)/REAL(jpk-1,wp) ) e3uw_0(ji,jj,jk) = ( (hbatu(ji,jj)-hifu(ji,jj))*z_esigw(jk) + hifu(ji,jj)/REAL(jpk-1,wp) ) e3vw_0(ji,jj,jk) = ( (hbatv(ji,jj)-hifv(ji,jj))*z_esigw(jk) + hifv(ji,jj)/REAL(jpk-1,wp) ) END DO END DO END DO DEALLOCATE( z_gsigw, z_gsigt, z_gsi3w ) DEALLOCATE( z_esigt, z_esigw ) END SUBROUTINE s_tanh FUNCTION fssig( pk ) RESULT( pf ) !!---------------------------------------------------------------------- !! *** ROUTINE fssig *** !! !! ** 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) !!---------------------------------------------------------------------- USE utils, ONLY : wp REAL(wp), INTENT(in) :: pk ! continuous "k" coordinate REAL(wp) :: pf ! sigma value !!---------------------------------------------------------------------- ! pf = ( TANH( rn_theta * ( -(pk-0.5) / REAL(jpk-1) + rn_thetb ) ) & & - TANH( rn_thetb * rn_theta ) ) & & * ( COSH( rn_theta ) & & + COSH( rn_theta * ( 2. * rn_thetb - 1. ) ) ) & & / ( 2. * SINH( rn_theta ) ) ! END FUNCTION fssig FUNCTION fssig1( pk1, pbb ) RESULT( pf1 ) !!---------------------------------------------------------------------- !! *** ROUTINE fssig1 *** !! !! ** 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) !!---------------------------------------------------------------------- USE utils, ONLY : wp REAL(wp), INTENT(in) :: pk1 ! continuous "k" coordinate REAL(wp), INTENT(in) :: pbb ! Stretching coefficient REAL(wp) :: pf1 ! sigma value !!---------------------------------------------------------------------- ! IF ( rn_theta == 0 ) then ! uniform sigma pf1 = - ( pk1 - 0.5 ) / REAL( jpk-1 ) ELSE ! stretched sigma pf1 = ( 1. - pbb ) * ( SINH( rn_theta*(-(pk1-0.5)/REAL(jpk-1)) ) ) / SINH( rn_theta ) & & + pbb * ( (TANH( rn_theta*( (-(pk1-0.5)/REAL(jpk-1)) + 0.5) ) - TANH( 0.5 * rn_theta ) ) & & / ( 2. * TANH( 0.5 * rn_theta ) ) ) ENDIF ! END FUNCTION fssig1 FUNCTION fgamma( pk1, pzb, pzs, psmth) RESULT( p_gamma ) !!---------------------------------------------------------------------- !! *** ROUTINE fgamma *** !! !! ** Purpose : provide analytical function for the 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) !! !! This method allows the maintenance of fixed surface and or !! bottom cell resolutions (cf. geopotential coordinates) !! within an analytically derived stretched S-coordinate framework. !! !! Reference : Siddorn and Furner, in prep !!---------------------------------------------------------------------- USE utils, ONLY : jpk,jk,wp REAL(wp), INTENT(in ) :: pk1(jpk) ! continuous "k" coordinate REAL(wp) :: p_gamma(jpk) ! stretched coordinate REAL(wp), INTENT(in ) :: pzb ! Bottom box depth REAL(wp), INTENT(in ) :: pzs ! surface box depth REAL(wp), INTENT(in ) :: psmth ! Smoothing parameter REAL(wp) :: za1,za2,za3 ! local variables REAL(wp) :: zn1,zn2 ! local variables REAL(wp) :: za,zb,zx ! local variables !!---------------------------------------------------------------------- ! zn1 = 1./(jpk-1.) zn2 = 1. - zn1 za1 = (rn_alpha+2.0)*zn1**(rn_alpha+1.0)-(rn_alpha+1.0)*zn1**(rn_alpha+2.0) za2 = (rn_alpha+2.0)*zn2**(rn_alpha+1.0)-(rn_alpha+1.0)*zn2**(rn_alpha+2.0) za3 = (zn2**3.0 - za2)/( zn1**3.0 - za1) za = pzb - za3*(pzs-za1)-za2 za = za/( zn2-0.5*(za2+zn2**2.0) - za3*(zn1-0.5*(za1+zn1**2.0) ) ) zb = (pzs - za1 - za*( zn1-0.5*(za1+zn1**2.0 ) ) ) / (zn1**3.0 - za1) zx = 1.0-za/2.0-zb DO jk = 1, jpk p_gamma(jk) = za*(pk1(jk)*(1.0-pk1(jk)/2.0))+zb*pk1(jk)**3.0 + & & zx*( (rn_alpha+2.0)*pk1(jk)**(rn_alpha+1.0)- & & (rn_alpha+1.0)*pk1(jk)**(rn_alpha+2.0) ) p_gamma(jk) = p_gamma(jk)*psmth+pk1(jk)*(1.0-psmth) ENDDO ! END FUNCTION fgamma