Changeset 14789 for NEMO/branches/2021/dev_r13747_HPC-11_mcastril_HPDAonline_DiagGPU/src/OCE/DYN/dynhpg.F90
- Timestamp:
- 2021-05-05T13:18:04+02:00 (3 years ago)
- Location:
- NEMO/branches/2021/dev_r13747_HPC-11_mcastril_HPDAonline_DiagGPU
- Files:
-
- 2 edited
-
. (modified) (1 prop)
-
src/OCE/DYN/dynhpg.F90 (modified) (27 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2021/dev_r13747_HPC-11_mcastril_HPDAonline_DiagGPU
- Property svn:externals
-
old new 3 3 ^/utils/build/mk@HEAD mk 4 4 ^/utils/tools@HEAD tools 5 ^/vendors/AGRIF/dev _r12970_AGRIF_CMEMSext/AGRIF5 ^/vendors/AGRIF/dev@HEAD ext/AGRIF 6 6 ^/vendors/FCM@HEAD ext/FCM 7 7 ^/vendors/IOIPSL@HEAD ext/IOIPSL 8 ^/vendors/PPR@HEAD ext/PPR 8 9 9 10 # SETTE 10 ^/utils/CI/sette@1 3559sette11 ^/utils/CI/sette@14244 sette
-
- Property svn:externals
-
NEMO/branches/2021/dev_r13747_HPC-11_mcastril_HPDAonline_DiagGPU/src/OCE/DYN/dynhpg.F90
r13295 r14789 17 17 !! ! (A. Coward) suppression of hel, wdj and rot options 18 18 !! 3.6 ! 2014-11 (P. Mathiot) hpg_isf: original code for ice shelf cavity 19 !! 4.2 ! 2020-12 (M. Bell, A. Young) hpg_djc: revised djc scheme 19 20 !!---------------------------------------------------------------------- 20 21 … … 72 73 INTEGER, PARAMETER :: np_isf = 5 ! s-coordinate similar to sco modify for isf 73 74 ! 74 INTEGER, PUBLIC :: nhpg !: type of pressure gradient scheme used ! (deduced from ln_hpg_... flags) (PUBLIC for TAM) 75 INTEGER, PUBLIC :: nhpg !: type of pressure gradient scheme used ! (deduced from ln_hpg_... flags) (PUBLIC for TAM) 76 ! 77 LOGICAL :: ln_hpg_djc_vnh, ln_hpg_djc_vnv ! flag to specify hpg_djc boundary condition type 78 REAL(wp), PUBLIC :: aco_bc_hor, bco_bc_hor, aco_bc_vrt, bco_bc_vrt !: coefficients for hpg_djc hor and vert boundary conditions 75 79 76 80 !! * Substitutions … … 150 154 !! 151 155 INTEGER :: ji, jj, jk, ikt ! dummy loop indices ISF 152 REAL(wp) :: znad153 156 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: zts_top, zrhd ! hypothesys on isf density 154 157 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zrhdtop_isf ! density at bottom of ISF … … 156 159 !! 157 160 NAMELIST/namdyn_hpg/ ln_hpg_zco, ln_hpg_zps, ln_hpg_sco, & 158 & ln_hpg_djc, ln_hpg_prj, ln_hpg_isf 161 & ln_hpg_djc, ln_hpg_prj, ln_hpg_isf, & 162 & ln_hpg_djc_vnh, ln_hpg_djc_vnv 159 163 !!---------------------------------------------------------------------- 160 164 ! … … 179 183 ENDIF 180 184 ! 181 IF( ln_hpg_djc ) & 182 & CALL ctl_stop('dyn_hpg_init : Density Jacobian: Cubic polynominal method', & 183 & ' currently disabled (bugs under investigation).' , & 184 & ' Please select either ln_hpg_sco or ln_hpg_prj instead' ) 185 ! 186 IF( .NOT.ln_linssh .AND. .NOT.(ln_hpg_sco.OR.ln_hpg_prj.OR.ln_hpg_isf) ) & 187 & CALL ctl_stop('dyn_hpg_init : non-linear free surface requires either ', & 188 & ' the standard jacobian formulation hpg_sco or ' , & 189 & ' the pressure jacobian formulation hpg_prj' ) 190 ! 185 IF( .NOT.ln_linssh .AND. (ln_hpg_zco.OR.ln_hpg_zps) ) & 186 & CALL ctl_stop( 'dyn_hpg_init : non-linear free surface incompatible with hpg_zco or hpg_zps' ) 187 ! 188 IF( (.NOT.ln_hpg_isf .AND. ln_isfcav) .OR. (ln_hpg_isf .AND. .NOT.ln_isfcav) ) & 189 & CALL ctl_stop( 'dyn_hpg_init : ln_hpg_isf=T requires ln_isfcav=T and vice versa' ) 190 ! 191 #if defined key_qco 191 192 IF( ln_hpg_isf ) THEN 192 IF( .NOT. ln_isfcav ) CALL ctl_stop( ' hpg_isf not available if ln_isfcav = false ' ) 193 ELSE 194 IF( ln_isfcav ) CALL ctl_stop( 'Only hpg_isf has been corrected to work with ice shelf cavity.' ) 193 CALL ctl_stop( 'dyn_hpg_init : key_qco and ln_hpg_isf not yet compatible' ) 195 194 ENDIF 195 #endif 196 196 ! 197 197 ! ! Set nhpg from ln_hpg_... flags & consistency check … … 220 220 ENDIF 221 221 ! 222 IF ( ln_hpg_djc ) THEN 223 IF (ln_hpg_djc_vnh) THEN ! Von Neumann boundary condition 224 IF(lwp) WRITE(numout,*) ' horizontal bc: von Neumann ' 225 aco_bc_hor = 6.0_wp/5.0_wp 226 bco_bc_hor = 7.0_wp/15.0_wp 227 ELSE ! Linear extrapolation 228 IF(lwp) WRITE(numout,*) ' horizontal bc: linear extrapolation' 229 aco_bc_hor = 3.0_wp/2.0_wp 230 bco_bc_hor = 1.0_wp/2.0_wp 231 END IF 232 IF (ln_hpg_djc_vnv) THEN ! Von Neumann boundary condition 233 IF(lwp) WRITE(numout,*) ' vertical bc: von Neumann ' 234 aco_bc_vrt = 6.0_wp/5.0_wp 235 bco_bc_vrt = 7.0_wp/15.0_wp 236 ELSE ! Linear extrapolation 237 IF(lwp) WRITE(numout,*) ' vertical bc: linear extrapolation' 238 aco_bc_vrt = 3.0_wp/2.0_wp 239 bco_bc_vrt = 1.0_wp/2.0_wp 240 END IF 241 END IF 242 ! 222 243 END SUBROUTINE dyn_hpg_init 223 244 … … 245 266 INTEGER :: ji, jj, jk ! dummy loop indices 246 267 REAL(wp) :: zcoef0, zcoef1 ! temporary scalars 247 REAL(wp), DIMENSION(jpi,jpj ,jpk) :: zhpi, zhpj268 REAL(wp), DIMENSION(jpi,jpj) :: zhpi, zhpj 248 269 !!---------------------------------------------------------------------- 249 270 ! … … 253 274 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ z-coordinate case ' 254 275 ENDIF 255 256 zcoef0 = - grav * 0.5_wp ! Local constant initialization 257 258 ! Surface value 259 DO_2D( 0, 0, 0, 0 ) 276 ! 277 zcoef0 = - grav * 0.5_wp ! Local constant initialization 278 ! 279 DO_2D( 0, 0, 0, 0 ) ! Surface value 260 280 zcoef1 = zcoef0 * e3w(ji,jj,1,Kmm) 261 ! hydrostatic pressure gradient262 zhpi(ji,jj ,1) = zcoef1 * ( rhd(ji+1,jj,1) - rhd(ji,jj,1) ) * r1_e1u(ji,jj)263 zhpj(ji,jj ,1) = zcoef1 * ( rhd(ji,jj+1,1) - rhd(ji,jj,1) ) * r1_e2v(ji,jj)264 ! add to the general momentum trend265 puu(ji,jj,1,Krhs) = puu(ji,jj,1,Krhs) + zhpi(ji,jj ,1)266 pvv(ji,jj,1,Krhs) = pvv(ji,jj,1,Krhs) + zhpj(ji,jj ,1)281 ! ! hydrostatic pressure gradient 282 zhpi(ji,jj) = zcoef1 * ( rhd(ji+1,jj,1) - rhd(ji,jj,1) ) * r1_e1u(ji,jj) 283 zhpj(ji,jj) = zcoef1 * ( rhd(ji,jj+1,1) - rhd(ji,jj,1) ) * r1_e2v(ji,jj) 284 ! ! add to the general momentum trend 285 puu(ji,jj,1,Krhs) = puu(ji,jj,1,Krhs) + zhpi(ji,jj) 286 pvv(ji,jj,1,Krhs) = pvv(ji,jj,1,Krhs) + zhpj(ji,jj) 267 287 END_2D 268 269 ! 270 ! interior value (2=<jk=<jpkm1) 271 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 288 ! 289 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) ! interior value (2=<jk=<jpkm1) 272 290 zcoef1 = zcoef0 * e3w(ji,jj,jk,Kmm) 273 ! hydrostatic pressure gradient 274 zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) & 275 & + zcoef1 * ( ( rhd(ji+1,jj,jk)+rhd(ji+1,jj,jk-1) ) & 276 & - ( rhd(ji ,jj,jk)+rhd(ji ,jj,jk-1) ) ) * r1_e1u(ji,jj) 277 278 zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1) & 279 & + zcoef1 * ( ( rhd(ji,jj+1,jk)+rhd(ji,jj+1,jk-1) ) & 280 & - ( rhd(ji,jj, jk)+rhd(ji,jj ,jk-1) ) ) * r1_e2v(ji,jj) 281 ! add to the general momentum trend 282 puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) + zhpi(ji,jj,jk) 283 pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) + zhpj(ji,jj,jk) 291 ! ! hydrostatic pressure gradient 292 zhpi(ji,jj) = zhpi(ji,jj) + zcoef1 * ( ( rhd(ji+1,jj,jk)+rhd(ji+1,jj,jk-1) ) & 293 & - ( rhd(ji ,jj,jk)+rhd(ji ,jj,jk-1) ) ) * r1_e1u(ji,jj) 294 295 zhpj(ji,jj) = zhpj(ji,jj) + zcoef1 * ( ( rhd(ji,jj+1,jk)+rhd(ji,jj+1,jk-1) ) & 296 & - ( rhd(ji,jj, jk)+rhd(ji,jj ,jk-1) ) ) * r1_e2v(ji,jj) 297 ! ! add to the general momentum trend 298 puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) + zhpi(ji,jj) 299 pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) + zhpj(ji,jj) 284 300 END_3D 285 301 ! … … 302 318 INTEGER :: iku, ikv ! temporary integers 303 319 REAL(wp) :: zcoef0, zcoef1, zcoef2, zcoef3 ! temporary scalars 304 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zhpi, zhpj 305 REAL(wp), DIMENSION(jpi,jpj) :: zgtsu, zgtsv, zgru, zgrv 320 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zhpi, zhpj 321 REAL(wp), DIMENSION(jpi,jpj,jpts) :: zgtsu, zgtsv 322 REAL(wp), DIMENSION(jpi,jpj) :: zgru, zgrv 306 323 !!---------------------------------------------------------------------- 307 324 ! … … 390 407 REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) :: puu, pvv ! ocean velocities and RHS of momentum equation 391 408 !! 392 INTEGER :: ji, jj, jk, jii, jjj 393 REAL(wp) :: zcoef0, zuap, zvap, z nad, ztmp ! temporaryscalars394 LOGICAL :: ll_tmp1, ll_tmp2 409 INTEGER :: ji, jj, jk, jii, jjj ! dummy loop indices 410 REAL(wp) :: zcoef0, zuap, zvap, ztmp ! local scalars 411 LOGICAL :: ll_tmp1, ll_tmp2 ! local logical variables 395 412 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zhpi, zhpj 396 413 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: zcpx, zcpy !W/D pressure filter … … 402 419 IF(lwp) WRITE(numout,*) 403 420 IF(lwp) WRITE(numout,*) 'dyn:hpg_sco : hydrostatic pressure gradient trend' 404 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ s-coordinate case, O PAoriginal scheme used'421 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ s-coordinate case, OCE original scheme used' 405 422 ENDIF 406 423 ! 407 424 zcoef0 = - grav * 0.5_wp 408 IF ( ln_linssh ) THEN ; znad = 0._wp ! Fixed volume: density anomaly409 ELSE ; znad = 1._wp ! Variable volume: density410 ENDIF411 425 ! 412 426 IF( ln_wd_il ) THEN … … 448 462 END IF 449 463 END_2D 450 CALL lbc_lnk _multi( 'dynhpg', zcpx, 'U', 1.0_wp, zcpy, 'V', 1.0_wp )464 CALL lbc_lnk( 'dynhpg', zcpx, 'U', 1.0_wp, zcpy, 'V', 1.0_wp ) 451 465 END IF 452 453 ! Surface value 454 DO_2D( 0, 0, 0, 0 ) 455 ! hydrostatic pressure gradient along s-surfaces 456 zhpi(ji,jj,1) = & 457 & zcoef0 * ( e3w(ji+1,jj ,1,Kmm) * ( znad + rhd(ji+1,jj ,1) ) & 458 & - e3w(ji ,jj ,1,Kmm) * ( znad + rhd(ji ,jj ,1) ) ) & 459 & * r1_e1u(ji,jj) 460 zhpj(ji,jj,1) = & 461 & zcoef0 * ( e3w(ji ,jj+1,1,Kmm) * ( znad + rhd(ji ,jj+1,1) ) & 462 & - e3w(ji ,jj ,1,Kmm) * ( znad + rhd(ji ,jj ,1) ) ) & 463 & * r1_e2v(ji,jj) 464 ! s-coordinate pressure gradient correction 465 zuap = -zcoef0 * ( rhd (ji+1,jj,1) + rhd (ji,jj,1) + 2._wp * znad ) & 466 ! 467 DO_2D( 0, 0, 0, 0 ) ! Surface value 468 ! ! hydrostatic pressure gradient along s-surfaces 469 zhpi(ji,jj,1) = zcoef0 * r1_e1u(ji,jj) & 470 & * ( e3w(ji+1,jj ,1,Kmm) * rhd(ji+1,jj ,1) & 471 & - e3w(ji ,jj ,1,Kmm) * rhd(ji ,jj ,1) ) 472 zhpj(ji,jj,1) = zcoef0 * r1_e2v(ji,jj) & 473 & * ( e3w(ji ,jj+1,1,Kmm) * rhd(ji ,jj+1,1) & 474 & - e3w(ji ,jj ,1,Kmm) * rhd(ji ,jj ,1) ) 475 ! ! s-coordinate pressure gradient correction 476 zuap = -zcoef0 * ( rhd (ji+1,jj,1) + rhd (ji,jj,1) ) & 466 477 & * ( gde3w(ji+1,jj,1) - gde3w(ji,jj,1) ) * r1_e1u(ji,jj) 467 zvap = -zcoef0 * ( rhd (ji,jj+1,1) + rhd (ji,jj,1) + 2._wp * znad) &478 zvap = -zcoef0 * ( rhd (ji,jj+1,1) + rhd (ji,jj,1) ) & 468 479 & * ( gde3w(ji,jj+1,1) - gde3w(ji,jj,1) ) * r1_e2v(ji,jj) 469 480 ! … … 474 485 zvap = zvap * zcpy(ji,jj) 475 486 ENDIF 476 ! 477 ! add to the general momentum trend 487 ! ! add to the general momentum trend 478 488 puu(ji,jj,1,Krhs) = puu(ji,jj,1,Krhs) + zhpi(ji,jj,1) + zuap 479 489 pvv(ji,jj,1,Krhs) = pvv(ji,jj,1,Krhs) + zhpj(ji,jj,1) + zvap 480 490 END_2D 481 482 ! interior value (2=<jk=<jpkm1) 483 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 484 ! hydrostatic pressure gradient along s-surfaces 485 zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) + zcoef0 * r1_e1u(ji,jj) & 486 & * ( e3w(ji+1,jj,jk,Kmm) * ( rhd(ji+1,jj,jk) + rhd(ji+1,jj,jk-1) + 2*znad ) & 487 & - e3w(ji ,jj,jk,Kmm) * ( rhd(ji ,jj,jk) + rhd(ji ,jj,jk-1) + 2*znad ) ) 488 zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1) + zcoef0 * r1_e2v(ji,jj) & 489 & * ( e3w(ji,jj+1,jk,Kmm) * ( rhd(ji,jj+1,jk) + rhd(ji,jj+1,jk-1) + 2*znad ) & 490 & - e3w(ji,jj ,jk,Kmm) * ( rhd(ji,jj, jk) + rhd(ji,jj ,jk-1) + 2*znad ) ) 491 ! s-coordinate pressure gradient correction 492 zuap = -zcoef0 * ( rhd (ji+1,jj ,jk) + rhd (ji,jj,jk) + 2._wp * znad ) & 491 ! 492 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) ! interior value (2=<jk=<jpkm1) 493 ! ! hydrostatic pressure gradient along s-surfaces 494 zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) + zcoef0 * r1_e1u(ji,jj) & 495 & * ( e3w(ji+1,jj,jk,Kmm) * ( rhd(ji+1,jj,jk) + rhd(ji+1,jj,jk-1) ) & 496 & - e3w(ji ,jj,jk,Kmm) * ( rhd(ji ,jj,jk) + rhd(ji ,jj,jk-1) ) ) 497 zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1) + zcoef0 * r1_e2v(ji,jj) & 498 & * ( e3w(ji,jj+1,jk,Kmm) * ( rhd(ji,jj+1,jk) + rhd(ji,jj+1,jk-1) ) & 499 & - e3w(ji,jj ,jk,Kmm) * ( rhd(ji,jj, jk) + rhd(ji,jj ,jk-1) ) ) 500 ! ! s-coordinate pressure gradient correction 501 zuap = -zcoef0 * ( rhd (ji+1,jj ,jk) + rhd (ji,jj,jk) ) & 493 502 & * ( gde3w(ji+1,jj ,jk) - gde3w(ji,jj,jk) ) * r1_e1u(ji,jj) 494 zvap = -zcoef0 * ( rhd (ji ,jj+1,jk) + rhd (ji,jj,jk) + 2._wp * znad )&503 zvap = -zcoef0 * ( rhd (ji ,jj+1,jk) + rhd (ji,jj,jk) ) & 495 504 & * ( gde3w(ji ,jj+1,jk) - gde3w(ji,jj,jk) ) * r1_e2v(ji,jj) 496 505 ! … … 535 544 REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) :: puu, pvv ! ocean velocities and RHS of momentum equation 536 545 !! 537 INTEGER :: ji, jj, jk, ikt, iktp1i, iktp1j ! dummy loop indices 538 REAL(wp) :: zcoef0, zuap, zvap, znad ! temporary scalars 546 INTEGER :: ji, jj, jk ! dummy loop indices 547 INTEGER :: ikt , ikti1, iktj1 ! local integer 548 REAL(wp) :: ze3w, ze3wi1, ze3wj1 ! local scalars 549 REAL(wp) :: zcoef0, zuap, zvap ! - - 539 550 REAL(wp), DIMENSION(jpi,jpj,jpk ) :: zhpi, zhpj 540 551 REAL(wp), DIMENSION(jpi,jpj,jpts) :: zts_top … … 543 554 ! 544 555 zcoef0 = - grav * 0.5_wp ! Local constant initialization 545 !546 znad=1._wp ! To use density and not density anomaly547 556 ! 548 557 ! ! iniitialised to 0. zhpi zhpi … … 560 569 CALL eos( zts_top, risfdep, zrhdtop_oce ) 561 570 562 !================================================================================== 563 !===== Compute surface value ===================================================== 564 !================================================================================== 571 ! !===========================! 572 ! !===== surface value =====! 573 ! !===========================! 565 574 DO_2D( 0, 0, 0, 0 ) 566 ikt = mikt(ji,jj) 567 iktp1i = mikt(ji+1,jj) 568 iktp1j = mikt(ji,jj+1) 569 ! hydrostatic pressure gradient along s-surfaces and ice shelf pressure 570 ! we assume ISF is in isostatic equilibrium 571 zhpi(ji,jj,1) = zcoef0 / e1u(ji,jj) * ( 0.5_wp * e3w(ji+1,jj,iktp1i,Kmm) & 572 & * ( 2._wp * znad + rhd(ji+1,jj,iktp1i) + zrhdtop_oce(ji+1,jj) ) & 573 & - 0.5_wp * e3w(ji,jj,ikt,Kmm) & 574 & * ( 2._wp * znad + rhd(ji,jj,ikt) + zrhdtop_oce(ji,jj) ) & 575 & + ( risfload(ji+1,jj) - risfload(ji,jj)) ) 576 zhpj(ji,jj,1) = zcoef0 / e2v(ji,jj) * ( 0.5_wp * e3w(ji,jj+1,iktp1j,Kmm) & 577 & * ( 2._wp * znad + rhd(ji,jj+1,iktp1j) + zrhdtop_oce(ji,jj+1) ) & 578 & - 0.5_wp * e3w(ji,jj,ikt,Kmm) & 579 & * ( 2._wp * znad + rhd(ji,jj,ikt) + zrhdtop_oce(ji,jj) ) & 580 & + ( risfload(ji,jj+1) - risfload(ji,jj)) ) 581 ! s-coordinate pressure gradient correction (=0 if z coordinate) 582 zuap = -zcoef0 * ( rhd (ji+1,jj,1) + rhd (ji,jj,1) + 2._wp * znad ) & 575 ikt = mikt(ji ,jj ) ; ze3w = e3w(ji ,jj ,ikt ,Kmm) 576 ikti1 = mikt(ji+1,jj ) ; ze3wi1 = e3w(ji+1,jj ,ikti1,Kmm) 577 iktj1 = mikt(ji ,jj+1) ; ze3wj1 = e3w(ji ,jj+1,iktj1,Kmm) 578 ! ! hydrostatic pressure gradient along s-surfaces and ice shelf pressure 579 ! ! we assume ISF is in isostatic equilibrium 580 zhpi(ji,jj,1) = zcoef0 * r1_e1u(ji,jj) * ( risfload(ji+1,jj) - risfload(ji,jj) & 581 & + 0.5_wp * ( ze3wi1 * ( rhd(ji+1,jj,ikti1) + zrhdtop_oce(ji+1,jj) ) & 582 & - ze3w * ( rhd(ji ,jj,ikt ) + zrhdtop_oce(ji ,jj) ) ) ) 583 zhpj(ji,jj,1) = zcoef0 * r1_e2v(ji,jj) * ( risfload(ji,jj+1) - risfload(ji,jj) & 584 & + 0.5_wp * ( ze3wj1 * ( rhd(ji,jj+1,iktj1) + zrhdtop_oce(ji,jj+1) ) & 585 & - ze3w * ( rhd(ji,jj ,ikt ) + zrhdtop_oce(ji,jj ) ) ) ) 586 ! ! s-coordinate pressure gradient correction (=0 if z coordinate) 587 zuap = -zcoef0 * ( rhd (ji+1,jj,1) + rhd (ji,jj,1) ) & 583 588 & * ( gde3w(ji+1,jj,1) - gde3w(ji,jj,1) ) * r1_e1u(ji,jj) 584 zvap = -zcoef0 * ( rhd (ji,jj+1,1) + rhd (ji,jj,1) + 2._wp * znad) &589 zvap = -zcoef0 * ( rhd (ji,jj+1,1) + rhd (ji,jj,1) ) & 585 590 & * ( gde3w(ji,jj+1,1) - gde3w(ji,jj,1) ) * r1_e2v(ji,jj) 586 ! add to the general momentum trend591 ! ! add to the general momentum trend 587 592 puu(ji,jj,1,Krhs) = puu(ji,jj,1,Krhs) + (zhpi(ji,jj,1) + zuap) * umask(ji,jj,1) 588 593 pvv(ji,jj,1,Krhs) = pvv(ji,jj,1,Krhs) + (zhpj(ji,jj,1) + zvap) * vmask(ji,jj,1) 589 594 END_2D 590 !==================================================================================591 !===== Compute interior value ===================================================== 592 !================================================================================== 593 ! interior value (2=<jk=<jpkm1)595 ! 596 ! !=============================! 597 ! !===== interior values =====! 598 ! !=============================! 594 599 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 595 ! hydrostatic pressure gradient along s-surfaces 600 ze3w = e3w(ji ,jj ,jk,Kmm) 601 ze3wi1 = e3w(ji+1,jj ,jk,Kmm) 602 ze3wj1 = e3w(ji ,jj+1,jk,Kmm) 603 ! ! hydrostatic pressure gradient along s-surfaces 596 604 zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) + zcoef0 / e1u(ji,jj) & 597 & * ( e3w(ji+1,jj,jk,Kmm) & 598 & * ( rhd(ji+1,jj,jk) + rhd(ji+1,jj,jk-1) + 2*znad ) * wmask(ji+1,jj,jk) & 599 & - e3w(ji ,jj,jk,Kmm) & 600 & * ( rhd(ji ,jj,jk) + rhd(ji ,jj,jk-1) + 2*znad ) * wmask(ji ,jj,jk) ) 605 & * ( ze3wi1 * ( rhd(ji+1,jj,jk) + rhd(ji+1,jj,jk-1) ) * wmask(ji+1,jj,jk) & 606 & - ze3w * ( rhd(ji ,jj,jk) + rhd(ji ,jj,jk-1) ) * wmask(ji ,jj,jk) ) 601 607 zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1) + zcoef0 / e2v(ji,jj) & 602 & * ( e3w(ji,jj+1,jk,Kmm) & 603 & * ( rhd(ji,jj+1,jk) + rhd(ji,jj+1,jk-1) + 2*znad ) * wmask(ji,jj+1,jk) & 604 & - e3w(ji,jj ,jk,Kmm) & 605 & * ( rhd(ji,jj, jk) + rhd(ji,jj ,jk-1) + 2*znad ) * wmask(ji,jj ,jk) ) 606 ! s-coordinate pressure gradient correction 607 zuap = -zcoef0 * ( rhd (ji+1,jj ,jk) + rhd (ji,jj,jk) + 2._wp * znad ) & 608 & * ( ze3wj1 * ( rhd(ji,jj+1,jk) + rhd(ji,jj+1,jk-1) ) * wmask(ji,jj+1,jk) & 609 & - ze3w * ( rhd(ji,jj, jk) + rhd(ji,jj ,jk-1) ) * wmask(ji,jj ,jk) ) 610 ! ! s-coordinate pressure gradient correction 611 zuap = -zcoef0 * ( rhd (ji+1,jj ,jk) + rhd (ji,jj,jk) ) & 608 612 & * ( gde3w(ji+1,jj ,jk) - gde3w(ji,jj,jk) ) / e1u(ji,jj) 609 zvap = -zcoef0 * ( rhd (ji ,jj+1,jk) + rhd (ji,jj,jk) + 2._wp * znad) &613 zvap = -zcoef0 * ( rhd (ji ,jj+1,jk) + rhd (ji,jj,jk) ) & 610 614 & * ( gde3w(ji ,jj+1,jk) - gde3w(ji,jj,jk) ) / e2v(ji,jj) 611 ! add to the general momentum trend615 ! ! add to the general momentum trend 612 616 puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) + (zhpi(ji,jj,jk) + zuap) * umask(ji,jj,jk) 613 617 pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) + (zhpj(ji,jj,jk) + zvap) * vmask(ji,jj,jk) … … 630 634 !! 631 635 INTEGER :: ji, jj, jk ! dummy loop indices 636 INTEGER :: iktb, iktt ! jk indices at tracer points for top and bottom points 632 637 REAL(wp) :: zcoef0, zep, cffw ! temporary scalars 633 REAL(wp) :: z1_10, cffu, cffx ! " " 634 REAL(wp) :: z1_12, cffv, cffy ! " " 638 REAL(wp) :: z_grav_10, z1_12 639 REAL(wp) :: cffu, cffx ! " " 640 REAL(wp) :: cffv, cffy ! " " 635 641 LOGICAL :: ll_tmp1, ll_tmp2 ! local logical variables 636 642 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zhpi, zhpj 637 REAL(wp), DIMENSION(jpi,jpj,jpk) :: dzx, dzy, dzz, dzu, dzv, dzw 638 REAL(wp), DIMENSION(jpi,jpj,jpk) :: drhox, drhoy, drhoz, drhou, drhov, drhow 639 REAL(wp), DIMENSION(jpi,jpj,jpk) :: rho_i, rho_j, rho_k 643 644 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdzx, zdzy, zdzz ! Primitive grid differences ('delta_xyz') 645 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdz_i, zdz_j, zdz_k ! Harmonic average of primitive grid differences ('d_xyz') 646 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdrhox, zdrhoy, zdrhoz ! Primitive rho differences ('delta_rho') 647 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdrho_i, zdrho_j, zdrho_k ! Harmonic average of primitive rho differences ('d_rho') 648 REAL(wp), DIMENSION(jpi,jpj,jpk) :: z_rho_i, z_rho_j, z_rho_k ! Face intergrals 649 REAL(wp), DIMENSION(jpi,jpj) :: zz_dz_i, zz_dz_j, zz_drho_i, zz_drho_j ! temporary arrays 640 650 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: zcpx, zcpy !W/D pressure filter 641 651 !!---------------------------------------------------------------------- … … 679 689 END IF 680 690 END_2D 681 CALL lbc_lnk _multi( 'dynhpg', zcpx, 'U', 1.0_wp, zcpy, 'V', 1.0_wp )691 CALL lbc_lnk( 'dynhpg', zcpx, 'U', 1.0_wp, zcpy, 'V', 1.0_wp ) 682 692 END IF 683 693 … … 690 700 ! Local constant initialization 691 701 zcoef0 = - grav * 0.5_wp 692 z 1_10 = 1._wp/ 10._wp693 z1_12 = 1. _wp / 12._wp702 z_grav_10 = grav / 10._wp 703 z1_12 = 1.0_wp / 12._wp 694 704 695 705 !---------------------------------------------------------------------------------------- 696 ! compute and store in provisional arrays elementary vertical and horizontal differences706 ! 1. compute and store elementary vertical differences in provisional arrays 697 707 !---------------------------------------------------------------------------------------- 698 708 699 !!bug gm Not a true bug, but... dzz=e3w for dzx, dzy verify what it is really 700 701 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 702 drhoz(ji,jj,jk) = rhd (ji ,jj ,jk) - rhd (ji,jj,jk-1) 703 dzz (ji,jj,jk) = gde3w(ji ,jj ,jk) - gde3w(ji,jj,jk-1) 704 drhox(ji,jj,jk) = rhd (ji+1,jj ,jk) - rhd (ji,jj,jk ) 705 dzx (ji,jj,jk) = gde3w(ji+1,jj ,jk) - gde3w(ji,jj,jk ) 706 drhoy(ji,jj,jk) = rhd (ji ,jj+1,jk) - rhd (ji,jj,jk ) 707 dzy (ji,jj,jk) = gde3w(ji ,jj+1,jk) - gde3w(ji,jj,jk ) 709 !!bug gm Not a true bug, but... zdzz=e3w for zdzx, zdzy verify what it is really 710 711 DO_3D( 1, 1, 1, 1, 2, jpkm1 ) 712 zdrhoz(ji,jj,jk) = rhd (ji ,jj ,jk) - rhd (ji,jj,jk-1) 713 zdzz (ji,jj,jk) = - gde3w(ji ,jj ,jk) + gde3w(ji,jj,jk-1) 708 714 END_3D 709 715 710 716 !------------------------------------------------------------------------- 711 ! compute harmonic averages using eq. 5.18717 ! 2. compute harmonic averages for vertical differences using eq. 5.18 712 718 !------------------------------------------------------------------------- 713 719 zep = 1.e-15 714 720 715 !!bug gm drhoz not defined at level 1 and used (jk-1 with jk=2) 716 !!bug gm idem for drhox, drhoy et ji=jpi and jj=jpj 717 718 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 719 cffw = 2._wp * drhoz(ji ,jj ,jk) * drhoz(ji,jj,jk-1) 720 721 cffu = 2._wp * drhox(ji+1,jj ,jk) * drhox(ji,jj,jk ) 722 cffx = 2._wp * dzx (ji+1,jj ,jk) * dzx (ji,jj,jk ) 723 724 cffv = 2._wp * drhoy(ji ,jj+1,jk) * drhoy(ji,jj,jk ) 725 cffy = 2._wp * dzy (ji ,jj+1,jk) * dzy (ji,jj,jk ) 726 721 !! mb zdrho_k, zdz_k, zdrho_i, zdz_i, zdrho_j, zdz_j re-centred about the point (ji,jj,jk) 722 zdrho_k(:,:,:) = 0._wp 723 zdz_k (:,:,:) = 0._wp 724 725 DO_3D( 1, 1, 1, 1, 2, jpk-2 ) 726 cffw = 2._wp * zdrhoz(ji ,jj ,jk) * zdrhoz(ji,jj,jk+1) 727 727 IF( cffw > zep) THEN 728 drhow(ji,jj,jk) = 2._wp * drhoz(ji,jj,jk) * drhoz(ji,jj,jk-1) & 729 & / ( drhoz(ji,jj,jk) + drhoz(ji,jj,jk-1) ) 728 zdrho_k(ji,jj,jk) = cffw / ( zdrhoz(ji,jj,jk) + zdrhoz(ji,jj,jk+1) ) 729 ENDIF 730 zdz_k(ji,jj,jk) = 2._wp * zdzz(ji,jj,jk) * zdzz(ji,jj,jk+1) & 731 & / ( zdzz(ji,jj,jk) + zdzz(ji,jj,jk+1) ) 732 END_3D 733 734 !---------------------------------------------------------------------------------- 735 ! 3. apply boundary conditions at top and bottom using 5.36-5.37 736 !---------------------------------------------------------------------------------- 737 738 ! mb for sea-ice shelves we will need to re-write this upper boundary condition in the same form as the lower boundary condition 739 zdrho_k(:,:,1) = aco_bc_vrt * ( rhd (:,:,2) - rhd (:,:,1) ) - bco_bc_vrt * zdrho_k(:,:,2) 740 zdz_k (:,:,1) = aco_bc_vrt * (-gde3w(:,:,2) + gde3w(:,:,1) ) - bco_bc_vrt * zdz_k (:,:,2) 741 742 DO_2D( 1, 1, 1, 1 ) 743 IF ( mbkt(ji,jj)>1 ) THEN 744 iktb = mbkt(ji,jj) 745 zdrho_k(ji,jj,iktb) = aco_bc_vrt * ( rhd(ji,jj,iktb) - rhd(ji,jj,iktb-1) ) - bco_bc_vrt * zdrho_k(ji,jj,iktb-1) 746 zdz_k (ji,jj,iktb) = aco_bc_vrt * (-gde3w(ji,jj,iktb) + gde3w(ji,jj,iktb-1) ) - bco_bc_vrt * zdz_k (ji,jj,iktb-1) 747 END IF 748 END_2D 749 750 !-------------------------------------------------------------- 751 ! 4. Compute side face integrals 752 !------------------------------------------------------------- 753 754 !! ssh replaces e3w_n ; gde3w is a depth; the formulae involve heights 755 !! rho_k stores grav * FX / rho_0 756 757 !-------------------------------------------------------------- 758 ! 4. a) Upper half of top-most grid box, compute and store 759 !------------------------------------------------------------- 760 ! *** AY note: ssh(ji,jj,Kmm) + gde3w(ji,jj,1) = e3w(ji,jj,1,Kmm) 761 DO_2D( 0, 1, 0, 1) 762 z_rho_k(ji,jj,1) = grav * ( ssh(ji,jj,Kmm) + gde3w(ji,jj,1) ) & 763 & * ( rhd(ji,jj,1) & 764 & + 0.5_wp * ( rhd (ji,jj,2) - rhd (ji,jj,1) ) & 765 & * ( ssh (ji,jj,Kmm) + gde3w(ji,jj,1) ) & 766 & / ( - gde3w(ji,jj,2) + gde3w(ji,jj,1) ) ) 767 END_2D 768 769 !-------------------------------------------------------------- 770 ! 4. b) Interior faces, compute and store 771 !------------------------------------------------------------- 772 773 DO_3D( 0, 1, 0, 1, 2, jpkm1 ) 774 z_rho_k(ji,jj,jk) = zcoef0 * ( rhd (ji,jj,jk) + rhd (ji,jj,jk-1) ) & 775 & * ( - gde3w(ji,jj,jk) + gde3w(ji,jj,jk-1) ) & 776 & + z_grav_10 * ( & 777 & ( zdrho_k (ji,jj,jk) - zdrho_k (ji,jj,jk-1) ) & 778 & * ( - gde3w(ji,jj,jk) + gde3w(ji,jj,jk-1) - z1_12 * ( zdz_k (ji,jj,jk) + zdz_k (ji,jj,jk-1) ) ) & 779 & - ( zdz_k (ji,jj,jk) - zdz_k (ji,jj,jk-1) ) & 780 & * ( rhd (ji,jj,jk) - rhd (ji,jj,jk-1) - z1_12 * ( zdrho_k(ji,jj,jk) + zdrho_k(ji,jj,jk-1) ) ) & 781 & ) 782 END_3D 783 784 !---------------------------------------------------------------------------------------- 785 ! 5. compute and store elementary horizontal differences in provisional arrays 786 !---------------------------------------------------------------------------------------- 787 788 DO_3D( 1, 0, 1, 0, 1, jpkm1 ) 789 zdrhox(ji,jj,jk) = rhd (ji+1,jj ,jk) - rhd (ji,jj,jk ) 790 zdzx (ji,jj,jk) = - gde3w(ji+1,jj ,jk) + gde3w(ji,jj,jk ) 791 zdrhoy(ji,jj,jk) = rhd (ji ,jj+1,jk) - rhd (ji,jj,jk ) 792 zdzy (ji,jj,jk) = - gde3w(ji ,jj+1,jk) + gde3w(ji,jj,jk ) 793 END_3D 794 795 CALL lbc_lnk( 'dynhpg', zdrhox, 'U', 1., zdzx, 'U', 1., zdrhoy, 'V', 1., zdzy, 'V', 1. ) 796 797 !------------------------------------------------------------------------- 798 ! 6. compute harmonic averages using eq. 5.18 799 !------------------------------------------------------------------------- 800 801 DO_3D( 0, 1, 0, 1, 1, jpkm1 ) 802 cffu = 2._wp * zdrhox(ji-1,jj ,jk) * zdrhox(ji,jj,jk ) 803 IF( cffu > zep ) THEN 804 zdrho_i(ji,jj,jk) = cffu / ( zdrhox(ji-1,jj,jk) + zdrhox(ji,jj,jk) ) 730 805 ELSE 731 drhow(ji,jj,jk) = 0._wp 732 ENDIF 733 734 dzw(ji,jj,jk) = 2._wp * dzz(ji,jj,jk) * dzz(ji,jj,jk-1) & 735 & / ( dzz(ji,jj,jk) + dzz(ji,jj,jk-1) ) 736 737 IF( cffu > zep ) THEN 738 drhou(ji,jj,jk) = 2._wp * drhox(ji+1,jj,jk) * drhox(ji,jj,jk) & 739 & / ( drhox(ji+1,jj,jk) + drhox(ji,jj,jk) ) 806 zdrho_i(ji,jj,jk ) = 0._wp 807 ENDIF 808 809 cffx = 2._wp * zdzx (ji-1,jj ,jk) * zdzx (ji,jj,jk ) 810 IF( cffx > zep ) THEN 811 zdz_i(ji,jj,jk) = cffx / ( zdzx(ji-1,jj,jk) + zdzx(ji,jj,jk) ) 740 812 ELSE 741 drhou(ji,jj,jk) = 0._wp742 ENDIF 743 744 IF( cffx > zep ) THEN745 dzu(ji,jj,jk) = 2._wp * dzx(ji+1,jj,jk) * dzx(ji,jj,jk) &746 & / ( dzx(ji+1,jj,jk) + dzx(ji,jj,jk) )813 zdz_i(ji,jj,jk) = 0._wp 814 ENDIF 815 816 cffv = 2._wp * zdrhoy(ji ,jj-1,jk) * zdrhoy(ji,jj,jk ) 817 IF( cffv > zep ) THEN 818 zdrho_j(ji,jj,jk) = cffv / ( zdrhoy(ji,jj-1,jk) + zdrhoy(ji,jj,jk) ) 747 819 ELSE 748 dzu(ji,jj,jk) = 0._wp749 ENDIF 750 751 IF( cffv > zep ) THEN752 drhov(ji,jj,jk) = 2._wp * drhoy(ji,jj+1,jk) * drhoy(ji,jj,jk) &753 & / ( drhoy(ji,jj+1,jk) + drhoy(ji,jj,jk) )820 zdrho_j(ji,jj,jk) = 0._wp 821 ENDIF 822 823 cffy = 2._wp * zdzy (ji ,jj-1,jk) * zdzy (ji,jj,jk ) 824 IF( cffy > zep ) THEN 825 zdz_j(ji,jj,jk) = cffy / ( zdzy(ji,jj-1,jk) + zdzy(ji,jj,jk) ) 754 826 ELSE 755 drhov(ji,jj,jk) = 0._wp 756 ENDIF 757 758 IF( cffy > zep ) THEN 759 dzv(ji,jj,jk) = 2._wp * dzy(ji,jj+1,jk) * dzy(ji,jj,jk) & 760 & / ( dzy(ji,jj+1,jk) + dzy(ji,jj,jk) ) 761 ELSE 762 dzv(ji,jj,jk) = 0._wp 763 ENDIF 764 765 END_3D 827 zdz_j(ji,jj,jk) = 0._wp 828 ENDIF 829 END_3D 830 831 !!! Note that zdzx, zdzy, zdzz, zdrhox, zdrhoy and zdrhoz should NOT be used beyond this point 766 832 767 833 !---------------------------------------------------------------------------------- 768 ! apply boundary conditions at top and bottomusing 5.36-5.37834 ! 6B. apply boundary conditions at side boundaries using 5.36-5.37 769 835 !---------------------------------------------------------------------------------- 770 drhow(:,:, 1 ) = 1.5_wp * ( drhoz(:,:, 2 ) - drhoz(:,:, 1 ) ) - 0.5_wp * drhow(:,:, 2 ) 771 drhou(:,:, 1 ) = 1.5_wp * ( drhox(:,:, 2 ) - drhox(:,:, 1 ) ) - 0.5_wp * drhou(:,:, 2 ) 772 drhov(:,:, 1 ) = 1.5_wp * ( drhoy(:,:, 2 ) - drhoy(:,:, 1 ) ) - 0.5_wp * drhov(:,:, 2 ) 773 774 drhow(:,:,jpk) = 1.5_wp * ( drhoz(:,:,jpk) - drhoz(:,:,jpkm1) ) - 0.5_wp * drhow(:,:,jpkm1) 775 drhou(:,:,jpk) = 1.5_wp * ( drhox(:,:,jpk) - drhox(:,:,jpkm1) ) - 0.5_wp * drhou(:,:,jpkm1) 776 drhov(:,:,jpk) = 1.5_wp * ( drhoy(:,:,jpk) - drhoy(:,:,jpkm1) ) - 0.5_wp * drhov(:,:,jpkm1) 777 836 837 DO jk = 1, jpkm1 838 zz_drho_i(:,:) = zdrho_i(:,:,jk) 839 zz_dz_i (:,:) = zdz_i (:,:,jk) 840 zz_drho_j(:,:) = zdrho_j(:,:,jk) 841 zz_dz_j (:,:) = zdz_j (:,:,jk) 842 DO_2D( 0, 1, 0, 1) 843 ! Walls coming from left: should check from 2 to jpi-1 (and jpj=2-jpj) 844 IF (ji < jpi) THEN 845 IF ( umask(ji,jj,jk) > 0.5_wp .AND. umask(ji-1,jj,jk) < 0.5_wp .AND. umask(ji+1,jj,jk) > 0.5_wp) THEN 846 zz_drho_i(ji,jj) = aco_bc_hor * ( rhd (ji+1,jj,jk) - rhd (ji,jj,jk) ) - bco_bc_hor * zdrho_i(ji+1,jj,jk) 847 zz_dz_i (ji,jj) = aco_bc_hor * (-gde3w(ji+1,jj,jk) + gde3w(ji,jj,jk) ) - bco_bc_hor * zdz_i (ji+1,jj,jk) 848 END IF 849 END IF 850 ! Walls coming from right: should check from 3 to jpi (and jpj=2-jpj) 851 IF (ji > 2) THEN 852 IF ( umask(ji,jj,jk) < 0.5_wp .AND. umask(ji-1,jj,jk) > 0.5_wp .AND. umask(ji-2,jj,jk) > 0.5_wp) THEN 853 zz_drho_i(ji,jj) = aco_bc_hor * ( rhd (ji,jj,jk) - rhd (ji-1,jj,jk) ) - bco_bc_hor * zdrho_i(ji-1,jj,jk) 854 zz_dz_i (ji,jj) = aco_bc_hor * (-gde3w(ji,jj,jk) + gde3w(ji-1,jj,jk) ) - bco_bc_hor * zdz_i (ji-1,jj,jk) 855 END IF 856 END IF 857 ! Walls coming from left: should check from 2 to jpj-1 (and jpi=2-jpi) 858 IF (jj < jpj) THEN 859 IF ( vmask(ji,jj,jk) > 0.5_wp .AND. vmask(ji,jj-1,jk) < 0.5_wp .AND. vmask(ji,jj+1,jk) > 0.5_wp) THEN 860 zz_drho_j(ji,jj) = aco_bc_hor * ( rhd (ji,jj+1,jk) - rhd (ji,jj,jk) ) - bco_bc_hor * zdrho_j(ji,jj+1,jk) 861 zz_dz_j (ji,jj) = aco_bc_hor * (-gde3w(ji,jj+1,jk) + gde3w(ji,jj,jk) ) - bco_bc_hor * zdz_j (ji,jj+1,jk) 862 END IF 863 END IF 864 ! Walls coming from right: should check from 3 to jpj (and jpi=2-jpi) 865 IF (jj > 2) THEN 866 IF ( vmask(ji,jj,jk) < 0.5_wp .AND. vmask(ji,jj-1,jk) > 0.5_wp .AND. vmask(ji,jj-2,jk) > 0.5_wp) THEN 867 zz_drho_j(ji,jj) = aco_bc_hor * ( rhd (ji,jj,jk) - rhd (ji,jj-1,jk) ) - bco_bc_hor * zdrho_j(ji,jj-1,jk) 868 zz_dz_j (ji,jj) = aco_bc_hor * (-gde3w(ji,jj,jk) + gde3w(ji,jj-1,jk) ) - bco_bc_hor * zdz_j (ji,jj-1,jk) 869 END IF 870 END IF 871 END_2D 872 zdrho_i(:,:,jk) = zz_drho_i(:,:) 873 zdz_i (:,:,jk) = zz_dz_i (:,:) 874 zdrho_j(:,:,jk) = zz_drho_j(:,:) 875 zdz_j (:,:,jk) = zz_dz_j (:,:) 876 END DO 778 877 779 878 !-------------------------------------------------------------- 780 ! Upper half of top-most grid box, compute and store879 ! 7. Calculate integrals on side faces 781 880 !------------------------------------------------------------- 782 881 783 !!bug gm : e3w-gde3w(:,:,:) = 0.5*e3w .... and gde3w(:,:,2)-gde3w(:,:,1)=e3w(:,:,2,Kmm) .... to be verified 784 ! true if gde3w(:,:,:) is really defined as the sum of the e3w scale factors as, it seems to me, it should be 785 786 DO_2D( 0, 0, 0, 0 ) 787 rho_k(ji,jj,1) = -grav * ( e3w(ji,jj,1,Kmm) - gde3w(ji,jj,1) ) & 788 & * ( rhd(ji,jj,1) & 789 & + 0.5_wp * ( rhd (ji,jj,2) - rhd (ji,jj,1) ) & 790 & * ( e3w (ji,jj,1,Kmm) - gde3w(ji,jj,1) ) & 791 & / ( gde3w(ji,jj,2) - gde3w(ji,jj,1) ) ) 792 END_2D 793 794 !!bug gm : here also, simplification is possible 795 !!bug gm : optimisation: 1/10 and 1/12 the division should be done before the loop 796 797 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 798 799 rho_k(ji,jj,jk) = zcoef0 * ( rhd (ji,jj,jk) + rhd (ji,jj,jk-1) ) & 800 & * ( gde3w(ji,jj,jk) - gde3w(ji,jj,jk-1) ) & 801 & - grav * z1_10 * ( & 802 & ( drhow (ji,jj,jk) - drhow (ji,jj,jk-1) ) & 803 & * ( gde3w(ji,jj,jk) - gde3w(ji,jj,jk-1) - z1_12 * ( dzw (ji,jj,jk) + dzw (ji,jj,jk-1) ) ) & 804 & - ( dzw (ji,jj,jk) - dzw (ji,jj,jk-1) ) & 805 & * ( rhd (ji,jj,jk) - rhd (ji,jj,jk-1) - z1_12 * ( drhow(ji,jj,jk) + drhow(ji,jj,jk-1) ) ) & 806 & ) 807 808 rho_i(ji,jj,jk) = zcoef0 * ( rhd (ji+1,jj,jk) + rhd (ji,jj,jk) ) & 809 & * ( gde3w(ji+1,jj,jk) - gde3w(ji,jj,jk) ) & 810 & - grav* z1_10 * ( & 811 & ( drhou (ji+1,jj,jk) - drhou (ji,jj,jk) ) & 812 & * ( gde3w(ji+1,jj,jk) - gde3w(ji,jj,jk) - z1_12 * ( dzu (ji+1,jj,jk) + dzu (ji,jj,jk) ) ) & 813 & - ( dzu (ji+1,jj,jk) - dzu (ji,jj,jk) ) & 814 & * ( rhd (ji+1,jj,jk) - rhd (ji,jj,jk) - z1_12 * ( drhou(ji+1,jj,jk) + drhou(ji,jj,jk) ) ) & 815 & ) 816 817 rho_j(ji,jj,jk) = zcoef0 * ( rhd (ji,jj+1,jk) + rhd (ji,jj,jk) ) & 818 & * ( gde3w(ji,jj+1,jk) - gde3w(ji,jj,jk) ) & 819 & - grav* z1_10 * ( & 820 & ( drhov (ji,jj+1,jk) - drhov (ji,jj,jk) ) & 821 & * ( gde3w(ji,jj+1,jk) - gde3w(ji,jj,jk) - z1_12 * ( dzv (ji,jj+1,jk) + dzv (ji,jj,jk) ) ) & 822 & - ( dzv (ji,jj+1,jk) - dzv (ji,jj,jk) ) & 823 & * ( rhd (ji,jj+1,jk) - rhd (ji,jj,jk) - z1_12 * ( drhov(ji,jj+1,jk) + drhov(ji,jj,jk) ) ) & 824 & ) 825 826 END_3D 827 CALL lbc_lnk_multi( 'dynhpg', rho_k, 'W', 1.0_wp, rho_i, 'U', 1.0_wp, rho_j, 'V', 1.0_wp ) 828 882 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 883 ! two -ve signs cancel in next two lines (within zcoef0 and because gde3w is a depth not a height) 884 z_rho_i(ji,jj,jk) = zcoef0 * ( rhd (ji+1,jj,jk) + rhd (ji,jj,jk) ) & 885 & * ( gde3w(ji+1,jj,jk) - gde3w(ji,jj,jk) ) 886 IF ( umask(ji-1, jj, jk) > 0.5 .OR. umask(ji+1, jj, jk) > 0.5 ) THEN 887 z_rho_i(ji,jj,jk) = z_rho_i(ji,jj,jk) - z_grav_10 * ( & 888 & ( zdrho_i (ji+1,jj,jk) - zdrho_i (ji,jj,jk) ) & 889 & * ( - gde3w(ji+1,jj,jk) + gde3w(ji,jj,jk) - z1_12 * ( zdz_i (ji+1,jj,jk) + zdz_i (ji,jj,jk) ) ) & 890 & - ( zdz_i (ji+1,jj,jk) - zdz_i (ji,jj,jk) ) & 891 & * ( rhd (ji+1,jj,jk) - rhd (ji,jj,jk) - z1_12 * ( zdrho_i(ji+1,jj,jk) + zdrho_i(ji,jj,jk) ) ) & 892 & ) 893 END IF 894 895 z_rho_j(ji,jj,jk) = zcoef0 * ( rhd (ji,jj+1,jk) + rhd (ji,jj,jk) ) & 896 & * ( gde3w(ji,jj+1,jk) - gde3w(ji,jj,jk) ) 897 IF ( vmask(ji, jj-1, jk) > 0.5 .OR. vmask(ji, jj+1, jk) > 0.5 ) THEN 898 z_rho_j(ji,jj,jk) = z_rho_j(ji,jj,jk) - z_grav_10 * ( & 899 & ( zdrho_j (ji,jj+1,jk) - zdrho_j (ji,jj,jk) ) & 900 & * ( - gde3w(ji,jj+1,jk) + gde3w(ji,jj,jk) - z1_12 * ( zdz_j (ji,jj+1,jk) + zdz_j (ji,jj,jk) ) ) & 901 & - ( zdz_j (ji,jj+1,jk) - zdz_j (ji,jj,jk) ) & 902 & * ( rhd (ji,jj+1,jk) - rhd (ji,jj,jk) - z1_12 * ( zdrho_j(ji,jj+1,jk) + zdrho_j(ji,jj,jk) ) ) & 903 & ) 904 END IF 905 END_3D 906 907 !-------------------------------------------------------------- 908 ! 8. Integrate in the vertical 909 !------------------------------------------------------------- 910 ! 829 911 ! --------------- 830 912 ! Surface value 831 913 ! --------------- 832 914 DO_2D( 0, 0, 0, 0 ) 833 zhpi(ji,jj,1) = ( rho_k(ji+1,jj ,1) - rho_k(ji,jj,1) -rho_i(ji,jj,1) ) * r1_e1u(ji,jj)834 zhpj(ji,jj,1) = ( rho_k(ji ,jj+1,1) - rho_k(ji,jj,1) -rho_j(ji,jj,1) ) * r1_e2v(ji,jj)915 zhpi(ji,jj,1) = ( z_rho_k(ji,jj,1) - z_rho_k(ji+1,jj ,1) - z_rho_i(ji,jj,1) ) * r1_e1u(ji,jj) 916 zhpj(ji,jj,1) = ( z_rho_k(ji,jj,1) - z_rho_k(ji ,jj+1,1) - z_rho_j(ji,jj,1) ) * r1_e2v(ji,jj) 835 917 IF( ln_wd_il ) THEN 836 918 zhpi(ji,jj,1) = zhpi(ji,jj,1) * zcpx(ji,jj) … … 847 929 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 848 930 ! hydrostatic pressure gradient along s-surfaces 849 zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) &850 & + ( ( rho_k(ji+1,jj,jk) - rho_k(ji,jj,jk ) )&851 & - ( rho_i(ji ,jj,jk) - rho_i(ji,jj,jk-1) ) ) * r1_e1u(ji,jj)852 zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1) &853 & + ( ( rho_k(ji,jj+1,jk) - rho_k(ji,jj,jk ) )&854 & -( rho_j(ji,jj ,jk) - rho_j(ji,jj,jk-1) ) ) * r1_e2v(ji,jj)931 zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) & 932 & + ( ( z_rho_k(ji,jj,jk) - z_rho_k(ji+1,jj,jk ) ) & 933 & - ( z_rho_i(ji,jj,jk) - z_rho_i(ji ,jj,jk-1) ) ) * r1_e1u(ji,jj) 934 zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1) & 935 & + ( ( z_rho_k(ji,jj,jk) - z_rho_k(ji,jj+1,jk ) ) & 936 & -( z_rho_j(ji,jj,jk) - z_rho_j(ji,jj ,jk-1) ) ) * r1_e2v(ji,jj) 855 937 IF( ln_wd_il ) THEN 856 938 zhpi(ji,jj,jk) = zhpi(ji,jj,jk) * zcpx(ji,jj) … … 892 974 REAL(wp) :: zrhdt1 893 975 REAL(wp) :: zdpdx1, zdpdx2, zdpdy1, zdpdy2 976 REAL(wp), DIMENSION(jpi,jpj) :: zpgu, zpgv ! 2D workspace 977 REAL(wp), DIMENSION(jpi,jpj) :: zsshu_n, zsshv_n 894 978 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdept, zrhh 895 979 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zhpi, zu, zv, fsp, xsp, asp, bsp, csp, dsp 896 REAL(wp), DIMENSION(jpi,jpj) :: zsshu_n, zsshv_n897 980 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: zcpx, zcpy !W/D pressure filter 898 981 !!---------------------------------------------------------------------- … … 907 990 zcoef0 = - grav 908 991 znad = 1._wp 909 IF( ln_linssh ) znad = 0._wp 910 992 IF( ln_linssh ) znad = 1._wp 993 ! 994 ! --------------- 995 ! Surface pressure gradient to be removed 996 ! --------------- 997 DO_2D( 0, 0, 0, 0 ) 998 zpgu(ji,jj) = - grav * ( ssh(ji+1,jj,Kmm) - ssh(ji,jj,Kmm) ) * r1_e1u(ji,jj) 999 zpgv(ji,jj) = - grav * ( ssh(ji,jj+1,Kmm) - ssh(ji,jj,Kmm) ) * r1_e2v(ji,jj) 1000 END_2D 1001 ! 911 1002 IF( ln_wd_il ) THEN 912 1003 ALLOCATE( zcpx(jpi,jpj) , zcpy(jpi,jpj) ) … … 952 1043 ENDIF 953 1044 END_2D 954 CALL lbc_lnk _multi( 'dynhpg', zcpx, 'U', 1.0_wp, zcpy, 'V', 1.0_wp )1045 CALL lbc_lnk( 'dynhpg', zcpx, 'U', 1.0_wp, zcpy, 'V', 1.0_wp ) 955 1046 ENDIF 956 1047 … … 974 1065 ! Transfer the depth of "T(:,:,:)" to vertical coordinate "zdept(:,:,:)" 975 1066 DO_2D( 1, 1, 1, 1 ) 976 zdept(ji,jj,1) = 0.5_wp * e3w(ji,jj,1,Kmm) - ssh(ji,jj,Kmm) * znad1067 zdept(ji,jj,1) = 0.5_wp * e3w(ji,jj,1,Kmm) - ssh(ji,jj,Kmm) 977 1068 END_2D 978 1069 … … 1022 1113 END_2D 1023 1114 1024 CALL lbc_lnk _multi('dynhpg', zsshu_n, 'U', 1.0_wp, zsshv_n, 'V', 1.0_wp )1115 CALL lbc_lnk ('dynhpg', zsshu_n, 'U', 1.0_wp, zsshv_n, 'V', 1.0_wp ) 1025 1116 1026 1117 DO_2D( 0, 0, 0, 0 ) 1027 zu(ji,jj,1) = - ( e3u(ji,jj,1,Kmm) - zsshu_n(ji,jj) * znad)1028 zv(ji,jj,1) = - ( e3v(ji,jj,1,Kmm) - zsshv_n(ji,jj) * znad)1118 zu(ji,jj,1) = - ( e3u(ji,jj,1,Kmm) - zsshu_n(ji,jj) ) 1119 zv(ji,jj,1) = - ( e3v(ji,jj,1,Kmm) - zsshv_n(ji,jj) ) 1029 1120 END_2D 1030 1121 … … 1108 1199 zdpdx2 = zdpdx2 * zcpx(ji,jj) * wdrampu(ji,jj) 1109 1200 ENDIF 1110 puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) + (zdpdx1 + zdpdx2 ) * umask(ji,jj,jk)1201 puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) + (zdpdx1 + zdpdx2 - zpgu(ji,jj)) * umask(ji,jj,jk) 1111 1202 ENDIF 1112 1203 … … 1168 1259 ENDIF 1169 1260 1170 pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) + (zdpdy1 + zdpdy2 ) * vmask(ji,jj,jk)1261 pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) + (zdpdy1 + zdpdy2 - zpgv(ji,jj)) * vmask(ji,jj,jk) 1171 1262 ENDIF 1172 1263 !
Note: See TracChangeset
for help on using the changeset viewer.
