Changeset 14770 for NEMO/branches/2021/dev_r13747_HPC-11_mcastril_HPDAonline_DiagGPU/tests/CANAL/MY_SRC/usrdef_istate.F90
- Timestamp:
- 2021-04-30T12:05:23+02:00 (3 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2021/dev_r13747_HPC-11_mcastril_HPDAonline_DiagGPU/tests/CANAL/MY_SRC/usrdef_istate.F90
r13472 r14770 16 16 USE dom_oce 17 17 USE phycst ! physical constants 18 USE eosbn2, ONLY: rn_a0 18 19 ! 19 20 USE in_out_manager ! I/O manager … … 26 27 PRIVATE 27 28 28 PUBLIC usr_def_istate ! called by istate.F90 29 PUBLIC usr_def_istate ! called by istate.F90 30 PUBLIC usr_def_istate_ssh ! called by sshwzv.F90 29 31 30 32 !! * Substitutions … … 37 39 CONTAINS 38 40 39 SUBROUTINE usr_def_istate( pdept, ptmask, pts, pu, pv , pssh)41 SUBROUTINE usr_def_istate( pdept, ptmask, pts, pu, pv ) 40 42 !!---------------------------------------------------------------------- 41 43 !! *** ROUTINE usr_def_istate *** … … 52 54 REAL(wp), DIMENSION(jpi,jpj,jpk) , INTENT( out) :: pu ! i-component of the velocity [m/s] 53 55 REAL(wp), DIMENSION(jpi,jpj,jpk) , INTENT( out) :: pv ! j-component of the velocity [m/s] 54 REAL(wp), DIMENSION(jpi,jpj) , INTENT( out) :: pssh ! sea-surface height55 56 ! 56 57 INTEGER :: ji, jj, jk, jl ! dummy loop indices … … 75 76 CASE(-1) ! stratif at rest 76 77 77 ! sea level:78 pssh(:,:) = 0.79 78 ! temperature: 80 79 pts(:,:,1,jp_tem) = 25. !!30._wp … … 87 86 88 87 CASE(0) ! rest 89 90 ! sea level: 91 pssh(:,:) = 0. 88 ! 92 89 ! temperature: 93 90 pts(:,:,:,jp_tem) = 10._wp … … 99 96 100 97 CASE(1) ! geostrophic zonal jet from -zjety to +zjety 101 102 ! sea level: 103 SELECT CASE( nn_fcase ) 104 CASE(0) ! f = f0 105 ! sea level: ssh = - fuy / g 106 WHERE( ABS(gphit) <= zjety ) 107 pssh(:,:) = - ff_t(:,:) * rn_uzonal * gphit(:,:) * 1.e3 / grav 108 ELSEWHERE 109 pssh(:,:) = - ff_t(:,:) * rn_uzonal * SIGN(zjety, gphit(:,:)) * 1.e3 / grav 110 END WHERE 111 CASE(1) ! f = f0 + beta*y 112 ! sea level: ssh = - u / g * ( fy + 0.5 * beta * y^2 ) 113 zbeta = 2._wp * omega * COS( rad * rn_ppgphi0 ) / ra 114 WHERE( ABS(gphit) <= zjety ) 115 pssh(:,:) = - rn_uzonal / grav * ( zf0 * gphit(:,:) * 1.e3 + 0.5 * zbeta * gphit(:,:) * gphit(:,:) * 1.e6 ) 116 ELSEWHERE 117 pssh(:,:) = - rn_uzonal / grav * ( zf0 * SIGN(zjety, gphit(:,:)) * 1.e3 & 118 & + 0.5 * zbeta * zjety * zjety * 1.e6 ) 119 END WHERE 120 END SELECT 98 ! 121 99 ! temperature: 122 100 pts(:,:,:,jp_tem) = 10._wp … … 139 117 ! 140 118 CASE(2) ! geostrophic zonal current shear 141 142 ! sea level: 143 SELECT CASE( nn_fcase ) 144 CASE(0) ! f = f0 145 ! sea level: ssh = - fuy / g 146 WHERE( ABS(gphit) <= zjety ) 147 pssh(:,:) = - ff_t(:,:) * rn_uzonal * ABS(gphit(:,:)) * 1.e3 / grav 148 ELSEWHERE 149 pssh(:,:) = - ff_t(:,:) * rn_uzonal * zjety * 1.e3 / grav 150 END WHERE 151 CASE(1) ! f = f0 + beta*y 152 ! sea level: ssh = - u / g * ( fy + 0.5 * beta * y^2 ) 153 zbeta = 2._wp * omega * COS( rad * rn_ppgphi0 ) / ra 154 WHERE( ABS(gphit) <= zjety ) 155 pssh(:,:) = - SIGN(rn_uzonal, gphit(:,:)) / grav & 156 & * ( zf0 * gphit(:,:) * 1.e3 + 0.5 * zbeta * gphit(:,:) * gphit(:,:) * 1.e6 ) 157 ELSEWHERE 158 pssh(:,:) = - SIGN(rn_uzonal, gphit(:,:)) / grav & 159 & * ( zf0 * SIGN(zjety, gphit(:,:)) * 1.e3 + 0.5 * zbeta * zjety * zjety * 1.e6 ) 160 END WHERE 161 END SELECT 119 ! 162 120 ! temperature: 163 121 pts(:,:,:,jp_tem) = 10._wp … … 176 134 ! 177 135 CASE(3) ! gaussian zonal currant 178 136 ! 179 137 ! zonal current 180 138 DO jk=1, jpkm1 … … 182 140 pu(:,:,jk) = rn_uzonal * EXP( - 0.5 * gphit(:,:)**2 / rn_lambda**2 ) 183 141 END DO 184 185 ! sea level:186 pssh(:,1) = - ff_t(:,1) / grav * pu(:,1,1) * e2t(:,1)187 DO jl=1, jpnj188 DO_2D( 0, 0, 0, 0 )189 pssh(ji,jj) = pssh(ji,jj-1) - ff_t(ji,jj) / grav * pu(ji,jj,1) * e2t(ji,jj)190 END_2D191 CALL lbc_lnk( 'usrdef_istate', pssh, 'T', 1. )192 END DO193 194 142 ! temperature: 195 143 pts(:,:,:,jp_tem) = 10._wp 196 144 ! salinity: 197 145 DO jk=1, jpkm1 198 pts(:,:,jk,jp_sal) = pssh(:,:)146 pts(:,:,jk,jp_sal) = gphit(:,:) 199 147 END DO 200 148 ! velocities: … … 202 150 ! 203 151 CASE(4) ! geostrophic zonal pulse 204 152 ! 205 153 DO_2D( 1, 1, 1, 1 ) 206 154 IF ( ABS(glamt(ji,jj)) <= zjetx ) THEN … … 210 158 ELSE 211 159 zdu = 0. 212 END 160 ENDIF 213 161 IF ( ABS(gphit(ji,jj)) <= zjety ) THEN 214 pssh(ji,jj) = - ff_t(ji,jj) * zdu * gphit(ji,jj) * 1.e3 / grav215 162 pu(ji,jj,:) = zdu 216 163 pts(ji,jj,:,jp_sal) = zdu / rn_uzonal + 1. 217 164 ELSE 218 pssh(ji,jj) = - ff_t(ji,jj) * zdu * SIGN(zjety,gphit(ji,jj)) * 1.e3 / grav219 165 pu(ji,jj,:) = 0. 220 166 pts(ji,jj,:,jp_sal) = 1. 221 END 222 END_2D 223 167 ENDIF 168 END_2D 169 ! 224 170 ! temperature: 225 171 pts(:,:,:,jp_tem) = 10._wp * ptmask(:,:,:) 226 172 pv(:,:,:) = 0. 227 228 229 173 ! 174 CASE(5) ! vortex 175 ! 230 176 zf0 = 2._wp * omega * SIN( rad * rn_ppgphi0 ) 231 zumax = rn_vtxmax * SIGN(1._wp, zf0) ! Here Anticyclonic: set zumax=-1 for cyclonic177 zumax = rn_vtxmax * SIGN(1._wp, zf0) ! Here Anticyclonic: set zumax=-1 for cyclonic 232 178 zlambda = SQRT(2._wp)*rn_lambda*1.e3 ! Horizontal scale in meters 233 179 zn2 = 3.e-3**2 … … 242 188 ! Surface pressure: P(x,y,z) = F(z) * Psurf(x,y) 243 189 zpsurf = zP0 * EXP(-(zx**2+zy**2)*zr_lambda2) - rho0 * ff_t(ji,jj) * rn_uzonal * zy 244 ! Sea level:245 pssh(ji,jj) = 0.246 DO jl=1,5247 zdt = pssh(ji,jj)248 zdzF = (1._wp - EXP(zdt-zH)) / (zH - 1._wp + EXP(-zH)) ! F'(z)249 zrho1 = rho0 * (1._wp + zn2*zdt/grav) - zdzF * zpsurf / grav ! -1/g Dz(P) = -1/g * F'(z) * Psurf(x,y)250 pssh(ji,jj) = zpsurf / (zrho1*grav) * ptmask(ji,jj,1) ! ssh = Psurf / (Rho*g)251 END DO252 190 ! temperature: 253 191 DO jk=1,jpk … … 258 196 zrho1 = zrho1 - zdzF * zpsurf / grav ! -1/g Dz(P) = -1/g * F'(z) * Psurf(x,y) 259 197 ENDIF 260 ! pts(ji,jj,jk,jp_tem) = (20._wp + (rho0-zrho1) / 0.28_wp) * ptmask(ji,jj,jk)261 pts(ji,jj,jk,jp_tem) = (10._wp + (rho0-zrho1) / 0.28_wp) * ptmask(ji,jj,jk)198 ! pts(ji,jj,jk,jp_tem) = (20._wp + (rho0-zrho1) / rn_a0) * ptmask(ji,jj,jk) 199 pts(ji,jj,jk,jp_tem) = (10._wp + (rho0-zrho1) / rn_a0) * ptmask(ji,jj,jk) 262 200 END DO 263 201 END_2D … … 299 237 ! 300 238 END SELECT 301 239 ! 240 CALL lbc_lnk( 'usrdef_istate', pts , 'T', 1. ) 241 CALL lbc_lnk( 'usrdef_istate', pu, 'U', -1., pv, 'V', -1. ) 242 243 END SUBROUTINE usr_def_istate 244 245 246 SUBROUTINE usr_def_istate_ssh( ptmask, pssh ) 247 !!---------------------------------------------------------------------- 248 !! *** ROUTINE usr_def_istate_ssh *** 249 !! 250 !! ** Purpose : Initialization of the dynamics and tracers 251 !! Here CANAL configuration 252 !! 253 !! ** Method : Set ssh 254 !!---------------------------------------------------------------------- 255 REAL(wp), DIMENSION(jpi,jpj,jpk) , INTENT(in ) :: ptmask ! t-point ocean mask [m] 256 REAL(wp), DIMENSION(jpi,jpj) , INTENT( out) :: pssh ! sea-surface height 257 ! 258 INTEGER :: ji, jj, jk, jl ! dummy loop indices 259 REAL(wp) :: zx, zy, zP0, zumax, zlambda, zr_lambda2, zn2, zf0, zH, zrho1, za, zf, zdzF 260 REAL(wp) :: zpsurf, zdyPs, zdxPs 261 REAL(wp) :: zdt, zdu, zdv 262 REAL(wp) :: zjetx, zjety, zbeta 263 REAL(wp), DIMENSION(jpi,jpj) :: zrandom 264 !!---------------------------------------------------------------------- 265 ! 266 IF(lwp) WRITE(numout,*) 267 IF(lwp) WRITE(numout,*) 'usr_def_istate_ssh : CANAL configuration, analytical definition of initial state' 268 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~ ' 269 ! 270 IF (ln_sshnoise) CALL RANDOM_NUMBER(zrandom) 271 zjetx = ABS(rn_ujetszx)/2. 272 zjety = ABS(rn_ujetszy)/2. 273 ! 274 SELECT CASE(nn_initcase) 275 CASE(0) !== rest ==! 276 ! 277 pssh(:,:) = 0. 278 ! 279 CASE(1) !== geostrophic zonal jet from -zjety to +zjety ==! 280 ! 281 SELECT CASE( nn_fcase ) 282 CASE(0) !* f = f0 : ssh = - fuy / g 283 WHERE( ABS(gphit) <= zjety ) 284 pssh(:,:) = - ff_t(:,:) * rn_uzonal * gphit(:,:) * 1.e3 / grav 285 ELSEWHERE 286 pssh(:,:) = - ff_t(:,:) * rn_uzonal * SIGN(zjety, gphit(:,:)) * 1.e3 / grav 287 END WHERE 288 CASE(1) !* f = f0 + beta*y : ssh = - u / g * ( fy + 0.5 * beta * y^2 ) 289 zbeta = 2._wp * omega * COS( rad * rn_ppgphi0 ) / ra 290 WHERE( ABS(gphit) <= zjety ) 291 pssh(:,:) = - rn_uzonal / grav * ( ff_t(:,:) * gphit(:,:) * 1.e3 + 0.5 * zbeta * gphit(:,:) * gphit(:,:) * 1.e6 ) 292 ELSEWHERE 293 pssh(:,:) = - rn_uzonal / grav * ( ff_t(:,:) * SIGN(zjety, gphit(:,:)) * 1.e3 & 294 & + 0.5 * zbeta * zjety * zjety * 1.e6 ) 295 END WHERE 296 END SELECT 297 ! 298 CASE(2) !== geostrophic zonal current shear ==! 299 ! 300 SELECT CASE( nn_fcase ) 301 CASE(0) !* f = f0 : ssh = - fuy / g 302 WHERE( ABS(gphit) <= zjety ) 303 pssh(:,:) = - ff_t(:,:) * rn_uzonal * ABS(gphit(:,:)) * 1.e3 / grav 304 ELSEWHERE 305 pssh(:,:) = - ff_t(:,:) * rn_uzonal * zjety * 1.e3 / grav 306 END WHERE 307 CASE(1) !* f = f0 + beta*y : ssh = - u / g * ( fy + 0.5 * beta * y^2 ) 308 zbeta = 2._wp * omega * COS( rad * rn_ppgphi0 ) / ra 309 WHERE( ABS(gphit) <= zjety ) 310 pssh(:,:) = - SIGN(rn_uzonal, gphit(:,:)) / grav & 311 & * ( ff_t(:,:) * gphit(:,:) * 1.e3 + 0.5 * zbeta * gphit(:,:) * gphit(:,:) * 1.e6 ) 312 ELSEWHERE 313 pssh(:,:) = - SIGN(rn_uzonal, gphit(:,:)) / grav & 314 & * ( ff_t(:,:) * SIGN(zjety, gphit(:,:)) * 1.e3 + 0.5 * zbeta * zjety * zjety * 1.e6 ) 315 END WHERE 316 END SELECT 317 ! 318 CASE(3) !== gaussian zonal currant ==! 319 ! 320 pssh(:,1) = - ff_t(:,1) / grav * e2t(:,1) * rn_uzonal * EXP( - 0.5 * gphit(:,1)**2 / rn_lambda**2 ) 321 DO jl=1, jpnj 322 DO_2D( 0, 0, 0, 0 ) 323 pssh(ji,jj) = pssh(ji,jj-1) - ff_t(ji,jj) / grav * rn_uzonal * EXP( - 0.5 * gphit(ji,jj)**2 / rn_lambda**2 ) * e2t(ji,jj) 324 END_2D 325 CALL lbc_lnk( 'usrdef_istate_ssh', pssh, 'T', 1. ) 326 END DO 327 ! 328 CASE(4) !== geostrophic zonal pulse !!st need to implement a way to separate ssh properly ==! 329 ! 330 DO_2D( 1, 1, 1, 1 ) 331 IF ( ABS(glamt(ji,jj)) <= zjetx ) THEN 332 zdu = rn_uzonal 333 ELSEIF ( ABS(glamt(ji,jj)) <= zjetx + 100. ) THEN 334 zdu = rn_uzonal * ( ( zjetx-ABS(glamt(ji,jj)) )/100. + 1. ) 335 ELSE 336 zdu = 0. 337 ENDIF 338 IF ( ABS(gphit(ji,jj)) <= zjety ) THEN 339 pssh(ji,jj) = - ff_t(ji,jj) * zdu * gphit(ji,jj) * 1.e3 / grav 340 ELSE 341 pssh(ji,jj) = - ff_t(ji,jj) * zdu * SIGN(zjety,gphit(ji,jj)) * 1.e3 / grav 342 ENDIF 343 END_2D 344 ! 345 CASE(5) !== vortex ==! 346 ! 347 zf0 = 2._wp * omega * SIN( rad * rn_ppgphi0 ) 348 zumax = rn_vtxmax * SIGN(1._wp, zf0) ! Here Anticyclonic: set zumax=-1 for cyclonic 349 zlambda = SQRT(2._wp)*rn_lambda ! Horizontal scale in meters 350 zn2 = 3.e-3**2 351 zH = 0.5_wp * 5000._wp 352 ! 353 zr_lambda2 = 1._wp / zlambda**2 354 zP0 = rho0 * zf0 * zumax * zlambda * SQRT(EXP(1._wp)/2._wp) 355 ! 356 DO_2D( 1, 1, 1, 1 ) 357 zx = glamt(ji,jj) * 1.e3 358 zy = gphit(ji,jj) * 1.e3 359 ! ! Surface pressure: P(x,y,z) = F(z) * Psurf(x,y) 360 zpsurf = zP0 * EXP(-(zx**2+zy**2)*zr_lambda2) - rho0 * ff_t(ji,jj) * rn_uzonal * zy 361 pssh(ji,jj) = 0. 362 DO jl=1,5 363 zdt = pssh(ji,jj) 364 zdzF = (1._wp - EXP(zdt-zH)) / (zH - 1._wp + EXP(-zH)) ! F'(z) 365 zrho1 = rho0 * (1._wp + zn2*zdt/grav) - zdzF * zpsurf / grav ! -1/g Dz(P) = -1/g * F'(z) * Psurf(x,y) 366 pssh(ji,jj) = zpsurf / (zrho1*grav) * ptmask(ji,jj,1) ! ssh = Psurf / (Rho*g) 367 END DO 368 END_2D 369 ! 370 END SELECT 371 ! !== add noise ==! 302 372 IF (ln_sshnoise) THEN 303 373 CALL RANDOM_SEED() 304 374 CALL RANDOM_NUMBER(zrandom) 305 375 pssh(:,:) = pssh(:,:) + ( 0.1 * zrandom(:,:) - 0.05 ) 306 END IF 307 CALL lbc_lnk( 'usrdef_istate', pssh, 'T', 1. ) 308 CALL lbc_lnk( 'usrdef_istate', pts , 'T', 1. ) 309 CALL lbc_lnk_multi( 'usrdef_istate', pu, 'U', -1., pv, 'V', -1. ) 310 311 END SUBROUTINE usr_def_istate 312 376 ENDIF 377 CALL lbc_lnk( 'usrdef_istate_ssh', pssh, 'T', 1. ) 378 ! 379 END SUBROUTINE usr_def_istate_ssh 380 313 381 !!====================================================================== 314 382 END MODULE usrdef_istate
Note: See TracChangeset
for help on using the changeset viewer.
