[5255] | 1 | PROGRAM SCOORD_GEN |
---|
| 2 | !!---------------------------------------------------------------------- |
---|
| 3 | !! *** PROGRAM scoord_gen *** |
---|
| 4 | !! |
---|
| 5 | !! ** Purpose : define the s-coordinate system |
---|
| 6 | !! |
---|
| 7 | !! ** Method : s-coordinate |
---|
| 8 | !! The depth of model levels is defined as the product of an |
---|
| 9 | !! analytical function by the local bathymetry, while the vertical |
---|
| 10 | !! scale factors are defined as the product of the first derivative |
---|
| 11 | !! of the analytical function by the bathymetry. |
---|
| 12 | !! (this solution save memory as depth and scale factors are not |
---|
| 13 | !! 3d fields) |
---|
| 14 | !! - Read bathymetry (in meters) at t-point and compute the |
---|
| 15 | !! bathymetry at u-, v-, and f-points. |
---|
| 16 | !! hbatu = mi( hbatt ) |
---|
| 17 | !! hbatv = mj( hbatt ) |
---|
| 18 | !! hbatf = mi( mj( hbatt ) ) |
---|
| 19 | !! - Compute z_gsigt, z_gsigw, z_esigt, z_esigw from an analytical |
---|
| 20 | !! function and its derivative given as function. |
---|
| 21 | !! z_gsigt(k) = fssig (k ) |
---|
| 22 | !! z_gsigw(k) = fssig (k-0.5) |
---|
| 23 | !! z_esigt(k) = fsdsig(k ) |
---|
| 24 | !! z_esigw(k) = fsdsig(k-0.5) |
---|
| 25 | !! Three options for stretching are give, and they can be modified |
---|
| 26 | !! following the users requirements. Nevertheless, the output as |
---|
| 27 | !! well as the way to compute the model levels and scale factors |
---|
| 28 | !! must be respected in order to insure second order accuracy |
---|
| 29 | !! schemes. |
---|
| 30 | !! |
---|
| 31 | !! The three methods for stretching available are: |
---|
| 32 | !! |
---|
| 33 | !! s_sh94 (Song and Haidvogel 1994) |
---|
| 34 | !! a sinh/tanh function that allows sigma and stretched sigma |
---|
| 35 | !! |
---|
| 36 | !! s_sf12 (Siddorn and Furner 2012?) |
---|
| 37 | !! allows the maintenance of fixed surface and or |
---|
| 38 | !! bottom cell resolutions (cf. geopotential coordinates) |
---|
| 39 | !! within an analytically derived stretched S-coordinate framework. |
---|
| 40 | !! |
---|
| 41 | !! s_tanh (Madec et al 1996) |
---|
| 42 | !! a cosh/tanh function that gives stretched coordinates |
---|
| 43 | !! |
---|
| 44 | !!---------------------------------------------------------------------- |
---|
| 45 | ! |
---|
[5257] | 46 | USE utils |
---|
[5255] | 47 | IMPLICIT NONE |
---|
| 48 | |
---|
| 49 | |
---|
| 50 | !!---------------------------------------------------------------------- |
---|
| 51 | ! |
---|
| 52 | ! |
---|
| 53 | OPEN( UNIT=numnam, FILE='namelist', FORM='FORMATTED', STATUS='OLD' ) |
---|
| 54 | READ( numnam, namzgr_sco ) |
---|
| 55 | CLOSE( numnam ) |
---|
[5257] | 56 | jpk = INT(rn_jpk) |
---|
[5255] | 57 | |
---|
[5257] | 58 | WRITE(*,*) |
---|
| 59 | WRITE(*,*) 'domzgr_sco : s-coordinate or hybrid z-s-coordinate' |
---|
| 60 | WRITE(*,*) '~~~~~~~~~~~' |
---|
| 61 | WRITE(*,*) ' Namelist namzgr_sco' |
---|
| 62 | WRITE(*,*) ' stretching coeffs ' |
---|
| 63 | WRITE(*,*) ' Number of levels rn_jpk = ',rn_jpk |
---|
| 64 | WRITE(*,*) ' maximum depth of s-bottom surface (>0) rn_sbot_max = ',rn_sbot_max |
---|
| 65 | WRITE(*,*) ' minimum depth of s-bottom surface (>0) rn_sbot_min = ',rn_sbot_min |
---|
| 66 | WRITE(*,*) ' Critical depth rn_hc = ',rn_hc |
---|
| 67 | WRITE(*,*) ' maximum cut-off r-value allowed rn_rmax = ',rn_rmax |
---|
| 68 | WRITE(*,*) ' Song and Haidvogel 1994 stretching ln_s_sh94 = ',ln_s_sh94 |
---|
| 69 | WRITE(*,*) ' Song and Haidvogel 1994 stretching coefficients' |
---|
| 70 | WRITE(*,*) ' surface control parameter (0<=rn_theta<=20) rn_theta = ',rn_theta |
---|
| 71 | WRITE(*,*) ' bottom control parameter (0<=rn_thetb<= 1) rn_thetb = ',rn_thetb |
---|
| 72 | WRITE(*,*) ' stretching parameter (song and haidvogel) rn_bb = ',rn_bb |
---|
| 73 | WRITE(*,*) ' Siddorn and Furner 2012 stretching ln_s_sf12 = ',ln_s_sf12 |
---|
| 74 | WRITE(*,*) ' switching to sigma (T) or Z (F) at H<Hc ln_sigcrit = ',ln_sigcrit |
---|
| 75 | WRITE(*,*) ' Siddorn and Furner 2012 stretching coefficients' |
---|
| 76 | WRITE(*,*) ' stretchin parameter ( >1 surface; <1 bottom) rn_alpha = ',rn_alpha |
---|
| 77 | WRITE(*,*) ' e-fold length scale for transition region rn_efold = ',rn_efold |
---|
| 78 | WRITE(*,*) ' Surface cell depth (Zs) (m) rn_zs = ',rn_zs |
---|
| 79 | WRITE(*,*) ' Bathymetry multiplier for Zb rn_zb_a = ',rn_zb_a |
---|
| 80 | WRITE(*,*) ' Offset for Zb rn_zb_b = ',rn_zb_b |
---|
| 81 | WRITE(*,*) ' Bottom cell (Zb) (m) = H*rn_zb_a + rn_zb_b' |
---|
[5255] | 82 | |
---|
| 83 | |
---|
| 84 | ! Read in bathy, jpi and jpj from bathy.nc |
---|
| 85 | CALL read_bathy() |
---|
| 86 | |
---|
| 87 | !Allocate all other allocatable variables |
---|
| 88 | ios = dom_oce_alloc() |
---|
| 89 | IF( ios .NE. 0) THEN |
---|
| 90 | WRITE(0,*) 'Unable to allocate all arrays' |
---|
| 91 | STOP 1 |
---|
| 92 | ENDIF |
---|
| 93 | |
---|
| 94 | hift(:,:) = rn_sbot_min ! set the minimum depth for the s-coordinate |
---|
| 95 | hifu(:,:) = rn_sbot_min |
---|
| 96 | hifv(:,:) = rn_sbot_min |
---|
| 97 | hiff(:,:) = rn_sbot_min |
---|
| 98 | |
---|
| 99 | ! ! set maximum ocean depth |
---|
| 100 | bathy(:,:) = MIN( rn_sbot_max, bathy(:,:) ) |
---|
| 101 | |
---|
| 102 | DO jj = 1, jpj |
---|
| 103 | DO ji = 1, jpi |
---|
[5257] | 104 | IF( bathy(ji,jj) > 0. ) bathy(ji,jj) = MAX( rn_sbot_min, bathy(ji,jj) ) |
---|
[5255] | 105 | END DO |
---|
| 106 | END DO |
---|
| 107 | ! ! ============================= |
---|
| 108 | ! ! Define the envelop bathymetry (hbatt) |
---|
| 109 | ! ! ============================= |
---|
| 110 | ! use r-value to create hybrid coordinates |
---|
| 111 | zenv(:,:) = bathy(:,:) |
---|
| 112 | ! |
---|
| 113 | ! set first land point adjacent to a wet cell to sbot_min as this needs to be included in smoothing |
---|
| 114 | DO jj = 1, jpj |
---|
| 115 | DO ji = 1, jpi |
---|
[5257] | 116 | IF( bathy(ji,jj) == 0. ) THEN |
---|
[5255] | 117 | iip1 = MIN( ji+1, jpi ) |
---|
| 118 | ijp1 = MIN( jj+1, jpj ) |
---|
| 119 | iim1 = MAX( ji-1, 1 ) |
---|
| 120 | ijm1 = MAX( jj-1, 1 ) |
---|
| 121 | IF( (bathy(iip1,jj) + bathy(iim1,jj) + bathy(ji,ijp1) + bathy(ji,ijm1) + & |
---|
[5257] | 122 | & bathy(iip1,ijp1) + bathy(iim1,ijm1) + bathy(iip1,ijp1) + bathy(iim1,ijm1)) > 0. ) THEN |
---|
[5255] | 123 | zenv(ji,jj) = rn_sbot_min |
---|
| 124 | ENDIF |
---|
| 125 | ENDIF |
---|
| 126 | END DO |
---|
| 127 | END DO |
---|
| 128 | ! |
---|
| 129 | ! smooth the bathymetry (if required) |
---|
[5257] | 130 | scosrf(:,:) = 0. ! ocean surface depth (here zero: no under ice-shelf sea) |
---|
[5255] | 131 | scobot(:,:) = bathy(:,:) ! ocean bottom depth |
---|
| 132 | ! |
---|
| 133 | jl = 0 |
---|
[5257] | 134 | zrmax = 1. |
---|
[5255] | 135 | ! |
---|
| 136 | ! |
---|
| 137 | ! set scaling factor used in reducing vertical gradients |
---|
[5257] | 138 | zrfact = ( 1. - rn_rmax ) / ( 1. + rn_rmax ) |
---|
[5255] | 139 | ! |
---|
| 140 | ! initialise temporary evelope depth arrays |
---|
| 141 | ztmpi1(:,:) = zenv(:,:) |
---|
| 142 | ztmpi2(:,:) = zenv(:,:) |
---|
| 143 | ztmpj1(:,:) = zenv(:,:) |
---|
| 144 | ztmpj2(:,:) = zenv(:,:) |
---|
| 145 | ! |
---|
| 146 | ! initialise temporary r-value arrays |
---|
[5257] | 147 | zri(:,:) = 1. |
---|
| 148 | zrj(:,:) = 1. |
---|
[5255] | 149 | ! ! ================ ! |
---|
[5257] | 150 | DO WHILE( jl <= 10000 .AND. ( zrmax - rn_rmax ) > 1.e-8 ) ! Iterative loop ! |
---|
[5255] | 151 | ! ! ================ ! |
---|
| 152 | jl = jl + 1 |
---|
[5257] | 153 | zrmax = 0. |
---|
[5255] | 154 | ! we set zrmax from previous r-values (zri and zrj) first |
---|
| 155 | ! if set after current r-value calculation (as previously) |
---|
| 156 | ! we could exit DO WHILE prematurely before checking r-value |
---|
| 157 | ! of current zenv |
---|
[5257] | 158 | ! DO jj = 1, nlcj |
---|
| 159 | ! DO ji = 1, nlci |
---|
| 160 | DO jj = 1, jpi !jpi or jpim1? |
---|
| 161 | DO ji = 1, jpj |
---|
[5255] | 162 | zrmax = MAX( zrmax, ABS(zri(ji,jj)), ABS(zrj(ji,jj)) ) |
---|
| 163 | END DO |
---|
| 164 | END DO |
---|
[5257] | 165 | zri(:,:) = 0. |
---|
| 166 | zrj(:,:) = 0. |
---|
| 167 | ! DO jj = 1, nlci |
---|
| 168 | ! DO ji = 1, nlcj |
---|
| 169 | ! iip1 = MIN( ji+1, nlci ) ! force zri = 0 on last line (ji=ncli+1 to jpi) |
---|
| 170 | ! ijp1 = MIN( jj+1, nlcj ) ! force zrj = 0 on last raw (jj=nclj+1 to jpj) |
---|
| 171 | ! IF( (zenv(ji,jj) > 0.) .AND. (zenv(iip1,jj) > 0.)) THEN |
---|
| 172 | ! zri(ji,jj) = ( zenv(iip1,jj ) - zenv(ji,jj) ) / ( zenv(iip1,jj ) + zenv(ji,jj) ) |
---|
| 173 | ! END IF |
---|
| 174 | ! IF( (zenv(ji,jj) > 0.) .AND. (zenv(ji,ijp1) > 0.)) THEN |
---|
| 175 | ! zrj(ji,jj) = ( zenv(ji ,ijp1) - zenv(ji,jj) ) / ( zenv(ji ,ijp1) + zenv(ji,jj) ) |
---|
| 176 | ! END IF |
---|
| 177 | ! IF( zri(ji,jj) > rn_rmax ) ztmpi1(ji ,jj ) = zenv(iip1,jj ) * zrfact |
---|
| 178 | ! IF( zri(ji,jj) < -rn_rmax ) ztmpi2(iip1,jj ) = zenv(ji ,jj ) * zrfact |
---|
| 179 | ! IF( zrj(ji,jj) > rn_rmax ) ztmpj1(ji ,jj ) = zenv(ji ,ijp1) * zrfact |
---|
| 180 | ! IF( zrj(ji,jj) < -rn_rmax ) ztmpj2(ji ,ijp1) = zenv(ji ,jj ) * zrfact |
---|
| 181 | ! END DO |
---|
| 182 | ! END DO |
---|
| 183 | DO jj = 1, jpi-1 |
---|
| 184 | DO ji = 1, jpj-1 |
---|
| 185 | iip1 = ji+1 |
---|
| 186 | ijp1 = jj+1 |
---|
| 187 | IF( (zenv(ji,jj) > 0.) .AND. (zenv(iip1,jj) > 0.)) THEN |
---|
[5255] | 188 | zri(ji,jj) = ( zenv(iip1,jj ) - zenv(ji,jj) ) / ( zenv(iip1,jj ) + zenv(ji,jj) ) |
---|
| 189 | END IF |
---|
[5257] | 190 | IF( (zenv(ji,jj) > 0.) .AND. (zenv(ji,ijp1) > 0.)) THEN |
---|
[5255] | 191 | zrj(ji,jj) = ( zenv(ji ,ijp1) - zenv(ji,jj) ) / ( zenv(ji ,ijp1) + zenv(ji,jj) ) |
---|
| 192 | END IF |
---|
| 193 | IF( zri(ji,jj) > rn_rmax ) ztmpi1(ji ,jj ) = zenv(iip1,jj ) * zrfact |
---|
| 194 | IF( zri(ji,jj) < -rn_rmax ) ztmpi2(iip1,jj ) = zenv(ji ,jj ) * zrfact |
---|
| 195 | IF( zrj(ji,jj) > rn_rmax ) ztmpj1(ji ,jj ) = zenv(ji ,ijp1) * zrfact |
---|
| 196 | IF( zrj(ji,jj) < -rn_rmax ) ztmpj2(ji ,ijp1) = zenv(ji ,jj ) * zrfact |
---|
| 197 | END DO |
---|
| 198 | END DO |
---|
| 199 | ! |
---|
| 200 | WRITE(*,*) 'zgr_sco : iter= ',jl, ' rmax= ', zrmax |
---|
| 201 | ! |
---|
[5257] | 202 | DO jj = 1, jpi |
---|
| 203 | DO ji = 1, jpj |
---|
[5255] | 204 | zenv(ji,jj) = MAX(zenv(ji,jj), ztmpi1(ji,jj), ztmpi2(ji,jj), ztmpj1(ji,jj), ztmpj2(ji,jj) ) |
---|
| 205 | END DO |
---|
| 206 | END DO |
---|
| 207 | ! ! ================ ! |
---|
| 208 | END DO ! End loop ! |
---|
| 209 | ! ! ================ ! |
---|
| 210 | DO jj = 1, jpj |
---|
| 211 | DO ji = 1, jpi |
---|
| 212 | zenv(ji,jj) = MAX( zenv(ji,jj), rn_sbot_min ) ! set all points to avoid undefined scale value warnings |
---|
| 213 | END DO |
---|
| 214 | END DO |
---|
| 215 | ! |
---|
| 216 | ! Envelope bathymetry saved in hbatt |
---|
[5257] | 217 | ! TODO - get this section to work |
---|
[5255] | 218 | hbatt(:,:) = zenv(:,:) |
---|
[5257] | 219 | ! IF( MINVAL( gphit(:,:) ) * MAXVAL( gphit(:,:) ) <= 0. ) THEN |
---|
| 220 | ! CALL ctl_warn( ' s-coordinates are tapered in vicinity of the Equator' ) |
---|
| 221 | ! DO jj = 1, jpj |
---|
| 222 | ! DO ji = 1, jpi |
---|
| 223 | ! ztaper = EXP( -(gphit(ji,jj)/8.)**2. ) |
---|
| 224 | ! hbatt(ji,jj) = rn_sbot_max * ztaper + hbatt(ji,jj) * ( 1. - ztaper ) |
---|
| 225 | ! END DO |
---|
| 226 | ! END DO |
---|
| 227 | ! ENDIF |
---|
[5255] | 228 | ! |
---|
| 229 | ! ! ============================== |
---|
| 230 | ! ! hbatu, hbatv, hbatf fields |
---|
| 231 | ! ! ============================== |
---|
[5257] | 232 | WRITE(*,*) |
---|
| 233 | WRITE(*,*) ' zgr_sco: minimum depth of the envelop topography set to : ', rn_sbot_min |
---|
[5255] | 234 | hbatu(:,:) = rn_sbot_min |
---|
| 235 | hbatv(:,:) = rn_sbot_min |
---|
| 236 | hbatf(:,:) = rn_sbot_min |
---|
[5257] | 237 | DO jj = 1, jpj-1 |
---|
| 238 | DO ji = 1, jpi-1 ! NO vector opt. |
---|
| 239 | hbatu(ji,jj) = 0.50 * ( hbatt(ji ,jj) + hbatt(ji+1,jj ) ) |
---|
| 240 | hbatv(ji,jj) = 0.50 * ( hbatt(ji ,jj) + hbatt(ji ,jj+1) ) |
---|
| 241 | hbatf(ji,jj) = 0.25 * ( hbatt(ji ,jj) + hbatt(ji ,jj+1) & |
---|
[5255] | 242 | & + hbatt(ji+1,jj) + hbatt(ji+1,jj+1) ) |
---|
| 243 | END DO |
---|
| 244 | END DO |
---|
| 245 | ! |
---|
| 246 | zhbat(:,:) = hbatu(:,:) |
---|
| 247 | DO jj = 1, jpj |
---|
| 248 | DO ji = 1, jpi |
---|
[5257] | 249 | IF( hbatu(ji,jj) == 0. ) THEN |
---|
| 250 | IF( zhbat(ji,jj) == 0. ) hbatu(ji,jj) = rn_sbot_min |
---|
| 251 | IF( zhbat(ji,jj) /= 0. ) hbatu(ji,jj) = zhbat(ji,jj) |
---|
[5255] | 252 | ENDIF |
---|
| 253 | END DO |
---|
| 254 | END DO |
---|
| 255 | zhbat(:,:) = hbatv(:,:) |
---|
| 256 | DO jj = 1, jpj |
---|
| 257 | DO ji = 1, jpi |
---|
[5257] | 258 | IF( hbatv(ji,jj) == 0. ) THEN |
---|
| 259 | IF( zhbat(ji,jj) == 0. ) hbatv(ji,jj) = rn_sbot_min |
---|
| 260 | IF( zhbat(ji,jj) /= 0. ) hbatv(ji,jj) = zhbat(ji,jj) |
---|
[5255] | 261 | ENDIF |
---|
| 262 | END DO |
---|
| 263 | END DO |
---|
| 264 | zhbat(:,:) = hbatf(:,:) |
---|
| 265 | DO jj = 1, jpj |
---|
| 266 | DO ji = 1, jpi |
---|
[5257] | 267 | IF( hbatf(ji,jj) == 0. ) THEN |
---|
| 268 | IF( zhbat(ji,jj) == 0. ) hbatf(ji,jj) = rn_sbot_min |
---|
| 269 | IF( zhbat(ji,jj) /= 0. ) hbatf(ji,jj) = zhbat(ji,jj) |
---|
[5255] | 270 | ENDIF |
---|
| 271 | END DO |
---|
| 272 | END DO |
---|
| 273 | |
---|
| 274 | !!bug: key_helsinki a verifer |
---|
| 275 | hift(:,:) = MIN( hift(:,:), hbatt(:,:) ) |
---|
| 276 | hifu(:,:) = MIN( hifu(:,:), hbatu(:,:) ) |
---|
| 277 | hifv(:,:) = MIN( hifv(:,:), hbatv(:,:) ) |
---|
| 278 | hiff(:,:) = MIN( hiff(:,:), hbatf(:,:) ) |
---|
| 279 | |
---|
[5257] | 280 | WRITE(*,*) ' MAX val hif t ', MAXVAL( hift (:,:) ), ' f ', MAXVAL( hiff (:,:) ), & |
---|
[5255] | 281 | & ' u ', MAXVAL( hifu (:,:) ), ' v ', MAXVAL( hifv (:,:) ) |
---|
[5257] | 282 | WRITE(*,*) ' MIN val hif t ', MINVAL( hift (:,:) ), ' f ', MINVAL( hiff (:,:) ), & |
---|
[5255] | 283 | & ' u ', MINVAL( hifu (:,:) ), ' v ', MINVAL( hifv (:,:) ) |
---|
[5257] | 284 | WRITE(*,*) ' MAX val hbat t ', MAXVAL( hbatt(:,:) ), ' f ', MAXVAL( hbatf(:,:) ), & |
---|
[5255] | 285 | & ' u ', MAXVAL( hbatu(:,:) ), ' v ', MAXVAL( hbatv(:,:) ) |
---|
[5257] | 286 | WRITE(*,*) ' MIN val hbat t ', MINVAL( hbatt(:,:) ), ' f ', MINVAL( hbatf(:,:) ), & |
---|
[5255] | 287 | & ' u ', MINVAL( hbatu(:,:) ), ' v ', MINVAL( hbatv(:,:) ) |
---|
[5269] | 288 | |
---|
| 289 | ! Create the output file |
---|
| 290 | CALL make_coord_file() |
---|
[5255] | 291 | |
---|
| 292 | ! ! ======================= |
---|
| 293 | ! ! s-coordinate fields (gdep., e3.) |
---|
| 294 | ! ! ======================= |
---|
| 295 | ! |
---|
| 296 | ! non-dimensional "sigma" for model level depth at w- and t-levels |
---|
| 297 | |
---|
[5269] | 298 | |
---|
[5255] | 299 | !======================================================================== |
---|
| 300 | ! Song and Haidvogel 1994 (ln_s_sh94=T) |
---|
| 301 | ! Siddorn and Furner 2012 (ln_sf12=T) |
---|
| 302 | ! or tanh function (both false) |
---|
[5269] | 303 | ! To reduce memory loop over jpk and write each level to file |
---|
[5255] | 304 | !======================================================================== |
---|
| 305 | IF ( ln_s_sh94 ) THEN |
---|
| 306 | CALL s_sh94() |
---|
| 307 | ELSE IF ( ln_s_sf12 ) THEN |
---|
| 308 | CALL s_sf12() |
---|
| 309 | ELSE |
---|
| 310 | CALL s_tanh() |
---|
| 311 | ENDIF |
---|
| 312 | |
---|
[5269] | 313 | CALL check_nf90( nf90_put_var(ncout, var_ids(11), mbathy) ) |
---|
| 314 | CALL check_nf90( nf90_close(ncout) ) |
---|
[5255] | 315 | |
---|
[5269] | 316 | |
---|
[5255] | 317 | ! |
---|
[5257] | 318 | END PROGRAM SCOORD_GEN |
---|
[5255] | 319 | |
---|
| 320 | !!====================================================================== |
---|
| 321 | SUBROUTINE s_sh94() |
---|
| 322 | |
---|
| 323 | !!---------------------------------------------------------------------- |
---|
| 324 | !! *** ROUTINE s_sh94 *** |
---|
| 325 | !! |
---|
| 326 | !! ** Purpose : stretch the s-coordinate system |
---|
| 327 | !! |
---|
| 328 | !! ** Method : s-coordinate stretch using the Song and Haidvogel 1994 |
---|
| 329 | !! mixed S/sigma coordinate |
---|
| 330 | !! |
---|
| 331 | !! Reference : Song and Haidvogel 1994. |
---|
| 332 | !!---------------------------------------------------------------------- |
---|
| 333 | ! |
---|
[5257] | 334 | USE utils |
---|
[5255] | 335 | REAL(wp) :: zcoeft, zcoefw ! temporary scalars |
---|
| 336 | ! |
---|
[5269] | 337 | REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: z_gsigw3, z_gsigt3, z_gsi3w3 |
---|
| 338 | REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: z_gsigw3p1, z_gsigt3m1, z_gsi3w3m1 |
---|
| 339 | REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: z_esigt3, z_esigw3, z_esigtu3 |
---|
| 340 | REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 |
---|
[5255] | 341 | |
---|
[5269] | 342 | ALLOCATE( z_gsigw3(jpi,jpj), z_gsigt3(jpi,jpj), z_gsi3w3(jpi,jpj) ) |
---|
| 343 | ALLOCATE( z_esigt3(jpi,jpj), z_esigw3(jpi,jpj), z_esigtu3(jpi,jpj)) |
---|
| 344 | ALLOCATE( z_esigtv3(jpi,jpj), z_esigtf3(jpi,jpj), z_esigwu3(jpi,jpj)) |
---|
| 345 | ALLOCATE( z_esigwv3(jpi,jpj), z_gsigw3p1(jpi,jpj), z_gsigt3m1(jpi,jpj) ) |
---|
| 346 | ALLOCATE( z_gsi3w3m1(jpi,jpj) ) |
---|
| 347 | |
---|
[5257] | 348 | z_gsigw3 = 0. ; z_gsigt3 = 0. ; z_gsi3w3 = 0. |
---|
| 349 | z_esigt3 = 0. ; z_esigw3 = 0. |
---|
[5269] | 350 | z_esigt3p1 = 0. ; z_esigw3p1 = 0. |
---|
[5257] | 351 | z_esigtu3 = 0. ; z_esigtv3 = 0. ; z_esigtf3 = 0. |
---|
| 352 | z_esigwu3 = 0. ; z_esigwv3 = 0. |
---|
[5255] | 353 | |
---|
[5269] | 354 | DO jk = 1,jpk |
---|
| 355 | DO ji = 1, jpi |
---|
| 356 | DO jj = 1, jpj |
---|
[5255] | 357 | |
---|
[5269] | 358 | IF( hbatt(ji,jj) > rn_hc ) THEN !deep water, stretched sigma |
---|
| 359 | z_gsigw3(ji,jj) = -fssig1( REAL(jk,wp)-0.5, rn_bb ) |
---|
| 360 | z_gsigw3p1(ji,jj) = -fssig1( REAL(jk+1,wp)-0.5, rn_bb ) |
---|
| 361 | z_gsigt3(ji,jj) = -fssig1( REAL(jk,wp) , rn_bb ) |
---|
| 362 | ELSE ! shallow water, uniform sigma |
---|
| 363 | z_gsigw3(ji,jj) = REAL(jk-1,wp) / REAL(jpk-1,wp) |
---|
| 364 | z_gsigw3p1(ji,jj) = REAL(jk,wp) / REAL(jpk-1,wp) |
---|
| 365 | z_gsigt3(ji,jj) = ( REAL(jk-1,wp) + 0.5 ) / REAL(jpk-1,wp) |
---|
| 366 | ENDIF |
---|
| 367 | ! |
---|
| 368 | !gsi3w3m1 & gsigt3m1 only used if jk /= 1 and is set at the end of the loop over jk |
---|
| 369 | IF( jk .EQ. 1) THEN |
---|
| 370 | z_esigw3(ji,jj ) = 2. * ( z_gsigt3(ji,jj ) - z_gsigw3(ji,jj ) ) |
---|
| 371 | z_gsi3w3(ji,jj) = 0.5 * z_esigw3(ji,jj) |
---|
| 372 | ELSE |
---|
| 373 | z_esigw3(ji,jj) = z_gsigt3(ji,jj) - z_gsigt3m1(ji,jj) |
---|
| 374 | z_gsi3w3(ji,jj) = z_gsi3w3m1(ji,jj) + z_esigw3(ji,jj) |
---|
| 375 | ENDIF |
---|
| 376 | IF( jk .EQ. jpk) THEN |
---|
| 377 | z_esigt3(ji,jj) = 2. * ( z_gsigt3(ji,jj) - z_gsigw3(ji,jj) ) |
---|
| 378 | ELSE |
---|
| 379 | z_esigt3(ji,jj ) = z_gsigw3p1(ji,jj) - z_gsigw3(ji,jj) |
---|
| 380 | ENDIF |
---|
| 381 | ! |
---|
| 382 | |
---|
[5257] | 383 | zcoeft = ( REAL(jk,wp) - 0.5 ) / REAL(jpk-1,wp) |
---|
| 384 | zcoefw = ( REAL(jk,wp) - 1.0 ) / REAL(jpk-1,wp) |
---|
[5269] | 385 | gdept_0(ji,jj) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*z_gsigt3(ji,jj)+rn_hc*zcoeft ) |
---|
| 386 | gdepw_0(ji,jj) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*z_gsigw3(ji,jj)+rn_hc*zcoefw ) |
---|
| 387 | gdep3w_0(ji,jj) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*z_gsi3w3(ji,jj)+rn_hc*zcoeft ) |
---|
| 388 | ! |
---|
| 389 | END DO ! for all jj's |
---|
| 390 | END DO ! for all ji's |
---|
[5255] | 391 | |
---|
[5269] | 392 | DO ji = 1, jpim1 |
---|
| 393 | DO jj = 1, jpjm1 |
---|
| 394 | z_esigtu3(ji,jj) = ( hbatt(ji,jj)*z_esigt3(ji,jj)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj) ) & |
---|
[5255] | 395 | & / ( hbatt(ji,jj)+hbatt(ji+1,jj) ) |
---|
[5269] | 396 | z_esigtv3(ji,jj) = ( hbatt(ji,jj)*z_esigt3(ji,jj)+hbatt(ji,jj+1)*z_esigt3(ji,jj+1) ) & |
---|
[5255] | 397 | & / ( hbatt(ji,jj)+hbatt(ji,jj+1) ) |
---|
[5269] | 398 | z_esigtf3(ji,jj) = ( hbatt(ji,jj)*z_esigt3(ji,jj)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj) & |
---|
| 399 | & + hbatt(ji,jj+1)*z_esigt3(ji,jj+1)+hbatt(ji+1,jj+1)*z_esigt3(ji+1,jj+1) ) & |
---|
[5255] | 400 | & / ( hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1) ) |
---|
[5269] | 401 | z_esigwu3(ji,jj) = ( hbatt(ji,jj)*z_esigw3(ji,jj)+hbatt(ji+1,jj)*z_esigw3(ji+1,jj) ) & |
---|
[5255] | 402 | & / ( hbatt(ji,jj)+hbatt(ji+1,jj) ) |
---|
[5269] | 403 | z_esigwv3(ji,jj) = ( hbatt(ji,jj)*z_esigw3(ji,jj)+hbatt(ji,jj+1)*z_esigw3(ji,jj+1) ) & |
---|
[5255] | 404 | & / ( hbatt(ji,jj)+hbatt(ji,jj+1) ) |
---|
| 405 | ! |
---|
[5269] | 406 | e3t_0(ji,jj) = ( (hbatt(ji,jj)-rn_hc)*z_esigt3 (ji,jj) + rn_hc/REAL(jpk-1,wp) ) |
---|
| 407 | e3u_0(ji,jj) = ( (hbatu(ji,jj)-rn_hc)*z_esigtu3(ji,jj) + rn_hc/REAL(jpk-1,wp) ) |
---|
| 408 | e3v_0(ji,jj) = ( (hbatv(ji,jj)-rn_hc)*z_esigtv3(ji,jj) + rn_hc/REAL(jpk-1,wp) ) |
---|
| 409 | e3f_0(ji,jj) = ( (hbatf(ji,jj)-rn_hc)*z_esigtf3(ji,jj) + rn_hc/REAL(jpk-1,wp) ) |
---|
[5255] | 410 | ! |
---|
[5269] | 411 | e3w_0 (ji,jj) = ( (hbatt(ji,jj)-rn_hc)*z_esigw3 (ji,jj) + rn_hc/REAL(jpk-1,wp) ) |
---|
| 412 | e3uw_0(ji,jj) = ( (hbatu(ji,jj)-rn_hc)*z_esigwu3(ji,jj) + rn_hc/REAL(jpk-1,wp) ) |
---|
| 413 | e3vw_0(ji,jj) = ( (hbatv(ji,jj)-rn_hc)*z_esigwv3(ji,jj) + rn_hc/REAL(jpk-1,wp) ) |
---|
[5255] | 414 | END DO |
---|
[5269] | 415 | END DO |
---|
[5255] | 416 | |
---|
[5269] | 417 | z_gsigt3m1 = z_gsigt3 |
---|
| 418 | z_gsiw3m1 = z_gsiw3 |
---|
| 419 | |
---|
| 420 | where (e3t_0 (:,:).eq.0.0) e3t_0(:,:) = 1.0 |
---|
| 421 | where (e3u_0 (:,:).eq.0.0) e3u_0(:,:) = 1.0 |
---|
| 422 | where (e3v_0 (:,:).eq.0.0) e3v_0(:,:) = 1.0 |
---|
| 423 | where (e3f_0 (:,:).eq.0.0) e3f_0(:,:) = 1.0 |
---|
| 424 | where (e3w_0 (:,:).eq.0.0) e3w_0(:,:) = 1.0 |
---|
| 425 | where (e3uw_0 (:,:).eq.0.0) e3uw_0(:,:) = 1.0 |
---|
| 426 | where (e3vw_0 (:,:).eq.0.0) e3vw_0(:,:) = 1.0 |
---|
| 427 | |
---|
| 428 | CALL write_netcdf_vars(jk) |
---|
| 429 | DO jj = 1, jpj |
---|
| 430 | DO ji = 1, jpi |
---|
| 431 | IF( scobot(ji,jj) >= gdept_0(ji,jj) ) mbathy(ji,jj) = MAX( 2, jk ) |
---|
| 432 | IF( scobot(ji,jj) == 0. ) mbathy(ji,jj) = 0 |
---|
| 433 | END DO |
---|
| 434 | END DO |
---|
| 435 | END DO !End of loop over jk |
---|
| 436 | |
---|
| 437 | DEALLOCATE( z_gsigw3, z_gsigt3, z_gsi3w3, z_gsigw3p1, z_gsigt3m1, z_gsi3w3m1) |
---|
[5255] | 438 | DEALLOCATE( z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 ) |
---|
| 439 | |
---|
| 440 | END SUBROUTINE s_sh94 |
---|
| 441 | |
---|
| 442 | SUBROUTINE s_sf12 |
---|
| 443 | |
---|
| 444 | !!---------------------------------------------------------------------- |
---|
| 445 | !! *** ROUTINE s_sf12 *** |
---|
| 446 | !! |
---|
| 447 | !! ** Purpose : stretch the s-coordinate system |
---|
| 448 | !! |
---|
| 449 | !! ** Method : s-coordinate stretch using the Siddorn and Furner 2012? |
---|
| 450 | !! mixed S/sigma/Z coordinate |
---|
| 451 | !! |
---|
| 452 | !! This method allows the maintenance of fixed surface and or |
---|
| 453 | !! bottom cell resolutions (cf. geopotential coordinates) |
---|
| 454 | !! within an analytically derived stretched S-coordinate framework. |
---|
| 455 | !! |
---|
| 456 | !! |
---|
| 457 | !! Reference : Siddorn and Furner 2012 (submitted Ocean modelling). |
---|
| 458 | !!---------------------------------------------------------------------- |
---|
| 459 | ! |
---|
[5257] | 460 | USE utils |
---|
[5255] | 461 | REAL(wp) :: zsmth ! smoothing around critical depth |
---|
| 462 | REAL(wp) :: zzs, zzb ! Surface and bottom cell thickness in sigma space |
---|
| 463 | ! |
---|
[5269] | 464 | REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: z_gsigw3, z_gsigt3, z_gsi3w3 |
---|
| 465 | REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: z_gsigw3p1, z_gsigt3m1, z_gsi3w3m1 |
---|
| 466 | REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: z_esigt3, z_esigw3, z_esigtu3 |
---|
| 467 | REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 |
---|
[5255] | 468 | |
---|
[5269] | 469 | ALLOCATE( z_gsigw3(jpi,jpj), z_gsigt3(jpi,jpj), z_gsi3w3(jpi,jpj) ) |
---|
| 470 | ALLOCATE( z_esigt3(jpi,jpj), z_esigw3(jpi,jpj), z_esigtu3(jpi,jpj)) |
---|
| 471 | ALLOCATE( z_esigtv3(jpi,jpj), z_esigtf3(jpi,jpj), z_esigwu3(jpi,jpj)) |
---|
| 472 | ALLOCATE( z_esigwv3(jpi,jpj), z_gsigw3p1(jpi,jpj), z_gsigt3m1(jpi,jpj) ) |
---|
| 473 | ALLOCATE( z_gsi3w3m1(jpi,jpj) ) |
---|
| 474 | |
---|
[5257] | 475 | z_gsigw3 = 0. ; z_gsigt3 = 0. ; z_gsi3w3 = 0. |
---|
| 476 | z_esigt3 = 0. ; z_esigw3 = 0. |
---|
| 477 | z_esigtu3 = 0. ; z_esigtv3 = 0. ; z_esigtf3 = 0. |
---|
| 478 | z_esigwu3 = 0. ; z_esigwv3 = 0. |
---|
[5269] | 479 | |
---|
[5255] | 480 | |
---|
[5269] | 481 | DO jk = 1,jpk |
---|
| 482 | DO ji = 1, jpi |
---|
| 483 | DO jj = 1, jpj |
---|
| 484 | IF (hbatt(ji,jj)>rn_hc) THEN !deep water, stretched sigma |
---|
[5255] | 485 | |
---|
[5269] | 486 | zzb = hbatt(ji,jj)*rn_zb_a + rn_zb_b ! this forces a linear bottom cell depth relationship with H,. |
---|
[5255] | 487 | ! could be changed by users but care must be taken to do so carefully |
---|
[5269] | 488 | zzb = 1.0-(zzb/hbatt(ji,jj)) |
---|
[5255] | 489 | |
---|
[5269] | 490 | zzs = rn_zs / hbatt(ji,jj) |
---|
[5255] | 491 | |
---|
[5269] | 492 | IF (rn_efold /= 0.0) THEN |
---|
| 493 | zsmth = tanh( (hbatt(ji,jj)- rn_hc ) / rn_efold ) |
---|
| 494 | ELSE |
---|
| 495 | zsmth = 1.0 |
---|
| 496 | ENDIF |
---|
| 497 | |
---|
| 498 | z_gsigw3(ji,jj) = REAL(jk-1,wp) /REAL(jpk-1,wp) |
---|
| 499 | z_gsigw3p1(ji,jj) = REAL(jk,wp) /REAL(jpk-1,wp) |
---|
| 500 | z_gsigt3(ji,jj) = (REAL(jk-1,wp)+0.5)/REAL(jpk-1,wp) |
---|
| 501 | z_gsigw3(ji,jj) = fgamma( z_gsigw3(ji,jj), zzb, zzs, zsmth ) |
---|
| 502 | z_gsigw3p1(ji,jj) = fgamma( z_gsigw3p1(ji,jj), zzb, zzs, zsmth ) |
---|
| 503 | z_gsigt3(ji,jj) = fgamma( z_gsigt3(ji,jj), zzb, zzs, zsmth ) |
---|
[5255] | 504 | |
---|
[5269] | 505 | ELSE IF(ln_sigcrit) THEN ! shallow water, uniform sigma |
---|
[5255] | 506 | |
---|
[5269] | 507 | z_gsigw3(ji,jj) = REAL(jk-1,wp) /REAL(jpk-1,wp) |
---|
| 508 | z_gsigw3p1(ji,jj) = REAL(jk,wp) /REAL(jpk-1,wp) |
---|
| 509 | z_gsigt3(ji,jj) = (REAL(jk-1,wp)+0.5)/REAL(jpk-1,wp) |
---|
[5255] | 510 | |
---|
[5269] | 511 | ELSE ! shallow water, z coordinates |
---|
[5255] | 512 | |
---|
[5269] | 513 | z_gsigw3(ji,jj) = REAL(jk-1,wp) /REAL(jpk-1,wp)*(rn_hc/hbatt(ji,jj)) |
---|
| 514 | z_gsigw3p1(ji,jj) = REAL(jk,wp) /REAL(jpk-1,wp)*(rn_hc/hbatt(ji,jj)) |
---|
| 515 | z_gsigt3(ji,jj) = (REAL(jk-1,wp)+0.5)/REAL(jpk-1,wp)*(rn_hc/hbatt(ji,jj)) |
---|
[5255] | 516 | |
---|
[5269] | 517 | ENDIF |
---|
[5255] | 518 | |
---|
[5269] | 519 | !gsi3w3m1 & z_gsigt3m1 only used if jk /= 1 and is set at the end of the loop over jk |
---|
| 520 | IF( jk .EQ. 1) THEN |
---|
| 521 | z_esigw3(ji,jj ) = 2.0 * (z_gsigt3(ji,jj ) - z_gsigw3(ji,jj )) |
---|
| 522 | z_gsi3w3(ji,jj) = 0.5 * z_esigw3(ji,jj) |
---|
| 523 | ELSE |
---|
| 524 | z_esigw3(ji,jj) = z_gsigt3(ji,jj) - z_gsigt3m1(ji,jj) |
---|
| 525 | z_gsi3w3(ji,jj) = z_gsi3w3m1(ji,jj) + z_esigw3(ji,jj) |
---|
| 526 | ENDIF |
---|
| 527 | IF( jk .EQ. jpk) THEN |
---|
| 528 | z_esigt3(ji,jj) = 2.0 * (z_gsigt3(ji,jj) - z_gsigw3(ji,jj)) |
---|
| 529 | ELSE |
---|
| 530 | z_esigt3(ji,jj) = z_gsigw3p1(ji,jj) - z_gsigw3(ji,jj) |
---|
| 531 | ENDIF |
---|
[5255] | 532 | |
---|
[5269] | 533 | gdept_0(ji,jj) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsigt3(ji,jj) |
---|
| 534 | gdepw_0(ji,jj) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsigw3(ji,jj) |
---|
| 535 | gdep3w_0(ji,jj) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsi3w3(ji,jj) |
---|
[5255] | 536 | |
---|
[5269] | 537 | ENDDO ! for all jj's |
---|
| 538 | ENDDO ! for all ji's |
---|
[5255] | 539 | |
---|
[5269] | 540 | DO ji=1,jpi-1 |
---|
| 541 | DO jj=1,jpj-1 |
---|
[5255] | 542 | |
---|
[5269] | 543 | z_esigtu3(ji,jj) = ( hbatt(ji,jj)*z_esigt3(ji,jj)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj) ) / & |
---|
| 544 | ( hbatt(ji,jj)+hbatt(ji+1,jj) ) |
---|
| 545 | z_esigtv3(ji,jj) = ( hbatt(ji,jj)*z_esigt3(ji,jj)+hbatt(ji,jj+1)*z_esigt3(ji,jj+1) ) / & |
---|
| 546 | ( hbatt(ji,jj)+hbatt(ji,jj+1) ) |
---|
| 547 | z_esigtf3(ji,jj) = ( hbatt(ji,jj)*z_esigt3(ji,jj)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj) + & |
---|
| 548 | hbatt(ji,jj+1)*z_esigt3(ji,jj+1)+hbatt(ji+1,jj+1)*z_esigt3(ji+1,jj+1) ) / & |
---|
| 549 | ( hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1) ) |
---|
| 550 | z_esigwu3(ji,jj) = ( hbatt(ji,jj)*z_esigw3(ji,jj)+hbatt(ji+1,jj)*z_esigw3(ji+1,jj) ) / & |
---|
| 551 | ( hbatt(ji,jj)+hbatt(ji+1,jj) ) |
---|
| 552 | z_esigwv3(ji,jj) = ( hbatt(ji,jj)*z_esigw3(ji,jj)+hbatt(ji,jj+1)*z_esigw3(ji,jj+1) ) / & |
---|
| 553 | ( hbatt(ji,jj)+hbatt(ji,jj+1) ) |
---|
[5255] | 554 | |
---|
[5269] | 555 | e3t_0(ji,jj)=(scosrf(ji,jj)+hbatt(ji,jj))*z_esigt3(ji,jj) |
---|
| 556 | e3u_0(ji,jj)=(scosrf(ji,jj)+hbatu(ji,jj))*z_esigtu3(ji,jj) |
---|
| 557 | e3v_0(ji,jj)=(scosrf(ji,jj)+hbatv(ji,jj))*z_esigtv3(ji,jj) |
---|
| 558 | e3f_0(ji,jj)=(scosrf(ji,jj)+hbatf(ji,jj))*z_esigtf3(ji,jj) |
---|
| 559 | ! |
---|
| 560 | e3w_0(ji,jj)=hbatt(ji,jj)*z_esigw3(ji,jj) |
---|
| 561 | e3uw_0(ji,jj)=hbatu(ji,jj)*z_esigwu3(ji,jj) |
---|
| 562 | e3vw_0(ji,jj)=hbatv(ji,jj)*z_esigwv3(ji,jj) |
---|
| 563 | ENDDO |
---|
| 564 | ENDDO |
---|
| 565 | ! Keep some arrays for next level |
---|
| 566 | z_gsigt3m1 = z_gsigt3 |
---|
| 567 | z_gsiw3m1 = z_gsiw3 |
---|
| 568 | |
---|
| 569 | where (e3t_0 (:,:).eq.0.0) e3t_0(:,:) = 1.0 |
---|
| 570 | where (e3u_0 (:,:).eq.0.0) e3u_0(:,:) = 1.0 |
---|
| 571 | where (e3v_0 (:,:).eq.0.0) e3v_0(:,:) = 1.0 |
---|
| 572 | where (e3f_0 (:,:).eq.0.0) e3f_0(:,:) = 1.0 |
---|
| 573 | where (e3w_0 (:,:).eq.0.0) e3w_0(:,:) = 1.0 |
---|
| 574 | where (e3uw_0 (:,:).eq.0.0) e3uw_0(:,:) = 1.0 |
---|
| 575 | where (e3vw_0 (:,:).eq.0.0) e3vw_0(:,:) = 1.0 |
---|
| 576 | |
---|
| 577 | WRITE(*,*) 'Writing level ',jk,' to file' |
---|
| 578 | CALL write_netcdf_vars(jk) |
---|
| 579 | WRITE(*,*) 'Written level ',jk,' to file' |
---|
| 580 | ENDDO ! End of loop over jk |
---|
[5255] | 581 | |
---|
[5269] | 582 | DEALLOCATE( z_gsigw3, z_gsigt3, z_gsi3w3, z_gsigw3p1, z_gsigt3m1, z_gsi3w3m1) |
---|
[5255] | 583 | DEALLOCATE( z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 ) |
---|
| 584 | |
---|
| 585 | END SUBROUTINE s_sf12 |
---|
| 586 | |
---|
| 587 | SUBROUTINE s_tanh() |
---|
| 588 | |
---|
| 589 | !!---------------------------------------------------------------------- |
---|
| 590 | !! *** ROUTINE s_tanh*** |
---|
| 591 | !! |
---|
| 592 | !! ** Purpose : stretch the s-coordinate system |
---|
| 593 | !! |
---|
| 594 | !! ** Method : s-coordinate stretch |
---|
| 595 | !! |
---|
| 596 | !! Reference : Madec, Lott, Delecluse and Crepon, 1996. JPO, 26, 1393-1408. |
---|
| 597 | !!---------------------------------------------------------------------- |
---|
| 598 | |
---|
[5257] | 599 | USE utils |
---|
[5255] | 600 | REAL(wp) :: zcoeft, zcoefw ! temporary scalars |
---|
| 601 | |
---|
| 602 | REAL(wp), ALLOCATABLE, DIMENSION(:) :: z_gsigw, z_gsigt, z_gsi3w |
---|
| 603 | REAL(wp), ALLOCATABLE, DIMENSION(:) :: z_esigt, z_esigw |
---|
| 604 | |
---|
| 605 | ALLOCATE( z_gsigw(jpk), z_gsigt(jpk), z_gsi3w(jpk) ) |
---|
| 606 | ALLOCATE( z_esigt(jpk), z_esigw(jpk) ) |
---|
| 607 | |
---|
[5257] | 608 | z_gsigw = 0. ; z_gsigt = 0. ; z_gsi3w = 0. |
---|
| 609 | z_esigt = 0. ; z_esigw = 0. |
---|
[5255] | 610 | |
---|
| 611 | DO jk = 1, jpk |
---|
[5257] | 612 | z_gsigw(jk) = -fssig( REAL(jk,wp)-0.5 ) |
---|
[5255] | 613 | z_gsigt(jk) = -fssig( REAL(jk,wp) ) |
---|
| 614 | END DO |
---|
| 615 | WRITE(*,*) 'z_gsigw 1 jpk ', z_gsigw(1), z_gsigw(jpk) |
---|
| 616 | ! |
---|
| 617 | ! Coefficients for vertical scale factors at w-, t- levels |
---|
| 618 | !!gm bug : define it from analytical function, not like juste bellow.... |
---|
| 619 | !!gm or betteroffer the 2 possibilities.... |
---|
| 620 | DO jk = 1, jpk-1 |
---|
| 621 | z_esigt(jk ) = z_gsigw(jk+1) - z_gsigw(jk) |
---|
| 622 | z_esigw(jk+1) = z_gsigt(jk+1) - z_gsigt(jk) |
---|
| 623 | END DO |
---|
[5257] | 624 | z_esigw( 1 ) = 2. * ( z_gsigt(1 ) - z_gsigw(1 ) ) |
---|
| 625 | z_esigt(jpk) = 2. * ( z_gsigt(jpk) - z_gsigw(jpk) ) |
---|
[5255] | 626 | ! |
---|
| 627 | ! Coefficients for vertical depth as the sum of e3w scale factors |
---|
[5257] | 628 | z_gsi3w(1) = 0.5 * z_esigw(1) |
---|
[5255] | 629 | DO jk = 2, jpk |
---|
| 630 | z_gsi3w(jk) = z_gsi3w(jk-1) + z_esigw(jk) |
---|
| 631 | END DO |
---|
| 632 | !!gm: depuw, depvw can be suppressed (modif in ldfslp) and depw=dep3w can be set (save 3 3D arrays) |
---|
| 633 | DO jk = 1, jpk |
---|
| 634 | END DO |
---|
| 635 | !!gm: e3uw, e3vw can be suppressed (modif in dynzdf, dynzdf_iso, zdfbfr) (save 2 3D arrays) |
---|
[5269] | 636 | DO jk = 1, jpk |
---|
| 637 | DO jj = 1, jpj |
---|
| 638 | DO ji = 1, jpi |
---|
| 639 | zcoeft = ( REAL(jk,wp) - 0.5 ) / REAL(jpk-1,wp) |
---|
| 640 | zcoefw = ( REAL(jk,wp) - 1.0 ) / REAL(jpk-1,wp) |
---|
| 641 | gdept_0(ji,jj) = ( scosrf(ji,jj) + (hbatt(ji,jj)-hift(ji,jj))*z_gsigt(jk) + hift(ji,jj)*zcoeft ) |
---|
| 642 | gdepw_0(:,:) = ( scosrf(ji,jj) + (hbatt(ji,jj)-hift(ji,jj))*z_gsigw(jk) + hift(ji,jj)*zcoefw ) |
---|
| 643 | gdep3w_0(:,:) = ( scosrf(ji,jj) + (hbatt(ji,jj)-hift(ji,jj))*z_gsi3w(jk) + hift(ji,jj)*zcoeft ) |
---|
| 644 | e3t_0(ji,jj) = ( (hbatt(ji,jj)-hift(ji,jj))*z_esigt(jk) + hift(ji,jj)/REAL(jpk-1,wp) ) |
---|
| 645 | e3u_0(ji,jj) = ( (hbatu(ji,jj)-hifu(ji,jj))*z_esigt(jk) + hifu(ji,jj)/REAL(jpk-1,wp) ) |
---|
| 646 | e3v_0(ji,jj) = ( (hbatv(ji,jj)-hifv(ji,jj))*z_esigt(jk) + hifv(ji,jj)/REAL(jpk-1,wp) ) |
---|
| 647 | e3f_0(ji,jj) = ( (hbatf(ji,jj)-hiff(ji,jj))*z_esigt(jk) + hiff(ji,jj)/REAL(jpk-1,wp) ) |
---|
[5255] | 648 | ! |
---|
[5269] | 649 | e3w_0(ji,jj) = ( (hbatt(ji,jj)-hift(ji,jj))*z_esigw(jk) + hift(ji,jj)/REAL(jpk-1,wp) ) |
---|
| 650 | e3uw_0(ji,jj) = ( (hbatu(ji,jj)-hifu(ji,jj))*z_esigw(jk) + hifu(ji,jj)/REAL(jpk-1,wp) ) |
---|
| 651 | e3vw_0(ji,jj) = ( (hbatv(ji,jj)-hifv(ji,jj))*z_esigw(jk) + hifv(ji,jj)/REAL(jpk-1,wp) ) |
---|
| 652 | IF( scobot(ji,jj) >= gdept_0(ji,jj) ) mbathy(ji,jj) = MAX( 2, jk ) |
---|
| 653 | IF( scobot(ji,jj) == 0. ) mbathy(ji,jj) = 0 |
---|
[5255] | 654 | END DO |
---|
| 655 | END DO |
---|
[5269] | 656 | where (e3t_0 (:,:).eq.0.0) e3t_0(:,:) = 1.0 |
---|
| 657 | where (e3u_0 (:,:).eq.0.0) e3u_0(:,:) = 1.0 |
---|
| 658 | where (e3v_0 (:,:).eq.0.0) e3v_0(:,:) = 1.0 |
---|
| 659 | where (e3f_0 (:,:).eq.0.0) e3f_0(:,:) = 1.0 |
---|
| 660 | where (e3w_0 (:,:).eq.0.0) e3w_0(:,:) = 1.0 |
---|
| 661 | where (e3uw_0 (:,:).eq.0.0) e3uw_0(:,:) = 1.0 |
---|
| 662 | where (e3vw_0 (:,:).eq.0.0) e3vw_0(:,:) = 1.0 |
---|
| 663 | |
---|
| 664 | CALL write_netcdf_vars(jk) |
---|
| 665 | ENDDO ! End of loop over jk |
---|
[5255] | 666 | |
---|
[5269] | 667 | |
---|
[5255] | 668 | DEALLOCATE( z_gsigw, z_gsigt, z_gsi3w ) |
---|
| 669 | DEALLOCATE( z_esigt, z_esigw ) |
---|
| 670 | |
---|
| 671 | END SUBROUTINE s_tanh |
---|
| 672 | |
---|
| 673 | FUNCTION fssig( pk ) RESULT( pf ) |
---|
| 674 | !!---------------------------------------------------------------------- |
---|
| 675 | !! *** ROUTINE fssig *** |
---|
| 676 | !! |
---|
| 677 | !! ** Purpose : provide the analytical function in s-coordinate |
---|
| 678 | !! |
---|
| 679 | !! ** Method : the function provide the non-dimensional position of |
---|
| 680 | !! T and W (i.e. between 0 and 1) |
---|
| 681 | !! T-points at integer values (between 1 and jpk) |
---|
| 682 | !! W-points at integer values - 1/2 (between 0.5 and jpk-0.5) |
---|
| 683 | !!---------------------------------------------------------------------- |
---|
[5257] | 684 | USE utils, ONLY : wp |
---|
[5255] | 685 | REAL(wp), INTENT(in) :: pk ! continuous "k" coordinate |
---|
| 686 | REAL(wp) :: pf ! sigma value |
---|
| 687 | !!---------------------------------------------------------------------- |
---|
| 688 | ! |
---|
[5257] | 689 | pf = ( TANH( rn_theta * ( -(pk-0.5) / REAL(jpk-1) + rn_thetb ) ) & |
---|
[5255] | 690 | & - TANH( rn_thetb * rn_theta ) ) & |
---|
| 691 | & * ( COSH( rn_theta ) & |
---|
[5257] | 692 | & + COSH( rn_theta * ( 2. * rn_thetb - 1. ) ) ) & |
---|
| 693 | & / ( 2. * SINH( rn_theta ) ) |
---|
[5255] | 694 | ! |
---|
| 695 | END FUNCTION fssig |
---|
| 696 | |
---|
| 697 | |
---|
| 698 | FUNCTION fssig1( pk1, pbb ) RESULT( pf1 ) |
---|
| 699 | !!---------------------------------------------------------------------- |
---|
| 700 | !! *** ROUTINE fssig1 *** |
---|
| 701 | !! |
---|
| 702 | !! ** Purpose : provide the Song and Haidvogel version of the analytical function in s-coordinate |
---|
| 703 | !! |
---|
| 704 | !! ** Method : the function provides the non-dimensional position of |
---|
| 705 | !! T and W (i.e. between 0 and 1) |
---|
| 706 | !! T-points at integer values (between 1 and jpk) |
---|
| 707 | !! W-points at integer values - 1/2 (between 0.5 and jpk-0.5) |
---|
| 708 | !!---------------------------------------------------------------------- |
---|
[5257] | 709 | USE utils, ONLY : wp |
---|
[5255] | 710 | REAL(wp), INTENT(in) :: pk1 ! continuous "k" coordinate |
---|
| 711 | REAL(wp), INTENT(in) :: pbb ! Stretching coefficient |
---|
| 712 | REAL(wp) :: pf1 ! sigma value |
---|
| 713 | !!---------------------------------------------------------------------- |
---|
| 714 | ! |
---|
| 715 | IF ( rn_theta == 0 ) then ! uniform sigma |
---|
[5257] | 716 | pf1 = - ( pk1 - 0.5 ) / REAL( jpk-1 ) |
---|
[5255] | 717 | ELSE ! stretched sigma |
---|
[5257] | 718 | pf1 = ( 1. - pbb ) * ( SINH( rn_theta*(-(pk1-0.5)/REAL(jpk-1)) ) ) / SINH( rn_theta ) & |
---|
| 719 | & + pbb * ( (TANH( rn_theta*( (-(pk1-0.5)/REAL(jpk-1)) + 0.5) ) - TANH( 0.5 * rn_theta ) ) & |
---|
| 720 | & / ( 2. * TANH( 0.5 * rn_theta ) ) ) |
---|
[5255] | 721 | ENDIF |
---|
| 722 | ! |
---|
| 723 | END FUNCTION fssig1 |
---|
| 724 | |
---|
| 725 | |
---|
| 726 | FUNCTION fgamma( pk1, pzb, pzs, psmth) RESULT( p_gamma ) |
---|
| 727 | !!---------------------------------------------------------------------- |
---|
| 728 | !! *** ROUTINE fgamma *** |
---|
| 729 | !! |
---|
| 730 | !! ** Purpose : provide analytical function for the s-coordinate |
---|
| 731 | !! |
---|
| 732 | !! ** Method : the function provides the non-dimensional position of |
---|
| 733 | !! T and W (i.e. between 0 and 1) |
---|
| 734 | !! T-points at integer values (between 1 and jpk) |
---|
| 735 | !! W-points at integer values - 1/2 (between 0.5 and jpk-0.5) |
---|
| 736 | !! |
---|
| 737 | !! This method allows the maintenance of fixed surface and or |
---|
| 738 | !! bottom cell resolutions (cf. geopotential coordinates) |
---|
| 739 | !! within an analytically derived stretched S-coordinate framework. |
---|
| 740 | !! |
---|
| 741 | !! Reference : Siddorn and Furner, in prep |
---|
| 742 | !!---------------------------------------------------------------------- |
---|
[5257] | 743 | USE utils, ONLY : jpk,jk,wp |
---|
[5269] | 744 | REAL(wp), INTENT(in ) :: pk1 ! continuous "k" coordinate |
---|
| 745 | REAL(wp) :: p_gamma ! stretched coordinate |
---|
[5255] | 746 | REAL(wp), INTENT(in ) :: pzb ! Bottom box depth |
---|
| 747 | REAL(wp), INTENT(in ) :: pzs ! surface box depth |
---|
[5269] | 748 | REAL(wp), INTENT(in ) :: psmth ! Smoothing parameter |
---|
| 749 | REAL(wp) :: za1,za2,za3 ! local variables |
---|
| 750 | REAL(wp) :: zn1,zn2 ! local variables |
---|
| 751 | REAL(wp) :: za,zb,zx ! local variables |
---|
[5255] | 752 | !!---------------------------------------------------------------------- |
---|
| 753 | ! |
---|
| 754 | |
---|
| 755 | zn1 = 1./(jpk-1.) |
---|
| 756 | zn2 = 1. - zn1 |
---|
| 757 | |
---|
[5257] | 758 | za1 = (rn_alpha+2.0)*zn1**(rn_alpha+1.0)-(rn_alpha+1.0)*zn1**(rn_alpha+2.0) |
---|
| 759 | za2 = (rn_alpha+2.0)*zn2**(rn_alpha+1.0)-(rn_alpha+1.0)*zn2**(rn_alpha+2.0) |
---|
| 760 | za3 = (zn2**3.0 - za2)/( zn1**3.0 - za1) |
---|
[5255] | 761 | |
---|
| 762 | za = pzb - za3*(pzs-za1)-za2 |
---|
[5257] | 763 | za = za/( zn2-0.5*(za2+zn2**2.0) - za3*(zn1-0.5*(za1+zn1**2.0) ) ) |
---|
| 764 | zb = (pzs - za1 - za*( zn1-0.5*(za1+zn1**2.0 ) ) ) / (zn1**3.0 - za1) |
---|
| 765 | zx = 1.0-za/2.0-zb |
---|
[5255] | 766 | |
---|
[5269] | 767 | p_gamma = za*(pk1*(1.0-pk1/2.0))+zb*pk1**3.0 + & |
---|
| 768 | & zx*( (rn_alpha+2.0)*pk1**(rn_alpha+1.0)- & |
---|
| 769 | & (rn_alpha+1.0)*pk1**(rn_alpha+2.0) ) |
---|
| 770 | p_gamma = p_gamma*psmth+pk1*(1.0-psmth) |
---|
[5255] | 771 | |
---|
| 772 | ! |
---|
| 773 | END FUNCTION fgamma |
---|
| 774 | |
---|