- Timestamp:
- 2020-09-15T12:49:18+02:00 (4 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2020/temporary_r4_trunk/tests/CANAL/MY_SRC/usrdef_istate.F90
r13466 r13469 184 184 pssh(:,1) = - ff_t(:,1) / grav * pu(:,1,1) * e2t(:,1) 185 185 DO jl=1, jpnj 186 DO jj=nldj, nlej 187 DO ji=nldi, nlei 188 pssh(ji,jj) = pssh(ji,jj-1) - ff_t(ji,jj) / grav * pu(ji,jj,1) * e2t(ji,jj) 189 END DO 190 END DO 186 DO_2D_nldj1_nldi1 187 pssh(ji,jj) = pssh(ji,jj-1) - ff_t(ji,jj) / grav * pu(ji,jj,1) * e2t(ji,jj) 188 END_2D 191 189 CALL lbc_lnk( 'usrdef_istate', pssh, 'T', 1. ) 192 190 END DO … … 203 201 CASE(4) ! geostrophic zonal pulse 204 202 205 DO jj=1, jpj 206 DO ji=1, jpi 207 IF ( ABS(glamt(ji,jj)) <= zjetx ) THEN 208 zdu = rn_uzonal 209 ELSEIF ( ABS(glamt(ji,jj)) <= zjetx + 100. ) THEN 210 zdu = rn_uzonal * ( ( zjetx-ABS(glamt(ji,jj)) )/100. + 1. ) 211 ELSE 212 zdu = 0. 213 END IF 214 IF ( ABS(gphit(ji,jj)) <= zjety ) THEN 215 pssh(ji,jj) = - ff_t(ji,jj) * zdu * gphit(ji,jj) * 1.e3 / grav 216 pu(ji,jj,:) = zdu 217 pts(ji,jj,:,jp_sal) = zdu / rn_uzonal + 1. 218 ELSE 219 pssh(ji,jj) = - ff_t(ji,jj) * zdu * SIGN(zjety,gphit(ji,jj)) * 1.e3 / grav 220 pu(ji,jj,:) = 0. 221 pts(ji,jj,:,jp_sal) = 1. 222 END IF 223 END DO 224 END DO 203 DO_2D_11_11 204 IF ( ABS(glamt(ji,jj)) <= zjetx ) THEN 205 zdu = rn_uzonal 206 ELSEIF ( ABS(glamt(ji,jj)) <= zjetx + 100. ) THEN 207 zdu = rn_uzonal * ( ( zjetx-ABS(glamt(ji,jj)) )/100. + 1. ) 208 ELSE 209 zdu = 0. 210 END IF 211 IF ( ABS(gphit(ji,jj)) <= zjety ) THEN 212 pssh(ji,jj) = - ff_t(ji,jj) * zdu * gphit(ji,jj) * 1.e3 / grav 213 pu(ji,jj,:) = zdu 214 pts(ji,jj,:,jp_sal) = zdu / rn_uzonal + 1. 215 ELSE 216 pssh(ji,jj) = - ff_t(ji,jj) * zdu * SIGN(zjety,gphit(ji,jj)) * 1.e3 / grav 217 pu(ji,jj,:) = 0. 218 pts(ji,jj,:,jp_sal) = 1. 219 END IF 220 END_2D 225 221 226 222 ! temperature: … … 240 236 zP0 = rau0 * zf0 * zumax * zlambda * SQRT(EXP(1._wp)/2._wp) 241 237 ! 242 DO jj=1, jpj 243 DO ji=1, jpi 244 zx = glamt(ji,jj) * 1.e3 245 zy = gphit(ji,jj) * 1.e3 246 ! Surface pressure: P(x,y,z) = F(z) * Psurf(x,y) 247 zpsurf = zP0 * EXP(-(zx**2+zy**2)*zr_lambda2) - rau0 * ff_t(ji,jj) * rn_uzonal * zy 248 ! Sea level: 249 pssh(ji,jj) = 0. 250 DO jl=1,5 251 zdt = pssh(ji,jj) 252 zdzF = (1._wp - EXP(zdt-zH)) / (zH - 1._wp + EXP(-zH)) ! F'(z) 253 zrho1 = rau0 * (1._wp + zn2*zdt/grav) - zdzF * zpsurf / grav ! -1/g Dz(P) = -1/g * F'(z) * Psurf(x,y) 254 pssh(ji,jj) = zpsurf / (zrho1*grav) * ptmask(ji,jj,1) ! ssh = Psurf / (Rho*g) 255 END DO 256 ! temperature: 257 DO jk=1,jpk 258 zdt = pdept(ji,jj,jk) 259 zrho1 = rau0 * (1._wp + zn2*zdt/grav) 260 IF (zdt < zH) THEN 261 zdzF = (1._wp-EXP(zdt-zH)) / (zH-1._wp + EXP(-zH)) ! F'(z) 262 zrho1 = zrho1 - zdzF * zpsurf / grav ! -1/g Dz(P) = -1/g * F'(z) * Psurf(x,y) 263 ENDIF 264 ! pts(ji,jj,jk,jp_tem) = (20._wp + (rau0-zrho1) / 0.28_wp) * ptmask(ji,jj,jk) 265 pts(ji,jj,jk,jp_tem) = (10._wp + (rau0-zrho1) / 0.28_wp) * ptmask(ji,jj,jk) 266 END DO 238 DO_2D_11_11 239 zx = glamt(ji,jj) * 1.e3 240 zy = gphit(ji,jj) * 1.e3 241 ! Surface pressure: P(x,y,z) = F(z) * Psurf(x,y) 242 zpsurf = zP0 * EXP(-(zx**2+zy**2)*zr_lambda2) - rau0 * ff_t(ji,jj) * rn_uzonal * zy 243 ! Sea level: 244 pssh(ji,jj) = 0. 245 DO jl=1,5 246 zdt = pssh(ji,jj) 247 zdzF = (1._wp - EXP(zdt-zH)) / (zH - 1._wp + EXP(-zH)) ! F'(z) 248 zrho1 = rau0 * (1._wp + zn2*zdt/grav) - zdzF * zpsurf / grav ! -1/g Dz(P) = -1/g * F'(z) * Psurf(x,y) 249 pssh(ji,jj) = zpsurf / (zrho1*grav) * ptmask(ji,jj,1) ! ssh = Psurf / (Rho*g) 267 250 END DO 268 END DO 251 ! temperature: 252 DO jk=1,jpk 253 zdt = pdept(ji,jj,jk) 254 zrho1 = rau0 * (1._wp + zn2*zdt/grav) 255 IF (zdt < zH) THEN 256 zdzF = (1._wp-EXP(zdt-zH)) / (zH-1._wp + EXP(-zH)) ! F'(z) 257 zrho1 = zrho1 - zdzF * zpsurf / grav ! -1/g Dz(P) = -1/g * F'(z) * Psurf(x,y) 258 ENDIF 259 ! pts(ji,jj,jk,jp_tem) = (20._wp + (rau0-zrho1) / 0.28_wp) * ptmask(ji,jj,jk) 260 pts(ji,jj,jk,jp_tem) = (10._wp + (rau0-zrho1) / 0.28_wp) * ptmask(ji,jj,jk) 261 END DO 262 END_2D 269 263 ! 270 264 ! salinity: … … 273 267 ! velocities: 274 268 za = 2._wp * zP0 / zlambda**2 275 DO jj = 2, jpjm1 276 DO ji = 2, jpim1 277 zx = glamu(ji,jj) * 1.e3 278 zy = gphiu(ji,jj) * 1.e3 279 DO jk=1, jpk 280 zdu = 0.5_wp * (pdept(ji,jj,jk) + pdept(ji+1,jj,jk)) 281 IF (zdu < zH) THEN 282 zf = (zH-1._wp-zdu+EXP(zdu-zH)) / (zH-1._wp+EXP(-zH)) 283 zdyPs = - za * zy * EXP(-(zx**2+zy**2)*zr_lambda2) - rau0 * ff_t(ji,jj) * rn_uzonal 284 pu(ji,jj,jk) = - zf / ( rau0 * ff_t(ji,jj) ) * zdyPs * ptmask(ji,jj,jk) * ptmask(ji+1,jj,jk) 285 ELSE 286 pu(ji,jj,jk) = 0._wp 287 ENDIF 288 END DO 269 DO_2D_00_00 270 zx = glamu(ji,jj) * 1.e3 271 zy = gphiu(ji,jj) * 1.e3 272 DO jk=1, jpk 273 zdu = 0.5_wp * (pdept(ji,jj,jk) + pdept(ji+1,jj,jk)) 274 IF (zdu < zH) THEN 275 zf = (zH-1._wp-zdu+EXP(zdu-zH)) / (zH-1._wp+EXP(-zH)) 276 zdyPs = - za * zy * EXP(-(zx**2+zy**2)*zr_lambda2) - rau0 * ff_t(ji,jj) * rn_uzonal 277 pu(ji,jj,jk) = - zf / ( rau0 * ff_t(ji,jj) ) * zdyPs * ptmask(ji,jj,jk) * ptmask(ji+1,jj,jk) 278 ELSE 279 pu(ji,jj,jk) = 0._wp 280 ENDIF 289 281 END DO 290 END DO 291 ! 292 DO jj = 2, jpjm1 293 DO ji = 2, jpim1 294 zx = glamv(ji,jj) * 1.e3 295 zy = gphiv(ji,jj) * 1.e3 296 DO jk=1, jpk 297 zdv = 0.5_wp * (pdept(ji,jj,jk) + pdept(ji,jj+1,jk)) 298 IF (zdv < zH) THEN 299 zf = (zH-1._wp-zdv+EXP(zdv-zH)) / (zH-1._wp+EXP(-zH)) 300 zdxPs = - za * zx * EXP(-(zx**2+zy**2)*zr_lambda2) 301 pv(ji,jj,jk) = zf / ( rau0 * ff_f(ji,jj) ) * zdxPs * ptmask(ji,jj,jk) * ptmask(ji,jj+1,jk) 302 ELSE 303 pv(ji,jj,jk) = 0._wp 304 ENDIF 305 END DO 282 END_2D 283 ! 284 DO_2D_00_00 285 zx = glamv(ji,jj) * 1.e3 286 zy = gphiv(ji,jj) * 1.e3 287 DO jk=1, jpk 288 zdv = 0.5_wp * (pdept(ji,jj,jk) + pdept(ji,jj+1,jk)) 289 IF (zdv < zH) THEN 290 zf = (zH-1._wp-zdv+EXP(zdv-zH)) / (zH-1._wp+EXP(-zH)) 291 zdxPs = - za * zx * EXP(-(zx**2+zy**2)*zr_lambda2) 292 pv(ji,jj,jk) = zf / ( rau0 * ff_f(ji,jj) ) * zdxPs * ptmask(ji,jj,jk) * ptmask(ji,jj+1,jk) 293 ELSE 294 pv(ji,jj,jk) = 0._wp 295 ENDIF 306 296 END DO 307 END DO297 END_2D 308 298 ! 309 299 CALL lbc_lnk_multi( 'usrdef_istate', pu, 'U', -1., pv, 'V', -1. )
Note: See TracChangeset
for help on using the changeset viewer.
