- Timestamp:
- 2019-11-22T15:29:17+01:00 (4 years ago)
- Location:
- NEMO/branches/2019/dev_r11943_MERGE_2019/src
- Files:
-
- 2 edited
-
. (modified) (1 prop)
-
OCE/DYN/dynhpg.F90 (modified) (64 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2019/dev_r11943_MERGE_2019/src
- Property svn:mergeinfo deleted
-
NEMO/branches/2019/dev_r11943_MERGE_2019/src/OCE/DYN/dynhpg.F90
r11536 r11949 81 81 CONTAINS 82 82 83 SUBROUTINE dyn_hpg( kt )83 SUBROUTINE dyn_hpg( kt, Kmm, puu, pvv, Krhs ) 84 84 !!--------------------------------------------------------------------- 85 85 !! *** ROUTINE dyn_hpg *** … … 88 88 !! using the scheme defined in the namelist 89 89 !! 90 !! ** Action : - Update ( ua,va) with the now hydrastatic pressure trend90 !! ** Action : - Update (puu(:,:,:,Krhs),pvv(:,:,:,Krhs)) with the now hydrastatic pressure trend 91 91 !! - send trends to trd_dyn for futher diagnostics (l_trddyn=T) 92 92 !!---------------------------------------------------------------------- 93 INTEGER, INTENT(in) :: kt ! ocean time-step index 93 INTEGER , INTENT( in ) :: kt ! ocean time-step index 94 INTEGER , INTENT( in ) :: Kmm, Krhs ! ocean time level indices 95 REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) :: puu, pvv ! ocean velocities and RHS of momentum equation 96 ! 94 97 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: ztrdu, ztrdv 95 98 !!---------------------------------------------------------------------- … … 97 100 IF( ln_timing ) CALL timing_start('dyn_hpg') 98 101 ! 99 IF( l_trddyn ) THEN ! Temporary saving of ua and vatrends (l_trddyn)102 IF( l_trddyn ) THEN ! Temporary saving of puu(:,:,:,Krhs) and pvv(:,:,:,Krhs) trends (l_trddyn) 100 103 ALLOCATE( ztrdu(jpi,jpj,jpk) , ztrdv(jpi,jpj,jpk) ) 101 ztrdu(:,:,:) = ua(:,:,:)102 ztrdv(:,:,:) = va(:,:,:)104 ztrdu(:,:,:) = puu(:,:,:,Krhs) 105 ztrdv(:,:,:) = pvv(:,:,:,Krhs) 103 106 ENDIF 104 107 ! 105 108 SELECT CASE ( nhpg ) ! Hydrostatic pressure gradient computation 106 CASE ( np_zco ) ; CALL hpg_zco ( kt )! z-coordinate107 CASE ( np_zps ) ; CALL hpg_zps ( kt )! z-coordinate plus partial steps (interpolation)108 CASE ( np_sco ) ; CALL hpg_sco ( kt )! s-coordinate (standard jacobian formulation)109 CASE ( np_djc ) ; CALL hpg_djc ( kt )! s-coordinate (Density Jacobian with Cubic polynomial)110 CASE ( np_prj ) ; CALL hpg_prj ( kt )! s-coordinate (Pressure Jacobian scheme)111 CASE ( np_isf ) ; CALL hpg_isf ( kt )! s-coordinate similar to sco modify for ice shelf109 CASE ( np_zco ) ; CALL hpg_zco ( kt, Kmm, puu, pvv, Krhs ) ! z-coordinate 110 CASE ( np_zps ) ; CALL hpg_zps ( kt, Kmm, puu, pvv, Krhs ) ! z-coordinate plus partial steps (interpolation) 111 CASE ( np_sco ) ; CALL hpg_sco ( kt, Kmm, puu, pvv, Krhs ) ! s-coordinate (standard jacobian formulation) 112 CASE ( np_djc ) ; CALL hpg_djc ( kt, Kmm, puu, pvv, Krhs ) ! s-coordinate (Density Jacobian with Cubic polynomial) 113 CASE ( np_prj ) ; CALL hpg_prj ( kt, Kmm, puu, pvv, Krhs ) ! s-coordinate (Pressure Jacobian scheme) 114 CASE ( np_isf ) ; CALL hpg_isf ( kt, Kmm, puu, pvv, Krhs ) ! s-coordinate similar to sco modify for ice shelf 112 115 END SELECT 113 116 ! 114 117 IF( l_trddyn ) THEN ! save the hydrostatic pressure gradient trends for momentum trend diagnostics 115 ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)116 ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)117 CALL trd_dyn( ztrdu, ztrdv, jpdyn_hpg, kt )118 ztrdu(:,:,:) = puu(:,:,:,Krhs) - ztrdu(:,:,:) 119 ztrdv(:,:,:) = pvv(:,:,:,Krhs) - ztrdv(:,:,:) 120 CALL trd_dyn( ztrdu, ztrdv, jpdyn_hpg, kt, Kmm ) 118 121 DEALLOCATE( ztrdu , ztrdv ) 119 122 ENDIF 120 123 ! 121 IF(ln_ctl) CALL prt_ctl( tab3d_1= ua, clinfo1=' hpg - Ua: ', mask1=umask, &122 & tab3d_2= va, clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' )124 IF(ln_ctl) CALL prt_ctl( tab3d_1=puu(:,:,:,Krhs), clinfo1=' hpg - Ua: ', mask1=umask, & 125 & tab3d_2=pvv(:,:,:,Krhs), clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' ) 123 126 ! 124 127 IF( ln_timing ) CALL timing_stop('dyn_hpg') … … 127 130 128 131 129 SUBROUTINE dyn_hpg_init 132 SUBROUTINE dyn_hpg_init( Kmm ) 130 133 !!---------------------------------------------------------------------- 131 134 !! *** ROUTINE dyn_hpg_init *** … … 137 140 !! with the type of vertical coordinate used (zco, zps, sco) 138 141 !!---------------------------------------------------------------------- 142 INTEGER, INTENT( in ) :: Kmm ! ocean time level index 143 ! 139 144 INTEGER :: ioptio = 0 ! temporary integer 140 145 INTEGER :: ios ! Local integer output status for namelist read … … 228 233 ! 229 234 DO jk = 1, jpk !- compute density of the water displaced by the ice shelf 230 CALL eos( zts_top(:,:,:), gdept _n(:,:,jk), zrhd(:,:,jk) )235 CALL eos( zts_top(:,:,:), gdept(:,:,jk,Kmm), zrhd(:,:,jk) ) 231 236 END DO 232 237 ! … … 239 244 DO ji = 1, jpi ! divided by 2 later 240 245 ikt = mikt(ji,jj) 241 ziceload(ji,jj) = ziceload(ji,jj) + (znad + zrhd(ji,jj,1) ) * e3w _n(ji,jj,1) * (1._wp - tmask(ji,jj,1))246 ziceload(ji,jj) = ziceload(ji,jj) + (znad + zrhd(ji,jj,1) ) * e3w(ji,jj,1,Kmm) * (1._wp - tmask(ji,jj,1)) 242 247 DO jk = 2, ikt-1 243 ziceload(ji,jj) = ziceload(ji,jj) + (2._wp * znad + zrhd(ji,jj,jk-1) + zrhd(ji,jj,jk)) * e3w _n(ji,jj,jk) &248 ziceload(ji,jj) = ziceload(ji,jj) + (2._wp * znad + zrhd(ji,jj,jk-1) + zrhd(ji,jj,jk)) * e3w(ji,jj,jk,Kmm) & 244 249 & * (1._wp - tmask(ji,jj,jk)) 245 250 END DO 246 251 IF (ikt >= 2) ziceload(ji,jj) = ziceload(ji,jj) + (2._wp * znad + zrhdtop_isf(ji,jj) + zrhd(ji,jj,ikt-1)) & 247 & * ( risfdep(ji,jj) - gdept _n(ji,jj,ikt-1) )252 & * ( risfdep(ji,jj) - gdept(ji,jj,ikt-1,Kmm) ) 248 253 END DO 249 254 END DO … … 256 261 257 262 258 SUBROUTINE hpg_zco( kt )263 SUBROUTINE hpg_zco( kt, Kmm, puu, pvv, Krhs ) 259 264 !!--------------------------------------------------------------------- 260 265 !! *** ROUTINE hpg_zco *** … … 266 271 !! level: zhpi = grav ..... 267 272 !! zhpj = grav ..... 268 !! add it to the general momentum trend (ua,va). 269 !! ua = ua - 1/e1u * zhpi 270 !! va = va - 1/e2v * zhpj 271 !! 272 !! ** Action : - Update (ua,va) with the now hydrastatic pressure trend 273 !!---------------------------------------------------------------------- 274 INTEGER, INTENT(in) :: kt ! ocean time-step index 273 !! add it to the general momentum trend (puu(:,:,:,Krhs),pvv(:,:,:,Krhs)). 274 !! puu(:,:,:,Krhs) = puu(:,:,:,Krhs) - 1/e1u * zhpi 275 !! pvv(:,:,:,Krhs) = pvv(:,:,:,Krhs) - 1/e2v * zhpj 276 !! 277 !! ** Action : - Update (puu(:,:,:,Krhs),pvv(:,:,:,Krhs)) with the now hydrastatic pressure trend 278 !!---------------------------------------------------------------------- 279 INTEGER , INTENT( in ) :: kt ! ocean time-step index 280 INTEGER , INTENT( in ) :: Kmm, Krhs ! ocean time level indices 281 REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) :: puu, pvv ! ocean velocities and RHS of momentum equation 275 282 ! 276 283 INTEGER :: ji, jj, jk ! dummy loop indices … … 290 297 DO jj = 2, jpjm1 291 298 DO ji = fs_2, fs_jpim1 ! vector opt. 292 zcoef1 = zcoef0 * e3w _n(ji,jj,1)299 zcoef1 = zcoef0 * e3w(ji,jj,1,Kmm) 293 300 ! hydrostatic pressure gradient 294 301 zhpi(ji,jj,1) = zcoef1 * ( rhd(ji+1,jj,1) - rhd(ji,jj,1) ) * r1_e1u(ji,jj) 295 302 zhpj(ji,jj,1) = zcoef1 * ( rhd(ji,jj+1,1) - rhd(ji,jj,1) ) * r1_e2v(ji,jj) 296 303 ! add to the general momentum trend 297 ua(ji,jj,1) = ua(ji,jj,1) + zhpi(ji,jj,1)298 va(ji,jj,1) = va(ji,jj,1) + zhpj(ji,jj,1)304 puu(ji,jj,1,Krhs) = puu(ji,jj,1,Krhs) + zhpi(ji,jj,1) 305 pvv(ji,jj,1,Krhs) = pvv(ji,jj,1,Krhs) + zhpj(ji,jj,1) 299 306 END DO 300 307 END DO … … 305 312 DO jj = 2, jpjm1 306 313 DO ji = fs_2, fs_jpim1 ! vector opt. 307 zcoef1 = zcoef0 * e3w _n(ji,jj,jk)314 zcoef1 = zcoef0 * e3w(ji,jj,jk,Kmm) 308 315 ! hydrostatic pressure gradient 309 316 zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) & … … 315 322 & - ( rhd(ji,jj, jk)+rhd(ji,jj ,jk-1) ) ) * r1_e2v(ji,jj) 316 323 ! add to the general momentum trend 317 ua(ji,jj,jk) = ua(ji,jj,jk) + zhpi(ji,jj,jk)318 va(ji,jj,jk) = va(ji,jj,jk) + zhpj(ji,jj,jk)324 puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) + zhpi(ji,jj,jk) 325 pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) + zhpj(ji,jj,jk) 319 326 END DO 320 327 END DO … … 324 331 325 332 326 SUBROUTINE hpg_zps( kt )333 SUBROUTINE hpg_zps( kt, Kmm, puu, pvv, Krhs ) 327 334 !!--------------------------------------------------------------------- 328 335 !! *** ROUTINE hpg_zps *** … … 330 337 !! ** Method : z-coordinate plus partial steps case. blahblah... 331 338 !! 332 !! ** Action : - Update (ua,va) with the now hydrastatic pressure trend 333 !!---------------------------------------------------------------------- 334 INTEGER, INTENT(in) :: kt ! ocean time-step index 339 !! ** Action : - Update (puu(:,:,:,Krhs),pvv(:,:,:,Krhs)) with the now hydrastatic pressure trend 340 !!---------------------------------------------------------------------- 341 INTEGER , INTENT( in ) :: kt ! ocean time-step index 342 INTEGER , INTENT( in ) :: Kmm, Krhs ! ocean time level indices 343 REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) :: puu, pvv ! ocean velocities and RHS of momentum equation 335 344 !! 336 345 INTEGER :: ji, jj, jk ! dummy loop indices … … 348 357 349 358 ! Partial steps: Compute NOW horizontal gradient of t, s, rd at the last ocean level 350 CALL zps_hde( kt, jpts, tsn, zgtsu, zgtsv, rhd, zgru , zgrv )359 CALL zps_hde( kt, Kmm, jpts, ts(:,:,:,:,Kmm), zgtsu, zgtsv, rhd, zgru , zgrv ) 351 360 352 361 ! Local constant initialization … … 356 365 DO jj = 2, jpjm1 357 366 DO ji = fs_2, fs_jpim1 ! vector opt. 358 zcoef1 = zcoef0 * e3w _n(ji,jj,1)367 zcoef1 = zcoef0 * e3w(ji,jj,1,Kmm) 359 368 ! hydrostatic pressure gradient 360 369 zhpi(ji,jj,1) = zcoef1 * ( rhd(ji+1,jj ,1) - rhd(ji,jj,1) ) * r1_e1u(ji,jj) 361 370 zhpj(ji,jj,1) = zcoef1 * ( rhd(ji ,jj+1,1) - rhd(ji,jj,1) ) * r1_e2v(ji,jj) 362 371 ! add to the general momentum trend 363 ua(ji,jj,1) = ua(ji,jj,1) + zhpi(ji,jj,1)364 va(ji,jj,1) = va(ji,jj,1) + zhpj(ji,jj,1)372 puu(ji,jj,1,Krhs) = puu(ji,jj,1,Krhs) + zhpi(ji,jj,1) 373 pvv(ji,jj,1,Krhs) = pvv(ji,jj,1,Krhs) + zhpj(ji,jj,1) 365 374 END DO 366 375 END DO … … 370 379 DO jj = 2, jpjm1 371 380 DO ji = fs_2, fs_jpim1 ! vector opt. 372 zcoef1 = zcoef0 * e3w _n(ji,jj,jk)381 zcoef1 = zcoef0 * e3w(ji,jj,jk,Kmm) 373 382 ! hydrostatic pressure gradient 374 383 zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) & … … 380 389 & - ( rhd(ji,jj, jk) + rhd(ji,jj ,jk-1) ) ) * r1_e2v(ji,jj) 381 390 ! add to the general momentum trend 382 ua(ji,jj,jk) = ua(ji,jj,jk) + zhpi(ji,jj,jk)383 va(ji,jj,jk) = va(ji,jj,jk) + zhpj(ji,jj,jk)391 puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) + zhpi(ji,jj,jk) 392 pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) + zhpj(ji,jj,jk) 384 393 END DO 385 394 END DO … … 391 400 iku = mbku(ji,jj) 392 401 ikv = mbkv(ji,jj) 393 zcoef2 = zcoef0 * MIN( e3w _n(ji,jj,iku), e3w_n(ji+1,jj ,iku) )394 zcoef3 = zcoef0 * MIN( e3w _n(ji,jj,ikv), e3w_n(ji ,jj+1,ikv) )402 zcoef2 = zcoef0 * MIN( e3w(ji,jj,iku,Kmm), e3w(ji+1,jj ,iku,Kmm) ) 403 zcoef3 = zcoef0 * MIN( e3w(ji,jj,ikv,Kmm), e3w(ji ,jj+1,ikv,Kmm) ) 395 404 IF( iku > 1 ) THEN ! on i-direction (level 2 or more) 396 ua (ji,jj,iku) = ua(ji,jj,iku) - zhpi(ji,jj,iku) ! subtract old value405 puu (ji,jj,iku,Krhs) = puu(ji,jj,iku,Krhs) - zhpi(ji,jj,iku) ! subtract old value 397 406 zhpi(ji,jj,iku) = zhpi(ji,jj,iku-1) & ! compute the new one 398 407 & + zcoef2 * ( rhd(ji+1,jj,iku-1) - rhd(ji,jj,iku-1) + zgru(ji,jj) ) * r1_e1u(ji,jj) 399 ua (ji,jj,iku) = ua(ji,jj,iku) + zhpi(ji,jj,iku) ! add the new one to the general momentum trend408 puu (ji,jj,iku,Krhs) = puu(ji,jj,iku,Krhs) + zhpi(ji,jj,iku) ! add the new one to the general momentum trend 400 409 ENDIF 401 410 IF( ikv > 1 ) THEN ! on j-direction (level 2 or more) 402 va (ji,jj,ikv) = va(ji,jj,ikv) - zhpj(ji,jj,ikv) ! subtract old value411 pvv (ji,jj,ikv,Krhs) = pvv(ji,jj,ikv,Krhs) - zhpj(ji,jj,ikv) ! subtract old value 403 412 zhpj(ji,jj,ikv) = zhpj(ji,jj,ikv-1) & ! compute the new one 404 413 & + zcoef3 * ( rhd(ji,jj+1,ikv-1) - rhd(ji,jj,ikv-1) + zgrv(ji,jj) ) * r1_e2v(ji,jj) 405 va (ji,jj,ikv) = va(ji,jj,ikv) + zhpj(ji,jj,ikv) ! add the new one to the general momentum trend414 pvv (ji,jj,ikv,Krhs) = pvv(ji,jj,ikv,Krhs) + zhpj(ji,jj,ikv) ! add the new one to the general momentum trend 406 415 ENDIF 407 416 END DO … … 411 420 412 421 413 SUBROUTINE hpg_sco( kt )422 SUBROUTINE hpg_sco( kt, Kmm, puu, pvv, Krhs ) 414 423 !!--------------------------------------------------------------------- 415 424 !! *** ROUTINE hpg_sco *** … … 423 432 !! zhpi = grav ..... + 1/e1u mi(rhd) di[ grav dep3w ] 424 433 !! zhpj = grav ..... + 1/e2v mj(rhd) dj[ grav dep3w ] 425 !! add it to the general momentum trend (ua,va). 426 !! ua = ua - 1/e1u * zhpi 427 !! va = va - 1/e2v * zhpj 428 !! 429 !! ** Action : - Update (ua,va) with the now hydrastatic pressure trend 430 !!---------------------------------------------------------------------- 431 INTEGER, INTENT(in) :: kt ! ocean time-step index 434 !! add it to the general momentum trend (puu(:,:,:,Krhs),pvv(:,:,:,Krhs)). 435 !! puu(:,:,:,Krhs) = puu(:,:,:,Krhs) - 1/e1u * zhpi 436 !! pvv(:,:,:,Krhs) = pvv(:,:,:,Krhs) - 1/e2v * zhpj 437 !! 438 !! ** Action : - Update (puu(:,:,:,Krhs),pvv(:,:,:,Krhs)) with the now hydrastatic pressure trend 439 !!---------------------------------------------------------------------- 440 INTEGER , INTENT( in ) :: kt ! ocean time-step index 441 INTEGER , INTENT( in ) :: Kmm, Krhs ! ocean time level indices 442 REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) :: puu, pvv ! ocean velocities and RHS of momentum equation 432 443 !! 433 444 INTEGER :: ji, jj, jk, jii, jjj ! dummy loop indices … … 454 465 DO jj = 2, jpjm1 455 466 DO ji = 2, jpim1 456 ll_tmp1 = MIN( ssh n(ji,jj) , sshn(ji+1,jj) ) > &467 ll_tmp1 = MIN( ssh(ji,jj,Kmm) , ssh(ji+1,jj,Kmm) ) > & 457 468 & MAX( -ht_0(ji,jj) , -ht_0(ji+1,jj) ) .AND. & 458 & MAX( ssh n(ji,jj) + ht_0(ji,jj), sshn(ji+1,jj) + ht_0(ji+1,jj) ) &469 & MAX( ssh(ji,jj,Kmm) + ht_0(ji,jj), ssh(ji+1,jj,Kmm) + ht_0(ji+1,jj) ) & 459 470 & > rn_wdmin1 + rn_wdmin2 460 ll_tmp2 = ( ABS( ssh n(ji,jj) - sshn(ji+1,jj) ) > 1.E-12 ) .AND. ( &461 & MAX( ssh n(ji,jj) , sshn(ji+1,jj) ) > &471 ll_tmp2 = ( ABS( ssh(ji,jj,Kmm) - ssh(ji+1,jj,Kmm) ) > 1.E-12 ) .AND. ( & 472 & MAX( ssh(ji,jj,Kmm) , ssh(ji+1,jj,Kmm) ) > & 462 473 & MAX( -ht_0(ji,jj) , -ht_0(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 ) 463 474 … … 465 476 zcpx(ji,jj) = 1.0_wp 466 477 ELSE IF(ll_tmp2) THEN 467 ! no worries about ssh n(ji+1,jj) - sshn(ji ,jj) = 0, it won't happen ! here468 zcpx(ji,jj) = ABS( (ssh n(ji+1,jj) + ht_0(ji+1,jj) - sshn(ji,jj) - ht_0(ji,jj)) &469 & / (ssh n(ji+1,jj) - sshn(ji ,jj)) )478 ! no worries about ssh(ji+1,jj,Kmm) - ssh(ji ,jj,Kmm) = 0, it won't happen ! here 479 zcpx(ji,jj) = ABS( (ssh(ji+1,jj,Kmm) + ht_0(ji+1,jj) - ssh(ji,jj,Kmm) - ht_0(ji,jj)) & 480 & / (ssh(ji+1,jj,Kmm) - ssh(ji ,jj,Kmm)) ) 470 481 ELSE 471 482 zcpx(ji,jj) = 0._wp 472 483 END IF 473 484 474 ll_tmp1 = MIN( ssh n(ji,jj) , sshn(ji,jj+1) ) > &485 ll_tmp1 = MIN( ssh(ji,jj,Kmm) , ssh(ji,jj+1,Kmm) ) > & 475 486 & MAX( -ht_0(ji,jj) , -ht_0(ji,jj+1) ) .AND. & 476 & MAX( ssh n(ji,jj) + ht_0(ji,jj), sshn(ji,jj+1) + ht_0(ji,jj+1) ) &487 & MAX( ssh(ji,jj,Kmm) + ht_0(ji,jj), ssh(ji,jj+1,Kmm) + ht_0(ji,jj+1) ) & 477 488 & > rn_wdmin1 + rn_wdmin2 478 ll_tmp2 = ( ABS( ssh n(ji,jj) - sshn(ji,jj+1) ) > 1.E-12 ) .AND. ( &479 & MAX( ssh n(ji,jj) , sshn(ji,jj+1) ) > &489 ll_tmp2 = ( ABS( ssh(ji,jj,Kmm) - ssh(ji,jj+1,Kmm) ) > 1.E-12 ) .AND. ( & 490 & MAX( ssh(ji,jj,Kmm) , ssh(ji,jj+1,Kmm) ) > & 480 491 & MAX( -ht_0(ji,jj) , -ht_0(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 ) 481 492 … … 483 494 zcpy(ji,jj) = 1.0_wp 484 495 ELSE IF(ll_tmp2) THEN 485 ! no worries about ssh n(ji,jj+1) - sshn(ji,jj) = 0, it won't happen ! here486 zcpy(ji,jj) = ABS( (ssh n(ji,jj+1) + ht_0(ji,jj+1) - sshn(ji,jj) - ht_0(ji,jj)) &487 & / (ssh n(ji,jj+1) - sshn(ji,jj)) )496 ! no worries about ssh(ji,jj+1,Kmm) - ssh(ji,jj ,Kmm) = 0, it won't happen ! here 497 zcpy(ji,jj) = ABS( (ssh(ji,jj+1,Kmm) + ht_0(ji,jj+1) - ssh(ji,jj,Kmm) - ht_0(ji,jj)) & 498 & / (ssh(ji,jj+1,Kmm) - ssh(ji,jj ,Kmm)) ) 488 499 ELSE 489 500 zcpy(ji,jj) = 0._wp … … 498 509 DO ji = fs_2, fs_jpim1 ! vector opt. 499 510 ! hydrostatic pressure gradient along s-surfaces 500 zhpi(ji,jj,1) = zcoef0 * ( e3w _n(ji+1,jj ,1) * ( znad + rhd(ji+1,jj ,1) ) &501 & - e3w _n(ji ,jj ,1) * ( znad + rhd(ji ,jj ,1) ) ) * r1_e1u(ji,jj)502 zhpj(ji,jj,1) = zcoef0 * ( e3w _n(ji ,jj+1,1) * ( znad + rhd(ji ,jj+1,1) ) &503 & - e3w _n(ji ,jj ,1) * ( znad + rhd(ji ,jj ,1) ) ) * r1_e2v(ji,jj)511 zhpi(ji,jj,1) = zcoef0 * ( e3w(ji+1,jj ,1,Kmm) * ( znad + rhd(ji+1,jj ,1) ) & 512 & - e3w(ji ,jj ,1,Kmm) * ( znad + rhd(ji ,jj ,1) ) ) * r1_e1u(ji,jj) 513 zhpj(ji,jj,1) = zcoef0 * ( e3w(ji ,jj+1,1,Kmm) * ( znad + rhd(ji ,jj+1,1) ) & 514 & - e3w(ji ,jj ,1,Kmm) * ( znad + rhd(ji ,jj ,1) ) ) * r1_e2v(ji,jj) 504 515 ! s-coordinate pressure gradient correction 505 516 zuap = -zcoef0 * ( rhd (ji+1,jj,1) + rhd (ji,jj,1) + 2._wp * znad ) & 506 & * ( gde3w _n(ji+1,jj,1) - gde3w_n(ji,jj,1) ) * r1_e1u(ji,jj)517 & * ( gde3w(ji+1,jj,1) - gde3w(ji,jj,1) ) * r1_e1u(ji,jj) 507 518 zvap = -zcoef0 * ( rhd (ji,jj+1,1) + rhd (ji,jj,1) + 2._wp * znad ) & 508 & * ( gde3w _n(ji,jj+1,1) - gde3w_n(ji,jj,1) ) * r1_e2v(ji,jj)519 & * ( gde3w(ji,jj+1,1) - gde3w(ji,jj,1) ) * r1_e2v(ji,jj) 509 520 ! 510 521 IF( ln_wd_il ) THEN … … 516 527 ! 517 528 ! add to the general momentum trend 518 ua(ji,jj,1) = ua(ji,jj,1) + zhpi(ji,jj,1) + zuap519 va(ji,jj,1) = va(ji,jj,1) + zhpj(ji,jj,1) + zvap529 puu(ji,jj,1,Krhs) = puu(ji,jj,1,Krhs) + zhpi(ji,jj,1) + zuap 530 pvv(ji,jj,1,Krhs) = pvv(ji,jj,1,Krhs) + zhpj(ji,jj,1) + zvap 520 531 END DO 521 532 END DO … … 527 538 ! hydrostatic pressure gradient along s-surfaces 528 539 zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) + zcoef0 * r1_e1u(ji,jj) & 529 & * ( e3w _n(ji+1,jj,jk) * ( rhd(ji+1,jj,jk) + rhd(ji+1,jj,jk-1) + 2*znad ) &530 & - e3w _n(ji ,jj,jk) * ( rhd(ji ,jj,jk) + rhd(ji ,jj,jk-1) + 2*znad ) )540 & * ( e3w(ji+1,jj,jk,Kmm) * ( rhd(ji+1,jj,jk) + rhd(ji+1,jj,jk-1) + 2*znad ) & 541 & - e3w(ji ,jj,jk,Kmm) * ( rhd(ji ,jj,jk) + rhd(ji ,jj,jk-1) + 2*znad ) ) 531 542 zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1) + zcoef0 * r1_e2v(ji,jj) & 532 & * ( e3w _n(ji,jj+1,jk) * ( rhd(ji,jj+1,jk) + rhd(ji,jj+1,jk-1) + 2*znad ) &533 & - e3w _n(ji,jj ,jk) * ( rhd(ji,jj, jk) + rhd(ji,jj ,jk-1) + 2*znad ) )543 & * ( e3w(ji,jj+1,jk,Kmm) * ( rhd(ji,jj+1,jk) + rhd(ji,jj+1,jk-1) + 2*znad ) & 544 & - e3w(ji,jj ,jk,Kmm) * ( rhd(ji,jj, jk) + rhd(ji,jj ,jk-1) + 2*znad ) ) 534 545 ! s-coordinate pressure gradient correction 535 546 zuap = -zcoef0 * ( rhd (ji+1,jj ,jk) + rhd (ji,jj,jk) + 2._wp * znad ) & 536 & * ( gde3w _n(ji+1,jj ,jk) - gde3w_n(ji,jj,jk) ) * r1_e1u(ji,jj)547 & * ( gde3w(ji+1,jj ,jk) - gde3w(ji,jj,jk) ) * r1_e1u(ji,jj) 537 548 zvap = -zcoef0 * ( rhd (ji ,jj+1,jk) + rhd (ji,jj,jk) + 2._wp * znad ) & 538 & * ( gde3w _n(ji ,jj+1,jk) - gde3w_n(ji,jj,jk) ) * r1_e2v(ji,jj)549 & * ( gde3w(ji ,jj+1,jk) - gde3w(ji,jj,jk) ) * r1_e2v(ji,jj) 539 550 ! 540 551 IF( ln_wd_il ) THEN … … 546 557 ! 547 558 ! add to the general momentum trend 548 ua(ji,jj,jk) = ua(ji,jj,jk) + zhpi(ji,jj,jk) + zuap549 va(ji,jj,jk) = va(ji,jj,jk) + zhpj(ji,jj,jk) + zvap559 puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) + zhpi(ji,jj,jk) + zuap 560 pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) + zhpj(ji,jj,jk) + zvap 550 561 END DO 551 562 END DO … … 557 568 558 569 559 SUBROUTINE hpg_isf( kt )570 SUBROUTINE hpg_isf( kt, Kmm, puu, pvv, Krhs ) 560 571 !!--------------------------------------------------------------------- 561 572 !! *** ROUTINE hpg_isf *** … … 569 580 !! zhpi = grav ..... + 1/e1u mi(rhd) di[ grav dep3w ] 570 581 !! zhpj = grav ..... + 1/e2v mj(rhd) dj[ grav dep3w ] 571 !! add it to the general momentum trend ( ua,va).572 !! ua = ua- 1/e1u * zhpi573 !! va = va- 1/e2v * zhpj582 !! add it to the general momentum trend (puu(:,:,:,Krhs),pvv(:,:,:,Krhs)). 583 !! puu(:,:,:,Krhs) = puu(:,:,:,Krhs) - 1/e1u * zhpi 584 !! pvv(:,:,:,Krhs) = pvv(:,:,:,Krhs) - 1/e2v * zhpj 574 585 !! iceload is added and partial cell case are added to the top and bottom 575 586 !! 576 !! ** Action : - Update (ua,va) with the now hydrastatic pressure trend 577 !!---------------------------------------------------------------------- 578 INTEGER, INTENT(in) :: kt ! ocean time-step index 587 !! ** Action : - Update (puu(:,:,:,Krhs),pvv(:,:,:,Krhs)) with the now hydrastatic pressure trend 588 !!---------------------------------------------------------------------- 589 INTEGER , INTENT( in ) :: kt ! ocean time-step index 590 INTEGER , INTENT( in ) :: Kmm, Krhs ! ocean time level indices 591 REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) :: puu, pvv ! ocean velocities and RHS of momentum equation 579 592 !! 580 593 INTEGER :: ji, jj, jk, ikt, iktp1i, iktp1j ! dummy loop indices … … 597 610 DO jj = 1, jpj 598 611 ikt = mikt(ji,jj) 599 zts_top(ji,jj,1) = ts n(ji,jj,ikt,1)600 zts_top(ji,jj,2) = ts n(ji,jj,ikt,2)612 zts_top(ji,jj,1) = ts(ji,jj,ikt,1,Kmm) 613 zts_top(ji,jj,2) = ts(ji,jj,ikt,2,Kmm) 601 614 END DO 602 615 END DO … … 613 626 ! hydrostatic pressure gradient along s-surfaces and ice shelf pressure 614 627 ! we assume ISF is in isostatic equilibrium 615 zhpi(ji,jj,1) = zcoef0 / e1u(ji,jj) * ( 0.5_wp * e3w _n(ji+1,jj,iktp1i) &628 zhpi(ji,jj,1) = zcoef0 / e1u(ji,jj) * ( 0.5_wp * e3w(ji+1,jj,iktp1i,Kmm) & 616 629 & * ( 2._wp * znad + rhd(ji+1,jj,iktp1i) + zrhdtop_oce(ji+1,jj) ) & 617 & - 0.5_wp * e3w _n(ji,jj,ikt) &630 & - 0.5_wp * e3w(ji,jj,ikt,Kmm) & 618 631 & * ( 2._wp * znad + rhd(ji,jj,ikt) + zrhdtop_oce(ji,jj) ) & 619 632 & + ( riceload(ji+1,jj) - riceload(ji,jj)) ) 620 zhpj(ji,jj,1) = zcoef0 / e2v(ji,jj) * ( 0.5_wp * e3w _n(ji,jj+1,iktp1j) &633 zhpj(ji,jj,1) = zcoef0 / e2v(ji,jj) * ( 0.5_wp * e3w(ji,jj+1,iktp1j,Kmm) & 621 634 & * ( 2._wp * znad + rhd(ji,jj+1,iktp1j) + zrhdtop_oce(ji,jj+1) ) & 622 & - 0.5_wp * e3w _n(ji,jj,ikt) &635 & - 0.5_wp * e3w(ji,jj,ikt,Kmm) & 623 636 & * ( 2._wp * znad + rhd(ji,jj,ikt) + zrhdtop_oce(ji,jj) ) & 624 637 & + ( riceload(ji,jj+1) - riceload(ji,jj)) ) 625 638 ! s-coordinate pressure gradient correction (=0 if z coordinate) 626 639 zuap = -zcoef0 * ( rhd (ji+1,jj,1) + rhd (ji,jj,1) + 2._wp * znad ) & 627 & * ( gde3w _n(ji+1,jj,1) - gde3w_n(ji,jj,1) ) * r1_e1u(ji,jj)640 & * ( gde3w(ji+1,jj,1) - gde3w(ji,jj,1) ) * r1_e1u(ji,jj) 628 641 zvap = -zcoef0 * ( rhd (ji,jj+1,1) + rhd (ji,jj,1) + 2._wp * znad ) & 629 & * ( gde3w _n(ji,jj+1,1) - gde3w_n(ji,jj,1) ) * r1_e2v(ji,jj)642 & * ( gde3w(ji,jj+1,1) - gde3w(ji,jj,1) ) * r1_e2v(ji,jj) 630 643 ! add to the general momentum trend 631 ua(ji,jj,1) = ua(ji,jj,1) + (zhpi(ji,jj,1) + zuap) * umask(ji,jj,1)632 va(ji,jj,1) = va(ji,jj,1) + (zhpj(ji,jj,1) + zvap) * vmask(ji,jj,1)644 puu(ji,jj,1,Krhs) = puu(ji,jj,1,Krhs) + (zhpi(ji,jj,1) + zuap) * umask(ji,jj,1) 645 pvv(ji,jj,1,Krhs) = pvv(ji,jj,1,Krhs) + (zhpj(ji,jj,1) + zvap) * vmask(ji,jj,1) 633 646 END DO 634 647 END DO … … 642 655 ! hydrostatic pressure gradient along s-surfaces 643 656 zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) + zcoef0 / e1u(ji,jj) & 644 & * ( e3w _n(ji+1,jj,jk) * ( rhd(ji+1,jj,jk) + rhd(ji+1,jj,jk-1) + 2*znad ) * wmask(ji+1,jj,jk) &645 & - e3w _n(ji ,jj,jk) * ( rhd(ji ,jj,jk) + rhd(ji ,jj,jk-1) + 2*znad ) * wmask(ji ,jj,jk) )657 & * ( e3w(ji+1,jj,jk,Kmm) * ( rhd(ji+1,jj,jk) + rhd(ji+1,jj,jk-1) + 2*znad ) * wmask(ji+1,jj,jk) & 658 & - e3w(ji ,jj,jk,Kmm) * ( rhd(ji ,jj,jk) + rhd(ji ,jj,jk-1) + 2*znad ) * wmask(ji ,jj,jk) ) 646 659 zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1) + zcoef0 / e2v(ji,jj) & 647 & * ( e3w _n(ji,jj+1,jk) * ( rhd(ji,jj+1,jk) + rhd(ji,jj+1,jk-1) + 2*znad ) * wmask(ji,jj+1,jk) &648 & - e3w _n(ji,jj ,jk) * ( rhd(ji,jj, jk) + rhd(ji,jj ,jk-1) + 2*znad ) * wmask(ji,jj ,jk) )660 & * ( e3w(ji,jj+1,jk,Kmm) * ( rhd(ji,jj+1,jk) + rhd(ji,jj+1,jk-1) + 2*znad ) * wmask(ji,jj+1,jk) & 661 & - e3w(ji,jj ,jk,Kmm) * ( rhd(ji,jj, jk) + rhd(ji,jj ,jk-1) + 2*znad ) * wmask(ji,jj ,jk) ) 649 662 ! s-coordinate pressure gradient correction 650 663 zuap = -zcoef0 * ( rhd (ji+1,jj ,jk) + rhd (ji,jj,jk) + 2._wp * znad ) & 651 & * ( gde3w _n(ji+1,jj ,jk) - gde3w_n(ji,jj,jk) ) / e1u(ji,jj)664 & * ( gde3w(ji+1,jj ,jk) - gde3w(ji,jj,jk) ) / e1u(ji,jj) 652 665 zvap = -zcoef0 * ( rhd (ji ,jj+1,jk) + rhd (ji,jj,jk) + 2._wp * znad ) & 653 & * ( gde3w _n(ji ,jj+1,jk) - gde3w_n(ji,jj,jk) ) / e2v(ji,jj)666 & * ( gde3w(ji ,jj+1,jk) - gde3w(ji,jj,jk) ) / e2v(ji,jj) 654 667 ! add to the general momentum trend 655 ua(ji,jj,jk) = ua(ji,jj,jk) + (zhpi(ji,jj,jk) + zuap) * umask(ji,jj,jk)656 va(ji,jj,jk) = va(ji,jj,jk) + (zhpj(ji,jj,jk) + zvap) * vmask(ji,jj,jk)668 puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) + (zhpi(ji,jj,jk) + zuap) * umask(ji,jj,jk) 669 pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) + (zhpj(ji,jj,jk) + zvap) * vmask(ji,jj,jk) 657 670 END DO 658 671 END DO … … 662 675 663 676 664 SUBROUTINE hpg_djc( kt )677 SUBROUTINE hpg_djc( kt, Kmm, puu, pvv, Krhs ) 665 678 !!--------------------------------------------------------------------- 666 679 !! *** ROUTINE hpg_djc *** … … 670 683 !! Reference: Shchepetkin and McWilliams, J. Geophys. Res., 108(C3), 3090, 2003 671 684 !!---------------------------------------------------------------------- 672 INTEGER, INTENT(in) :: kt ! ocean time-step index 685 INTEGER , INTENT( in ) :: kt ! ocean time-step index 686 INTEGER , INTENT( in ) :: Kmm, Krhs ! ocean time level indices 687 REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) :: puu, pvv ! ocean velocities and RHS of momentum equation 673 688 !! 674 689 INTEGER :: ji, jj, jk ! dummy loop indices … … 688 703 DO jj = 2, jpjm1 689 704 DO ji = 2, jpim1 690 ll_tmp1 = MIN( ssh n(ji,jj) , sshn(ji+1,jj) ) > &705 ll_tmp1 = MIN( ssh(ji,jj,Kmm) , ssh(ji+1,jj,Kmm) ) > & 691 706 & MAX( -ht_0(ji,jj) , -ht_0(ji+1,jj) ) .AND. & 692 & MAX( ssh n(ji,jj) + ht_0(ji,jj), sshn(ji+1,jj) + ht_0(ji+1,jj) ) &707 & MAX( ssh(ji,jj,Kmm) + ht_0(ji,jj), ssh(ji+1,jj,Kmm) + ht_0(ji+1,jj) ) & 693 708 & > rn_wdmin1 + rn_wdmin2 694 ll_tmp2 = ( ABS( ssh n(ji,jj) - sshn(ji+1,jj) ) > 1.E-12 ) .AND. ( &695 & MAX( ssh n(ji,jj) , sshn(ji+1,jj) ) > &709 ll_tmp2 = ( ABS( ssh(ji,jj,Kmm) - ssh(ji+1,jj,Kmm) ) > 1.E-12 ) .AND. ( & 710 & MAX( ssh(ji,jj,Kmm) , ssh(ji+1,jj,Kmm) ) > & 696 711 & MAX( -ht_0(ji,jj) , -ht_0(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 ) 697 712 IF(ll_tmp1) THEN 698 713 zcpx(ji,jj) = 1.0_wp 699 714 ELSE IF(ll_tmp2) THEN 700 ! no worries about ssh n(ji+1,jj) - sshn(ji ,jj) = 0, it won't happen ! here701 zcpx(ji,jj) = ABS( (ssh n(ji+1,jj) + ht_0(ji+1,jj) - sshn(ji,jj) - ht_0(ji,jj)) &702 & / (ssh n(ji+1,jj) - sshn(ji ,jj)) )715 ! no worries about ssh(ji+1,jj,Kmm) - ssh(ji ,jj,Kmm) = 0, it won't happen ! here 716 zcpx(ji,jj) = ABS( (ssh(ji+1,jj,Kmm) + ht_0(ji+1,jj) - ssh(ji,jj,Kmm) - ht_0(ji,jj)) & 717 & / (ssh(ji+1,jj,Kmm) - ssh(ji ,jj,Kmm)) ) 703 718 ELSE 704 719 zcpx(ji,jj) = 0._wp 705 720 END IF 706 721 707 ll_tmp1 = MIN( ssh n(ji,jj) , sshn(ji,jj+1) ) > &722 ll_tmp1 = MIN( ssh(ji,jj,Kmm) , ssh(ji,jj+1,Kmm) ) > & 708 723 & MAX( -ht_0(ji,jj) , -ht_0(ji,jj+1) ) .AND. & 709 & MAX( ssh n(ji,jj) + ht_0(ji,jj), sshn(ji,jj+1) + ht_0(ji,jj+1) ) &724 & MAX( ssh(ji,jj,Kmm) + ht_0(ji,jj), ssh(ji,jj+1,Kmm) + ht_0(ji,jj+1) ) & 710 725 & > rn_wdmin1 + rn_wdmin2 711 ll_tmp2 = ( ABS( ssh n(ji,jj) - sshn(ji,jj+1) ) > 1.E-12 ) .AND. ( &712 & MAX( ssh n(ji,jj) , sshn(ji,jj+1) ) > &726 ll_tmp2 = ( ABS( ssh(ji,jj,Kmm) - ssh(ji,jj+1,Kmm) ) > 1.E-12 ) .AND. ( & 727 & MAX( ssh(ji,jj,Kmm) , ssh(ji,jj+1,Kmm) ) > & 713 728 & MAX( -ht_0(ji,jj) , -ht_0(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 ) 714 729 … … 716 731 zcpy(ji,jj) = 1.0_wp 717 732 ELSE IF(ll_tmp2) THEN 718 ! no worries about ssh n(ji,jj+1) - sshn(ji,jj) = 0, it won't happen ! here719 zcpy(ji,jj) = ABS( (ssh n(ji,jj+1) + ht_0(ji,jj+1) - sshn(ji,jj) - ht_0(ji,jj)) &720 & / (ssh n(ji,jj+1) - sshn(ji,jj)) )733 ! no worries about ssh(ji,jj+1,Kmm) - ssh(ji,jj ,Kmm) = 0, it won't happen ! here 734 zcpy(ji,jj) = ABS( (ssh(ji,jj+1,Kmm) + ht_0(ji,jj+1) - ssh(ji,jj,Kmm) - ht_0(ji,jj)) & 735 & / (ssh(ji,jj+1,Kmm) - ssh(ji,jj ,Kmm)) ) 721 736 ELSE 722 737 zcpy(ji,jj) = 0._wp … … 748 763 DO ji = fs_2, fs_jpim1 ! vector opt. 749 764 drhoz(ji,jj,jk) = rhd (ji ,jj ,jk) - rhd (ji,jj,jk-1) 750 dzz (ji,jj,jk) = gde3w _n(ji ,jj ,jk) - gde3w_n(ji,jj,jk-1)765 dzz (ji,jj,jk) = gde3w(ji ,jj ,jk) - gde3w(ji,jj,jk-1) 751 766 drhox(ji,jj,jk) = rhd (ji+1,jj ,jk) - rhd (ji,jj,jk ) 752 dzx (ji,jj,jk) = gde3w _n(ji+1,jj ,jk) - gde3w_n(ji,jj,jk )767 dzx (ji,jj,jk) = gde3w(ji+1,jj ,jk) - gde3w(ji,jj,jk ) 753 768 drhoy(ji,jj,jk) = rhd (ji ,jj+1,jk) - rhd (ji,jj,jk ) 754 dzy (ji,jj,jk) = gde3w _n(ji ,jj+1,jk) - gde3w_n(ji,jj,jk )769 dzy (ji,jj,jk) = gde3w(ji ,jj+1,jk) - gde3w(ji,jj,jk ) 755 770 END DO 756 771 END DO … … 839 854 DO jj = 2, jpjm1 840 855 DO ji = fs_2, fs_jpim1 ! vector opt. 841 rho_k(ji,jj,1) = -grav * ( e3w _n(ji,jj,1) - gde3w_n(ji,jj,1) ) &856 rho_k(ji,jj,1) = -grav * ( e3w(ji,jj,1,Kmm) - gde3w(ji,jj,1) ) & 842 857 & * ( rhd(ji,jj,1) & 843 858 & + 0.5_wp * ( rhd (ji,jj,2) - rhd (ji,jj,1) ) & 844 & * ( e3w _n (ji,jj,1) - gde3w_n(ji,jj,1) ) &845 & / ( gde3w _n(ji,jj,2) - gde3w_n(ji,jj,1) ) )859 & * ( e3w (ji,jj,1,Kmm) - gde3w(ji,jj,1) ) & 860 & / ( gde3w(ji,jj,2) - gde3w(ji,jj,1) ) ) 846 861 END DO 847 862 END DO … … 855 870 856 871 rho_k(ji,jj,jk) = zcoef0 * ( rhd (ji,jj,jk) + rhd (ji,jj,jk-1) ) & 857 & * ( gde3w _n(ji,jj,jk) - gde3w_n(ji,jj,jk-1) ) &872 & * ( gde3w(ji,jj,jk) - gde3w(ji,jj,jk-1) ) & 858 873 & - grav * z1_10 * ( & 859 874 & ( drhow (ji,jj,jk) - drhow (ji,jj,jk-1) ) & 860 & * ( gde3w _n(ji,jj,jk) - gde3w_n(ji,jj,jk-1) - z1_12 * ( dzw (ji,jj,jk) + dzw (ji,jj,jk-1) ) ) &875 & * ( gde3w(ji,jj,jk) - gde3w(ji,jj,jk-1) - z1_12 * ( dzw (ji,jj,jk) + dzw (ji,jj,jk-1) ) ) & 861 876 & - ( dzw (ji,jj,jk) - dzw (ji,jj,jk-1) ) & 862 877 & * ( rhd (ji,jj,jk) - rhd (ji,jj,jk-1) - z1_12 * ( drhow(ji,jj,jk) + drhow(ji,jj,jk-1) ) ) & … … 864 879 865 880 rho_i(ji,jj,jk) = zcoef0 * ( rhd (ji+1,jj,jk) + rhd (ji,jj,jk) ) & 866 & * ( gde3w _n(ji+1,jj,jk) - gde3w_n(ji,jj,jk) ) &881 & * ( gde3w(ji+1,jj,jk) - gde3w(ji,jj,jk) ) & 867 882 & - grav* z1_10 * ( & 868 883 & ( drhou (ji+1,jj,jk) - drhou (ji,jj,jk) ) & 869 & * ( gde3w _n(ji+1,jj,jk) - gde3w_n(ji,jj,jk) - z1_12 * ( dzu (ji+1,jj,jk) + dzu (ji,jj,jk) ) ) &884 & * ( gde3w(ji+1,jj,jk) - gde3w(ji,jj,jk) - z1_12 * ( dzu (ji+1,jj,jk) + dzu (ji,jj,jk) ) ) & 870 885 & - ( dzu (ji+1,jj,jk) - dzu (ji,jj,jk) ) & 871 886 & * ( rhd (ji+1,jj,jk) - rhd (ji,jj,jk) - z1_12 * ( drhou(ji+1,jj,jk) + drhou(ji,jj,jk) ) ) & … … 873 888 874 889 rho_j(ji,jj,jk) = zcoef0 * ( rhd (ji,jj+1,jk) + rhd (ji,jj,jk) ) & 875 & * ( gde3w _n(ji,jj+1,jk) - gde3w_n(ji,jj,jk) ) &890 & * ( gde3w(ji,jj+1,jk) - gde3w(ji,jj,jk) ) & 876 891 & - grav* z1_10 * ( & 877 892 & ( drhov (ji,jj+1,jk) - drhov (ji,jj,jk) ) & 878 & * ( gde3w _n(ji,jj+1,jk) - gde3w_n(ji,jj,jk) - z1_12 * ( dzv (ji,jj+1,jk) + dzv (ji,jj,jk) ) ) &893 & * ( gde3w(ji,jj+1,jk) - gde3w(ji,jj,jk) - z1_12 * ( dzv (ji,jj+1,jk) + dzv (ji,jj,jk) ) ) & 879 894 & - ( dzv (ji,jj+1,jk) - dzv (ji,jj,jk) ) & 880 895 & * ( rhd (ji,jj+1,jk) - rhd (ji,jj,jk) - z1_12 * ( drhov(ji,jj+1,jk) + drhov(ji,jj,jk) ) ) & … … 898 913 ENDIF 899 914 ! add to the general momentum trend 900 ua(ji,jj,1) = ua(ji,jj,1) + zhpi(ji,jj,1)901 va(ji,jj,1) = va(ji,jj,1) + zhpj(ji,jj,1)915 puu(ji,jj,1,Krhs) = puu(ji,jj,1,Krhs) + zhpi(ji,jj,1) 916 pvv(ji,jj,1,Krhs) = pvv(ji,jj,1,Krhs) + zhpj(ji,jj,1) 902 917 END DO 903 918 END DO … … 921 936 ENDIF 922 937 ! add to the general momentum trend 923 ua(ji,jj,jk) = ua(ji,jj,jk) + zhpi(ji,jj,jk)924 va(ji,jj,jk) = va(ji,jj,jk) + zhpj(ji,jj,jk)938 puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) + zhpi(ji,jj,jk) 939 pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) + zhpj(ji,jj,jk) 925 940 END DO 926 941 END DO … … 932 947 933 948 934 SUBROUTINE hpg_prj( kt )949 SUBROUTINE hpg_prj( kt, Kmm, puu, pvv, Krhs ) 935 950 !!--------------------------------------------------------------------- 936 951 !! *** ROUTINE hpg_prj *** … … 941 956 !! all vertical coordinate systems 942 957 !! 943 !! ** Action : - Update ( ua,va) with the now hydrastatic pressure trend958 !! ** Action : - Update (puu(:,:,:,Krhs),pvv(:,:,:,Krhs)) with the now hydrastatic pressure trend 944 959 !!---------------------------------------------------------------------- 945 960 INTEGER, PARAMETER :: polynomial_type = 1 ! 1: cubic spline, 2: linear 946 INTEGER, INTENT(in) :: kt ! ocean time-step index 961 INTEGER , INTENT( in ) :: kt ! ocean time-step index 962 INTEGER , INTENT( in ) :: Kmm, Krhs ! ocean time level indices 963 REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) :: puu, pvv ! ocean velocities and RHS of momentum equation 947 964 !! 948 965 INTEGER :: ji, jj, jk, jkk ! dummy loop indices … … 976 993 DO jj = 2, jpjm1 977 994 DO ji = 2, jpim1 978 ll_tmp1 = MIN( ssh n(ji,jj) , sshn(ji+1,jj) ) > &995 ll_tmp1 = MIN( ssh(ji,jj,Kmm) , ssh(ji+1,jj,Kmm) ) > & 979 996 & MAX( -ht_0(ji,jj) , -ht_0(ji+1,jj) ) .AND. & 980 & MAX( ssh n(ji,jj) + ht_0(ji,jj), sshn(ji+1,jj) + ht_0(ji+1,jj) ) &997 & MAX( ssh(ji,jj,Kmm) + ht_0(ji,jj), ssh(ji+1,jj,Kmm) + ht_0(ji+1,jj) ) & 981 998 & > rn_wdmin1 + rn_wdmin2 982 ll_tmp2 = ( ABS( ssh n(ji,jj) - sshn(ji+1,jj) ) > 1.E-12 ) .AND. ( &983 & MAX( ssh n(ji,jj) , sshn(ji+1,jj) ) > &999 ll_tmp2 = ( ABS( ssh(ji,jj,Kmm) - ssh(ji+1,jj,Kmm) ) > 1.E-12 ) .AND. ( & 1000 & MAX( ssh(ji,jj,Kmm) , ssh(ji+1,jj,Kmm) ) > & 984 1001 & MAX( -ht_0(ji,jj) , -ht_0(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 ) 985 1002 … … 987 1004 zcpx(ji,jj) = 1.0_wp 988 1005 ELSE IF(ll_tmp2) THEN 989 ! no worries about ssh n(ji+1,jj) - sshn(ji ,jj) = 0, it won't happen ! here990 zcpx(ji,jj) = ABS( (ssh n(ji+1,jj) + ht_0(ji+1,jj) - sshn(ji,jj) - ht_0(ji,jj)) &991 & / (ssh n(ji+1,jj) - sshn(ji ,jj)) )1006 ! no worries about ssh(ji+1,jj,Kmm) - ssh(ji ,jj,Kmm) = 0, it won't happen ! here 1007 zcpx(ji,jj) = ABS( (ssh(ji+1,jj,Kmm) + ht_0(ji+1,jj) - ssh(ji,jj,Kmm) - ht_0(ji,jj)) & 1008 & / (ssh(ji+1,jj,Kmm) - ssh(ji ,jj,Kmm)) ) 992 1009 993 1010 zcpx(ji,jj) = max(min( zcpx(ji,jj) , 1.0_wp),0.0_wp) … … 996 1013 END IF 997 1014 998 ll_tmp1 = MIN( ssh n(ji,jj) , sshn(ji,jj+1) ) > &1015 ll_tmp1 = MIN( ssh(ji,jj,Kmm) , ssh(ji,jj+1,Kmm) ) > & 999 1016 & MAX( -ht_0(ji,jj) , -ht_0(ji,jj+1) ) .AND. & 1000 & MAX( ssh n(ji,jj) + ht_0(ji,jj), sshn(ji,jj+1) + ht_0(ji,jj+1) ) &1017 & MAX( ssh(ji,jj,Kmm) + ht_0(ji,jj), ssh(ji,jj+1,Kmm) + ht_0(ji,jj+1) ) & 1001 1018 & > rn_wdmin1 + rn_wdmin2 1002 ll_tmp2 = ( ABS( ssh n(ji,jj) - sshn(ji,jj+1) ) > 1.E-12 ) .AND. ( &1003 & MAX( ssh n(ji,jj) , sshn(ji,jj+1) ) > &1019 ll_tmp2 = ( ABS( ssh(ji,jj,Kmm) - ssh(ji,jj+1,Kmm) ) > 1.E-12 ) .AND. ( & 1020 & MAX( ssh(ji,jj,Kmm) , ssh(ji,jj+1,Kmm) ) > & 1004 1021 & MAX( -ht_0(ji,jj) , -ht_0(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 ) 1005 1022 … … 1007 1024 zcpy(ji,jj) = 1.0_wp 1008 1025 ELSE IF(ll_tmp2) THEN 1009 ! no worries about ssh n(ji,jj+1) - sshn(ji,jj) = 0, it won't happen ! here1010 zcpy(ji,jj) = ABS( (ssh n(ji,jj+1) + ht_0(ji,jj+1) - sshn(ji,jj) - ht_0(ji,jj)) &1011 & / (ssh n(ji,jj+1) - sshn(ji,jj)) )1026 ! no worries about ssh(ji,jj+1,Kmm) - ssh(ji,jj ,Kmm) = 0, it won't happen ! here 1027 zcpy(ji,jj) = ABS( (ssh(ji,jj+1,Kmm) + ht_0(ji,jj+1) - ssh(ji,jj,Kmm) - ht_0(ji,jj)) & 1028 & / (ssh(ji,jj+1,Kmm) - ssh(ji,jj ,Kmm)) ) 1012 1029 zcpy(ji,jj) = max(min( zcpy(ji,jj) , 1.0_wp),0.0_wp) 1013 1030 … … 1032 1049 ELSEIF( jk < jpkm1 ) THEN 1033 1050 DO jkk = jk+1, jpk 1034 zrhh(ji,jj,jkk) = interp1(gde3w _n(ji,jj,jkk ), gde3w_n(ji,jj,jkk-1), &1035 & gde3w _n(ji,jj,jkk-2), rhd (ji,jj,jkk-1), rhd(ji,jj,jkk-2))1051 zrhh(ji,jj,jkk) = interp1(gde3w(ji,jj,jkk ), gde3w(ji,jj,jkk-1), & 1052 & gde3w(ji,jj,jkk-2), rhd (ji,jj,jkk-1), rhd(ji,jj,jkk-2)) 1036 1053 END DO 1037 1054 ENDIF … … 1042 1059 DO jj = 1, jpj 1043 1060 DO ji = 1, jpi 1044 zdept(ji,jj,1) = 0.5_wp * e3w _n(ji,jj,1) - sshn(ji,jj) * znad1061 zdept(ji,jj,1) = 0.5_wp * e3w(ji,jj,1,Kmm) - ssh(ji,jj,Kmm) * znad 1045 1062 END DO 1046 1063 END DO … … 1049 1066 DO jj = 1, jpj 1050 1067 DO ji = 1, jpi 1051 zdept(ji,jj,jk) = zdept(ji,jj,jk-1) + e3w _n(ji,jj,jk)1068 zdept(ji,jj,jk) = zdept(ji,jj,jk-1) + e3w(ji,jj,jk,Kmm) 1052 1069 END DO 1053 1070 END DO … … 1066 1083 DO ji = 2, jpi 1067 1084 zrhdt1 = zrhh(ji,jj,1) - interp3( zdept(ji,jj,1), asp(ji,jj,1), bsp(ji,jj,1), & 1068 & csp(ji,jj,1), dsp(ji,jj,1) ) * 0.25_wp * e3w _n(ji,jj,1)1085 & csp(ji,jj,1), dsp(ji,jj,1) ) * 0.25_wp * e3w(ji,jj,1,Kmm) 1069 1086 1070 1087 ! assuming linear profile across the top half surface layer 1071 zhpi(ji,jj,1) = 0.5_wp * e3w _n(ji,jj,1) * zrhdt11088 zhpi(ji,jj,1) = 0.5_wp * e3w(ji,jj,1,Kmm) * zrhdt1 1072 1089 END DO 1073 1090 END DO … … 1091 1108 DO ji = 2, jpim1 1092 1109 !!gm BUG ? if it is ssh at u- & v-point then it should be: 1093 ! zsshu_n(ji,jj) = (e1e2t(ji,jj) * ssh n(ji,jj) + e1e2t(ji+1,jj) * sshn(ji+1,jj)) * &1110 ! zsshu_n(ji,jj) = (e1e2t(ji,jj) * ssh(ji,jj,Kmm) + e1e2t(ji+1,jj) * ssh(ji+1,jj,Kmm)) * & 1094 1111 ! & r1_e1e2u(ji,jj) * umask(ji,jj,1) * 0.5_wp 1095 ! zsshv_n(ji,jj) = (e1e2t(ji,jj) * ssh n(ji,jj) + e1e2t(ji,jj+1) * sshn(ji,jj+1)) * &1112 ! zsshv_n(ji,jj) = (e1e2t(ji,jj) * ssh(ji,jj,Kmm) + e1e2t(ji,jj+1) * ssh(ji,jj+1,Kmm)) * & 1096 1113 ! & r1_e1e2v(ji,jj) * vmask(ji,jj,1) * 0.5_wp 1097 1114 !!gm not this: 1098 zsshu_n(ji,jj) = (e1e2u(ji,jj) * ssh n(ji,jj) + e1e2u(ji+1, jj) * sshn(ji+1,jj)) * &1115 zsshu_n(ji,jj) = (e1e2u(ji,jj) * ssh(ji,jj,Kmm) + e1e2u(ji+1, jj) * ssh(ji+1,jj,Kmm)) * & 1099 1116 & r1_e1e2u(ji,jj) * umask(ji,jj,1) * 0.5_wp 1100 zsshv_n(ji,jj) = (e1e2v(ji,jj) * ssh n(ji,jj) + e1e2v(ji+1, jj) * sshn(ji,jj+1)) * &1117 zsshv_n(ji,jj) = (e1e2v(ji,jj) * ssh(ji,jj,Kmm) + e1e2v(ji+1, jj) * ssh(ji,jj+1,Kmm)) * & 1101 1118 & r1_e1e2v(ji,jj) * vmask(ji,jj,1) * 0.5_wp 1102 1119 END DO … … 1107 1124 DO jj = 2, jpjm1 1108 1125 DO ji = 2, jpim1 1109 zu(ji,jj,1) = - ( e3u _n(ji,jj,1) - zsshu_n(ji,jj) * znad)1110 zv(ji,jj,1) = - ( e3v _n(ji,jj,1) - zsshv_n(ji,jj) * znad)1126 zu(ji,jj,1) = - ( e3u(ji,jj,1,Kmm) - zsshu_n(ji,jj) * znad) 1127 zv(ji,jj,1) = - ( e3v(ji,jj,1,Kmm) - zsshv_n(ji,jj) * znad) 1111 1128 END DO 1112 1129 END DO … … 1115 1132 DO jj = 2, jpjm1 1116 1133 DO ji = 2, jpim1 1117 zu(ji,jj,jk) = zu(ji,jj,jk-1) - e3u _n(ji,jj,jk)1118 zv(ji,jj,jk) = zv(ji,jj,jk-1) - e3v _n(ji,jj,jk)1134 zu(ji,jj,jk) = zu(ji,jj,jk-1) - e3u(ji,jj,jk,Kmm) 1135 zv(ji,jj,jk) = zv(ji,jj,jk-1) - e3v(ji,jj,jk,Kmm) 1119 1136 END DO 1120 1137 END DO … … 1124 1141 DO jj = 2, jpjm1 1125 1142 DO ji = 2, jpim1 1126 zu(ji,jj,jk) = zu(ji,jj,jk) + 0.5_wp * e3u _n(ji,jj,jk)1127 zv(ji,jj,jk) = zv(ji,jj,jk) + 0.5_wp * e3v _n(ji,jj,jk)1143 zu(ji,jj,jk) = zu(ji,jj,jk) + 0.5_wp * e3u(ji,jj,jk,Kmm) 1144 zv(ji,jj,jk) = zv(ji,jj,jk) + 0.5_wp * e3v(ji,jj,jk,Kmm) 1128 1145 END DO 1129 1146 END DO … … 1177 1194 DO WHILE ( -zdept(jid,jj,jk1) < zuijk ) 1178 1195 IF( jk1 == 1 ) THEN 1179 zdeps = zdept(jid,jj,1) + MIN(zuijk, ssh n(jid,jj)*znad)1196 zdeps = zdept(jid,jj,1) + MIN(zuijk, ssh(jid,jj,Kmm)*znad) 1180 1197 zrhdt1 = zrhh(jid,jj,1) - interp3(zdept(jid,jj,1), asp(jid,jj,1), & 1181 1198 bsp(jid,jj,1), csp(jid,jj,1), & … … 1197 1214 IF( .NOT.ln_linssh ) THEN 1198 1215 zdpdx2 = zcoef0 * r1_e1u(ji,jj) * & 1199 & ( REAL(jis-jid, wp) * (zpwes + zpwed) + (ssh n(ji+1,jj)-sshn(ji,jj)) )1216 & ( REAL(jis-jid, wp) * (zpwes + zpwed) + (ssh(ji+1,jj,Kmm)-ssh(ji,jj,Kmm)) ) 1200 1217 ELSE 1201 1218 zdpdx2 = zcoef0 * r1_e1u(ji,jj) * REAL(jis-jid, wp) * (zpwes + zpwed) … … 1205 1222 zdpdx2 = zdpdx2 * zcpx(ji,jj) * wdrampu(ji,jj) 1206 1223 ENDIF 1207 ua(ji,jj,jk) = ua(ji,jj,jk) + (zdpdx1 + zdpdx2) * umask(ji,jj,jk)1224 puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) + (zdpdx1 + zdpdx2) * umask(ji,jj,jk) 1208 1225 ENDIF 1209 1226 … … 1235 1252 DO WHILE ( -zdept(ji,jjd,jk1) < zvijk ) 1236 1253 IF( jk1 == 1 ) THEN 1237 zdeps = zdept(ji,jjd,1) + MIN(zvijk, ssh n(ji,jjd)*znad)1254 zdeps = zdept(ji,jjd,1) + MIN(zvijk, ssh(ji,jjd,Kmm)*znad) 1238 1255 zrhdt1 = zrhh(ji,jjd,1) - interp3(zdept(ji,jjd,1), asp(ji,jjd,1), & 1239 1256 bsp(ji,jjd,1), csp(ji,jjd,1), & … … 1256 1273 IF( .NOT.ln_linssh ) THEN 1257 1274 zdpdy2 = zcoef0 * r1_e2v(ji,jj) * & 1258 ( REAL(jjs-jjd, wp) * (zpnss + zpnsd) + (ssh n(ji,jj+1)-sshn(ji,jj)) )1275 ( REAL(jjs-jjd, wp) * (zpnss + zpnsd) + (ssh(ji,jj+1,Kmm)-ssh(ji,jj,Kmm)) ) 1259 1276 ELSE 1260 1277 zdpdy2 = zcoef0 * r1_e2v(ji,jj) * REAL(jjs-jjd, wp) * (zpnss + zpnsd ) … … 1265 1282 ENDIF 1266 1283 1267 va(ji,jj,jk) = va(ji,jj,jk) + (zdpdy1 + zdpdy2) * vmask(ji,jj,jk)1284 pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) + (zdpdy1 + zdpdy2) * vmask(ji,jj,jk) 1268 1285 ENDIF 1269 1286 !
Note: See TracChangeset
for help on using the changeset viewer.
