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(:,:) ) ! Create the output file CALL make_coord_file() ! ! ======================= ! ! 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) ! To reduce memory loop over jpk and write each level to file !======================================================================== IF ( ln_s_sh94 ) THEN CALL s_sh94() ELSE IF ( ln_s_sf12 ) THEN CALL s_sf12() ELSE CALL s_tanh() ENDIF CALL check_nf90( nf90_put_var(ncout, var_ids(11), mbathy) ) CALL check_nf90( nf90_close(ncout) ) ! 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_gsigw3p1, z_gsigt3m1, z_gsi3w3m1 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), z_gsigt3(jpi,jpj), z_gsi3w3(jpi,jpj) ) ALLOCATE( z_esigt3(jpi,jpj), z_esigw3(jpi,jpj), z_esigtu3(jpi,jpj)) ALLOCATE( z_esigtv3(jpi,jpj), z_esigtf3(jpi,jpj), z_esigwu3(jpi,jpj)) ALLOCATE( z_esigwv3(jpi,jpj), z_gsigw3p1(jpi,jpj), z_gsigt3m1(jpi,jpj) ) ALLOCATE( z_gsi3w3m1(jpi,jpj) ) z_gsigw3 = 0. ; z_gsigt3 = 0. ; z_gsi3w3 = 0. z_esigt3 = 0. ; z_esigw3 = 0. z_esigt3p1 = 0. ; z_esigw3p1 = 0. z_esigtu3 = 0. ; z_esigtv3 = 0. ; z_esigtf3 = 0. z_esigwu3 = 0. ; z_esigwv3 = 0. DO jk = 1,jpk DO ji = 1, jpi DO jj = 1, jpj IF( hbatt(ji,jj) > rn_hc ) THEN !deep water, stretched sigma z_gsigw3(ji,jj) = -fssig1( REAL(jk,wp)-0.5, rn_bb ) z_gsigw3p1(ji,jj) = -fssig1( REAL(jk+1,wp)-0.5, rn_bb ) z_gsigt3(ji,jj) = -fssig1( REAL(jk,wp) , rn_bb ) ELSE ! shallow water, uniform sigma z_gsigw3(ji,jj) = REAL(jk-1,wp) / REAL(jpk-1,wp) z_gsigw3p1(ji,jj) = REAL(jk,wp) / REAL(jpk-1,wp) z_gsigt3(ji,jj) = ( REAL(jk-1,wp) + 0.5 ) / REAL(jpk-1,wp) ENDIF ! !gsi3w3m1 & gsigt3m1 only used if jk /= 1 and is set at the end of the loop over jk IF( jk .EQ. 1) THEN z_esigw3(ji,jj ) = 2. * ( z_gsigt3(ji,jj ) - z_gsigw3(ji,jj ) ) z_gsi3w3(ji,jj) = 0.5 * z_esigw3(ji,jj) ELSE z_esigw3(ji,jj) = z_gsigt3(ji,jj) - z_gsigt3m1(ji,jj) z_gsi3w3(ji,jj) = z_gsi3w3m1(ji,jj) + z_esigw3(ji,jj) ENDIF IF( jk .EQ. jpk) THEN z_esigt3(ji,jj) = 2. * ( z_gsigt3(ji,jj) - z_gsigw3(ji,jj) ) ELSE z_esigt3(ji,jj ) = z_gsigw3p1(ji,jj) - z_gsigw3(ji,jj) ENDIF ! zcoeft = ( REAL(jk,wp) - 0.5 ) / REAL(jpk-1,wp) zcoefw = ( REAL(jk,wp) - 1.0 ) / REAL(jpk-1,wp) gdept_0(ji,jj) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*z_gsigt3(ji,jj)+rn_hc*zcoeft ) gdepw_0(ji,jj) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*z_gsigw3(ji,jj)+rn_hc*zcoefw ) gdep3w_0(ji,jj) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*z_gsi3w3(ji,jj)+rn_hc*zcoeft ) ! END DO ! for all jj's END DO ! for all ji's DO ji = 1, jpim1 DO jj = 1, jpjm1 z_esigtu3(ji,jj) = ( hbatt(ji,jj)*z_esigt3(ji,jj)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj) ) & & / ( hbatt(ji,jj)+hbatt(ji+1,jj) ) z_esigtv3(ji,jj) = ( hbatt(ji,jj)*z_esigt3(ji,jj)+hbatt(ji,jj+1)*z_esigt3(ji,jj+1) ) & & / ( hbatt(ji,jj)+hbatt(ji,jj+1) ) z_esigtf3(ji,jj) = ( hbatt(ji,jj)*z_esigt3(ji,jj)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj) & & + hbatt(ji,jj+1)*z_esigt3(ji,jj+1)+hbatt(ji+1,jj+1)*z_esigt3(ji+1,jj+1) ) & & / ( hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1) ) z_esigwu3(ji,jj) = ( hbatt(ji,jj)*z_esigw3(ji,jj)+hbatt(ji+1,jj)*z_esigw3(ji+1,jj) ) & & / ( hbatt(ji,jj)+hbatt(ji+1,jj) ) z_esigwv3(ji,jj) = ( hbatt(ji,jj)*z_esigw3(ji,jj)+hbatt(ji,jj+1)*z_esigw3(ji,jj+1) ) & & / ( hbatt(ji,jj)+hbatt(ji,jj+1) ) ! e3t_0(ji,jj) = ( (hbatt(ji,jj)-rn_hc)*z_esigt3 (ji,jj) + rn_hc/REAL(jpk-1,wp) ) e3u_0(ji,jj) = ( (hbatu(ji,jj)-rn_hc)*z_esigtu3(ji,jj) + rn_hc/REAL(jpk-1,wp) ) e3v_0(ji,jj) = ( (hbatv(ji,jj)-rn_hc)*z_esigtv3(ji,jj) + rn_hc/REAL(jpk-1,wp) ) e3f_0(ji,jj) = ( (hbatf(ji,jj)-rn_hc)*z_esigtf3(ji,jj) + rn_hc/REAL(jpk-1,wp) ) ! e3w_0 (ji,jj) = ( (hbatt(ji,jj)-rn_hc)*z_esigw3 (ji,jj) + rn_hc/REAL(jpk-1,wp) ) e3uw_0(ji,jj) = ( (hbatu(ji,jj)-rn_hc)*z_esigwu3(ji,jj) + rn_hc/REAL(jpk-1,wp) ) e3vw_0(ji,jj) = ( (hbatv(ji,jj)-rn_hc)*z_esigwv3(ji,jj) + rn_hc/REAL(jpk-1,wp) ) END DO END DO z_gsigt3m1 = z_gsigt3 z_gsiw3m1 = z_gsiw3 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 CALL write_netcdf_vars(jk) DO jj = 1, jpj DO ji = 1, jpi IF( scobot(ji,jj) >= gdept_0(ji,jj) ) mbathy(ji,jj) = MAX( 2, jk ) IF( scobot(ji,jj) == 0. ) mbathy(ji,jj) = 0 END DO END DO END DO !End of loop over jk DEALLOCATE( z_gsigw3, z_gsigt3, z_gsi3w3, z_gsigw3p1, z_gsigt3m1, z_gsi3w3m1) 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_gsigw3p1, z_gsigt3m1, z_gsi3w3m1 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), z_gsigt3(jpi,jpj), z_gsi3w3(jpi,jpj) ) ALLOCATE( z_esigt3(jpi,jpj), z_esigw3(jpi,jpj), z_esigtu3(jpi,jpj)) ALLOCATE( z_esigtv3(jpi,jpj), z_esigtf3(jpi,jpj), z_esigwu3(jpi,jpj)) ALLOCATE( z_esigwv3(jpi,jpj), z_gsigw3p1(jpi,jpj), z_gsigt3m1(jpi,jpj) ) ALLOCATE( z_gsi3w3m1(jpi,jpj) ) 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 jk = 1,jpk 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 z_gsigw3(ji,jj) = REAL(jk-1,wp) /REAL(jpk-1,wp) z_gsigw3p1(ji,jj) = REAL(jk,wp) /REAL(jpk-1,wp) z_gsigt3(ji,jj) = (REAL(jk-1,wp)+0.5)/REAL(jpk-1,wp) z_gsigw3(ji,jj) = fgamma( z_gsigw3(ji,jj), zzb, zzs, zsmth ) z_gsigw3p1(ji,jj) = fgamma( z_gsigw3p1(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 z_gsigw3(ji,jj) = REAL(jk-1,wp) /REAL(jpk-1,wp) z_gsigw3p1(ji,jj) = REAL(jk,wp) /REAL(jpk-1,wp) z_gsigt3(ji,jj) = (REAL(jk-1,wp)+0.5)/REAL(jpk-1,wp) ELSE ! shallow water, z coordinates z_gsigw3(ji,jj) = REAL(jk-1,wp) /REAL(jpk-1,wp)*(rn_hc/hbatt(ji,jj)) z_gsigw3p1(ji,jj) = REAL(jk,wp) /REAL(jpk-1,wp)*(rn_hc/hbatt(ji,jj)) z_gsigt3(ji,jj) = (REAL(jk-1,wp)+0.5)/REAL(jpk-1,wp)*(rn_hc/hbatt(ji,jj)) ENDIF !gsi3w3m1 & z_gsigt3m1 only used if jk /= 1 and is set at the end of the loop over jk IF( jk .EQ. 1) THEN z_esigw3(ji,jj ) = 2.0 * (z_gsigt3(ji,jj ) - z_gsigw3(ji,jj )) z_gsi3w3(ji,jj) = 0.5 * z_esigw3(ji,jj) ELSE z_esigw3(ji,jj) = z_gsigt3(ji,jj) - z_gsigt3m1(ji,jj) z_gsi3w3(ji,jj) = z_gsi3w3m1(ji,jj) + z_esigw3(ji,jj) ENDIF IF( jk .EQ. jpk) THEN z_esigt3(ji,jj) = 2.0 * (z_gsigt3(ji,jj) - z_gsigw3(ji,jj)) ELSE z_esigt3(ji,jj) = z_gsigw3p1(ji,jj) - z_gsigw3(ji,jj) ENDIF gdept_0(ji,jj) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsigt3(ji,jj) gdepw_0(ji,jj) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsigw3(ji,jj) gdep3w_0(ji,jj) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsi3w3(ji,jj) ENDDO ! for all jj's ENDDO ! for all ji's DO ji=1,jpi-1 DO jj=1,jpj-1 z_esigtu3(ji,jj) = ( hbatt(ji,jj)*z_esigt3(ji,jj)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj) ) / & ( hbatt(ji,jj)+hbatt(ji+1,jj) ) z_esigtv3(ji,jj) = ( hbatt(ji,jj)*z_esigt3(ji,jj)+hbatt(ji,jj+1)*z_esigt3(ji,jj+1) ) / & ( hbatt(ji,jj)+hbatt(ji,jj+1) ) z_esigtf3(ji,jj) = ( hbatt(ji,jj)*z_esigt3(ji,jj)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj) + & hbatt(ji,jj+1)*z_esigt3(ji,jj+1)+hbatt(ji+1,jj+1)*z_esigt3(ji+1,jj+1) ) / & ( hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1) ) z_esigwu3(ji,jj) = ( hbatt(ji,jj)*z_esigw3(ji,jj)+hbatt(ji+1,jj)*z_esigw3(ji+1,jj) ) / & ( hbatt(ji,jj)+hbatt(ji+1,jj) ) z_esigwv3(ji,jj) = ( hbatt(ji,jj)*z_esigw3(ji,jj)+hbatt(ji,jj+1)*z_esigw3(ji,jj+1) ) / & ( hbatt(ji,jj)+hbatt(ji,jj+1) ) e3t_0(ji,jj)=(scosrf(ji,jj)+hbatt(ji,jj))*z_esigt3(ji,jj) e3u_0(ji,jj)=(scosrf(ji,jj)+hbatu(ji,jj))*z_esigtu3(ji,jj) e3v_0(ji,jj)=(scosrf(ji,jj)+hbatv(ji,jj))*z_esigtv3(ji,jj) e3f_0(ji,jj)=(scosrf(ji,jj)+hbatf(ji,jj))*z_esigtf3(ji,jj) ! e3w_0(ji,jj)=hbatt(ji,jj)*z_esigw3(ji,jj) e3uw_0(ji,jj)=hbatu(ji,jj)*z_esigwu3(ji,jj) e3vw_0(ji,jj)=hbatv(ji,jj)*z_esigwv3(ji,jj) ENDDO ENDDO ! Keep some arrays for next level z_gsigt3m1 = z_gsigt3 z_gsiw3m1 = z_gsiw3 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 WRITE(*,*) 'Writing level ',jk,' to file' CALL write_netcdf_vars(jk) WRITE(*,*) 'Written level ',jk,' to file' ENDDO ! End of loop over jk DEALLOCATE( z_gsigw3, z_gsigt3, z_gsi3w3, z_gsigw3p1, z_gsigt3m1, z_gsi3w3m1) 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 END DO !!gm: e3uw, e3vw can be suppressed (modif in dynzdf, dynzdf_iso, zdfbfr) (save 2 3D arrays) DO jk = 1, jpk DO jj = 1, jpj DO ji = 1, jpi zcoeft = ( REAL(jk,wp) - 0.5 ) / REAL(jpk-1,wp) zcoefw = ( REAL(jk,wp) - 1.0 ) / REAL(jpk-1,wp) gdept_0(ji,jj) = ( scosrf(ji,jj) + (hbatt(ji,jj)-hift(ji,jj))*z_gsigt(jk) + hift(ji,jj)*zcoeft ) gdepw_0(:,:) = ( scosrf(ji,jj) + (hbatt(ji,jj)-hift(ji,jj))*z_gsigw(jk) + hift(ji,jj)*zcoefw ) gdep3w_0(:,:) = ( scosrf(ji,jj) + (hbatt(ji,jj)-hift(ji,jj))*z_gsi3w(jk) + hift(ji,jj)*zcoeft ) e3t_0(ji,jj) = ( (hbatt(ji,jj)-hift(ji,jj))*z_esigt(jk) + hift(ji,jj)/REAL(jpk-1,wp) ) e3u_0(ji,jj) = ( (hbatu(ji,jj)-hifu(ji,jj))*z_esigt(jk) + hifu(ji,jj)/REAL(jpk-1,wp) ) e3v_0(ji,jj) = ( (hbatv(ji,jj)-hifv(ji,jj))*z_esigt(jk) + hifv(ji,jj)/REAL(jpk-1,wp) ) e3f_0(ji,jj) = ( (hbatf(ji,jj)-hiff(ji,jj))*z_esigt(jk) + hiff(ji,jj)/REAL(jpk-1,wp) ) ! e3w_0(ji,jj) = ( (hbatt(ji,jj)-hift(ji,jj))*z_esigw(jk) + hift(ji,jj)/REAL(jpk-1,wp) ) e3uw_0(ji,jj) = ( (hbatu(ji,jj)-hifu(ji,jj))*z_esigw(jk) + hifu(ji,jj)/REAL(jpk-1,wp) ) e3vw_0(ji,jj) = ( (hbatv(ji,jj)-hifv(ji,jj))*z_esigw(jk) + hifv(ji,jj)/REAL(jpk-1,wp) ) IF( scobot(ji,jj) >= gdept_0(ji,jj) ) mbathy(ji,jj) = MAX( 2, jk ) IF( scobot(ji,jj) == 0. ) mbathy(ji,jj) = 0 END DO END DO 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 CALL write_netcdf_vars(jk) ENDDO ! End of loop over jk 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 ! continuous "k" coordinate REAL(wp) :: p_gamma ! 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 p_gamma = za*(pk1*(1.0-pk1/2.0))+zb*pk1**3.0 + & & zx*( (rn_alpha+2.0)*pk1**(rn_alpha+1.0)- & & (rn_alpha+1.0)*pk1**(rn_alpha+2.0) ) p_gamma = p_gamma*psmth+pk1*(1.0-psmth) ! END FUNCTION fgamma