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