Changeset 10928 for NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN/dynhpg.F90
- Timestamp:
- 2019-05-03T17:44:56+02:00 (5 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN/dynhpg.F90
r10491 r10928 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(:,:,:)118 ztrdu(:,:,:) = puu(:,:,:,Krhs) - ztrdu(:,:,:) 119 ztrdv(:,:,:) = pvv(:,:,:,Krhs) - ztrdv(:,:,:) 117 120 CALL trd_dyn( ztrdu, ztrdv, jpdyn_hpg, kt ) 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 … … 347 356 348 357 ! Partial steps: bottom before horizontal gradient of t, s, rd at the last ocean level 349 !jc CALL zps_hde ( kt, jpts, ts n, gtsu, gtsv, rhd, gru , grv )358 !jc CALL zps_hde ( kt, jpts, ts(:,:,:,:,Kmm), gtsu, gtsv, rhd, gru , grv ) 350 359 351 360 ! Local constant initialization … … 355 364 DO jj = 2, jpjm1 356 365 DO ji = fs_2, fs_jpim1 ! vector opt. 357 zcoef1 = zcoef0 * e3w _n(ji,jj,1)366 zcoef1 = zcoef0 * e3w(ji,jj,1,Kmm) 358 367 ! hydrostatic pressure gradient 359 368 zhpi(ji,jj,1) = zcoef1 * ( rhd(ji+1,jj ,1) - rhd(ji,jj,1) ) * r1_e1u(ji,jj) 360 369 zhpj(ji,jj,1) = zcoef1 * ( rhd(ji ,jj+1,1) - rhd(ji,jj,1) ) * r1_e2v(ji,jj) 361 370 ! add to the general momentum trend 362 ua(ji,jj,1) = ua(ji,jj,1) + zhpi(ji,jj,1)363 va(ji,jj,1) = va(ji,jj,1) + zhpj(ji,jj,1)371 puu(ji,jj,1,Krhs) = puu(ji,jj,1,Krhs) + zhpi(ji,jj,1) 372 pvv(ji,jj,1,Krhs) = pvv(ji,jj,1,Krhs) + zhpj(ji,jj,1) 364 373 END DO 365 374 END DO … … 369 378 DO jj = 2, jpjm1 370 379 DO ji = fs_2, fs_jpim1 ! vector opt. 371 zcoef1 = zcoef0 * e3w _n(ji,jj,jk)380 zcoef1 = zcoef0 * e3w(ji,jj,jk,Kmm) 372 381 ! hydrostatic pressure gradient 373 382 zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) & … … 379 388 & - ( rhd(ji,jj, jk) + rhd(ji,jj ,jk-1) ) ) * r1_e2v(ji,jj) 380 389 ! add to the general momentum trend 381 ua(ji,jj,jk) = ua(ji,jj,jk) + zhpi(ji,jj,jk)382 va(ji,jj,jk) = va(ji,jj,jk) + zhpj(ji,jj,jk)390 puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) + zhpi(ji,jj,jk) 391 pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) + zhpj(ji,jj,jk) 383 392 END DO 384 393 END DO … … 390 399 iku = mbku(ji,jj) 391 400 ikv = mbkv(ji,jj) 392 zcoef2 = zcoef0 * MIN( e3w _n(ji,jj,iku), e3w_n(ji+1,jj ,iku) )393 zcoef3 = zcoef0 * MIN( e3w _n(ji,jj,ikv), e3w_n(ji ,jj+1,ikv) )401 zcoef2 = zcoef0 * MIN( e3w(ji,jj,iku,Kmm), e3w(ji+1,jj ,iku,Kmm) ) 402 zcoef3 = zcoef0 * MIN( e3w(ji,jj,ikv,Kmm), e3w(ji ,jj+1,ikv,Kmm) ) 394 403 IF( iku > 1 ) THEN ! on i-direction (level 2 or more) 395 ua (ji,jj,iku) = ua(ji,jj,iku) - zhpi(ji,jj,iku) ! subtract old value404 puu (ji,jj,iku,Krhs) = puu(ji,jj,iku,Krhs) - zhpi(ji,jj,iku) ! subtract old value 396 405 zhpi(ji,jj,iku) = zhpi(ji,jj,iku-1) & ! compute the new one 397 406 & + zcoef2 * ( rhd(ji+1,jj,iku-1) - rhd(ji,jj,iku-1) + gru(ji,jj) ) * r1_e1u(ji,jj) 398 ua (ji,jj,iku) = ua(ji,jj,iku) + zhpi(ji,jj,iku) ! add the new one to the general momentum trend407 puu (ji,jj,iku,Krhs) = puu(ji,jj,iku,Krhs) + zhpi(ji,jj,iku) ! add the new one to the general momentum trend 399 408 ENDIF 400 409 IF( ikv > 1 ) THEN ! on j-direction (level 2 or more) 401 va (ji,jj,ikv) = va(ji,jj,ikv) - zhpj(ji,jj,ikv) ! subtract old value410 pvv (ji,jj,ikv,Krhs) = pvv(ji,jj,ikv,Krhs) - zhpj(ji,jj,ikv) ! subtract old value 402 411 zhpj(ji,jj,ikv) = zhpj(ji,jj,ikv-1) & ! compute the new one 403 412 & + zcoef3 * ( rhd(ji,jj+1,ikv-1) - rhd(ji,jj,ikv-1) + grv(ji,jj) ) * r1_e2v(ji,jj) 404 va (ji,jj,ikv) = va(ji,jj,ikv) + zhpj(ji,jj,ikv) ! add the new one to the general momentum trend413 pvv (ji,jj,ikv,Krhs) = pvv(ji,jj,ikv,Krhs) + zhpj(ji,jj,ikv) ! add the new one to the general momentum trend 405 414 ENDIF 406 415 END DO … … 410 419 411 420 412 SUBROUTINE hpg_sco( kt )421 SUBROUTINE hpg_sco( kt, Kmm, puu, pvv, Krhs ) 413 422 !!--------------------------------------------------------------------- 414 423 !! *** ROUTINE hpg_sco *** … … 422 431 !! zhpi = grav ..... + 1/e1u mi(rhd) di[ grav dep3w ] 423 432 !! zhpj = grav ..... + 1/e2v mj(rhd) dj[ grav dep3w ] 424 !! add it to the general momentum trend (ua,va). 425 !! ua = ua - 1/e1u * zhpi 426 !! va = va - 1/e2v * zhpj 427 !! 428 !! ** Action : - Update (ua,va) with the now hydrastatic pressure trend 429 !!---------------------------------------------------------------------- 430 INTEGER, INTENT(in) :: kt ! ocean time-step index 433 !! add it to the general momentum trend (puu(:,:,:,Krhs),pvv(:,:,:,Krhs)). 434 !! puu(:,:,:,Krhs) = puu(:,:,:,Krhs) - 1/e1u * zhpi 435 !! pvv(:,:,:,Krhs) = pvv(:,:,:,Krhs) - 1/e2v * zhpj 436 !! 437 !! ** Action : - Update (puu(:,:,:,Krhs),pvv(:,:,:,Krhs)) with the now hydrastatic pressure trend 438 !!---------------------------------------------------------------------- 439 INTEGER , INTENT( in ) :: kt ! ocean time-step index 440 INTEGER , INTENT( in ) :: Kmm, Krhs ! ocean time level indices 441 REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) :: puu, pvv ! ocean velocities and RHS of momentum equation 431 442 !! 432 443 INTEGER :: ji, jj, jk, jii, jjj ! dummy loop indices … … 453 464 DO jj = 2, jpjm1 454 465 DO ji = 2, jpim1 455 ll_tmp1 = MIN( ssh n(ji,jj) , sshn(ji+1,jj) ) > &466 ll_tmp1 = MIN( ssh(ji,jj,Kmm) , ssh(ji+1,jj,Kmm) ) > & 456 467 & MAX( -ht_0(ji,jj) , -ht_0(ji+1,jj) ) .AND. & 457 & MAX( ssh n(ji,jj) + ht_0(ji,jj), sshn(ji+1,jj) + ht_0(ji+1,jj) ) &468 & MAX( ssh(ji,jj,Kmm) + ht_0(ji,jj), ssh(ji+1,jj,Kmm) + ht_0(ji+1,jj) ) & 458 469 & > rn_wdmin1 + rn_wdmin2 459 ll_tmp2 = ( ABS( ssh n(ji,jj) - sshn(ji+1,jj) ) > 1.E-12 ) .AND. ( &460 & MAX( ssh n(ji,jj) , sshn(ji+1,jj) ) > &470 ll_tmp2 = ( ABS( ssh(ji,jj,Kmm) - ssh(ji+1,jj,Kmm) ) > 1.E-12 ) .AND. ( & 471 & MAX( ssh(ji,jj,Kmm) , ssh(ji+1,jj,Kmm) ) > & 461 472 & MAX( -ht_0(ji,jj) , -ht_0(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 ) 462 473 … … 464 475 zcpx(ji,jj) = 1.0_wp 465 476 ELSE IF(ll_tmp2) THEN 466 ! no worries about ssh n(ji+1,jj) - sshn(ji ,jj) = 0, it won't happen ! here467 zcpx(ji,jj) = ABS( (ssh n(ji+1,jj) + ht_0(ji+1,jj) - sshn(ji,jj) - ht_0(ji,jj)) &468 & / (ssh n(ji+1,jj) - sshn(ji ,jj)) )477 ! no worries about ssh(ji+1,jj,Kmm) - ssh(ji ,jj,Kmm) = 0, it won't happen ! here 478 zcpx(ji,jj) = ABS( (ssh(ji+1,jj,Kmm) + ht_0(ji+1,jj) - ssh(ji,jj,Kmm) - ht_0(ji,jj)) & 479 & / (ssh(ji+1,jj,Kmm) - ssh(ji ,jj,Kmm)) ) 469 480 ELSE 470 481 zcpx(ji,jj) = 0._wp 471 482 END IF 472 483 473 ll_tmp1 = MIN( ssh n(ji,jj) , sshn(ji,jj+1) ) > &484 ll_tmp1 = MIN( ssh(ji,jj,Kmm) , ssh(ji,jj+1,Kmm) ) > & 474 485 & MAX( -ht_0(ji,jj) , -ht_0(ji,jj+1) ) .AND. & 475 & MAX( ssh n(ji,jj) + ht_0(ji,jj), sshn(ji,jj+1) + ht_0(ji,jj+1) ) &486 & MAX( ssh(ji,jj,Kmm) + ht_0(ji,jj), ssh(ji,jj+1,Kmm) + ht_0(ji,jj+1) ) & 476 487 & > rn_wdmin1 + rn_wdmin2 477 ll_tmp2 = ( ABS( ssh n(ji,jj) - sshn(ji,jj+1) ) > 1.E-12 ) .AND. ( &478 & MAX( ssh n(ji,jj) , sshn(ji,jj+1) ) > &488 ll_tmp2 = ( ABS( ssh(ji,jj,Kmm) - ssh(ji,jj+1,Kmm) ) > 1.E-12 ) .AND. ( & 489 & MAX( ssh(ji,jj,Kmm) , ssh(ji,jj+1,Kmm) ) > & 479 490 & MAX( -ht_0(ji,jj) , -ht_0(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 ) 480 491 … … 482 493 zcpy(ji,jj) = 1.0_wp 483 494 ELSE IF(ll_tmp2) THEN 484 ! no worries about ssh n(ji,jj+1) - sshn(ji,jj) = 0, it won't happen ! here485 zcpy(ji,jj) = ABS( (ssh n(ji,jj+1) + ht_0(ji,jj+1) - sshn(ji,jj) - ht_0(ji,jj)) &486 & / (ssh n(ji,jj+1) - sshn(ji,jj)) )495 ! no worries about ssh(ji,jj+1,Kmm) - ssh(ji,jj ,Kmm) = 0, it won't happen ! here 496 zcpy(ji,jj) = ABS( (ssh(ji,jj+1,Kmm) + ht_0(ji,jj+1) - ssh(ji,jj,Kmm) - ht_0(ji,jj)) & 497 & / (ssh(ji,jj+1,Kmm) - ssh(ji,jj ,Kmm)) ) 487 498 ELSE 488 499 zcpy(ji,jj) = 0._wp … … 497 508 DO ji = fs_2, fs_jpim1 ! vector opt. 498 509 ! hydrostatic pressure gradient along s-surfaces 499 zhpi(ji,jj,1) = zcoef0 * ( e3w _n(ji+1,jj ,1) * ( znad + rhd(ji+1,jj ,1) ) &500 & - e3w _n(ji ,jj ,1) * ( znad + rhd(ji ,jj ,1) ) ) * r1_e1u(ji,jj)501 zhpj(ji,jj,1) = zcoef0 * ( e3w _n(ji ,jj+1,1) * ( znad + rhd(ji ,jj+1,1) ) &502 & - e3w _n(ji ,jj ,1) * ( znad + rhd(ji ,jj ,1) ) ) * r1_e2v(ji,jj)510 zhpi(ji,jj,1) = zcoef0 * ( e3w(ji+1,jj ,1,Kmm) * ( znad + rhd(ji+1,jj ,1) ) & 511 & - e3w(ji ,jj ,1,Kmm) * ( znad + rhd(ji ,jj ,1) ) ) * r1_e1u(ji,jj) 512 zhpj(ji,jj,1) = zcoef0 * ( e3w(ji ,jj+1,1,Kmm) * ( znad + rhd(ji ,jj+1,1) ) & 513 & - e3w(ji ,jj ,1,Kmm) * ( znad + rhd(ji ,jj ,1) ) ) * r1_e2v(ji,jj) 503 514 ! s-coordinate pressure gradient correction 504 515 zuap = -zcoef0 * ( rhd (ji+1,jj,1) + rhd (ji,jj,1) + 2._wp * znad ) & 505 & * ( gde3w _n(ji+1,jj,1) - gde3w_n(ji,jj,1) ) * r1_e1u(ji,jj)516 & * ( gde3w(ji+1,jj,1) - gde3w(ji,jj,1) ) * r1_e1u(ji,jj) 506 517 zvap = -zcoef0 * ( rhd (ji,jj+1,1) + rhd (ji,jj,1) + 2._wp * znad ) & 507 & * ( gde3w _n(ji,jj+1,1) - gde3w_n(ji,jj,1) ) * r1_e2v(ji,jj)518 & * ( gde3w(ji,jj+1,1) - gde3w(ji,jj,1) ) * r1_e2v(ji,jj) 508 519 ! 509 520 IF( ln_wd_il ) THEN … … 515 526 ! 516 527 ! add to the general momentum trend 517 ua(ji,jj,1) = ua(ji,jj,1) + zhpi(ji,jj,1) + zuap518 va(ji,jj,1) = va(ji,jj,1) + zhpj(ji,jj,1) + zvap528 puu(ji,jj,1,Krhs) = puu(ji,jj,1,Krhs) + zhpi(ji,jj,1) + zuap 529 pvv(ji,jj,1,Krhs) = pvv(ji,jj,1,Krhs) + zhpj(ji,jj,1) + zvap 519 530 END DO 520 531 END DO … … 526 537 ! hydrostatic pressure gradient along s-surfaces 527 538 zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) + zcoef0 * r1_e1u(ji,jj) & 528 & * ( e3w _n(ji+1,jj,jk) * ( rhd(ji+1,jj,jk) + rhd(ji+1,jj,jk-1) + 2*znad ) &529 & - e3w _n(ji ,jj,jk) * ( rhd(ji ,jj,jk) + rhd(ji ,jj,jk-1) + 2*znad ) )539 & * ( e3w(ji+1,jj,jk,Kmm) * ( rhd(ji+1,jj,jk) + rhd(ji+1,jj,jk-1) + 2*znad ) & 540 & - e3w(ji ,jj,jk,Kmm) * ( rhd(ji ,jj,jk) + rhd(ji ,jj,jk-1) + 2*znad ) ) 530 541 zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1) + zcoef0 * r1_e2v(ji,jj) & 531 & * ( e3w _n(ji,jj+1,jk) * ( rhd(ji,jj+1,jk) + rhd(ji,jj+1,jk-1) + 2*znad ) &532 & - e3w _n(ji,jj ,jk) * ( rhd(ji,jj, jk) + rhd(ji,jj ,jk-1) + 2*znad ) )542 & * ( e3w(ji,jj+1,jk,Kmm) * ( rhd(ji,jj+1,jk) + rhd(ji,jj+1,jk-1) + 2*znad ) & 543 & - e3w(ji,jj ,jk,Kmm) * ( rhd(ji,jj, jk) + rhd(ji,jj ,jk-1) + 2*znad ) ) 533 544 ! s-coordinate pressure gradient correction 534 545 zuap = -zcoef0 * ( rhd (ji+1,jj ,jk) + rhd (ji,jj,jk) + 2._wp * znad ) & 535 & * ( gde3w _n(ji+1,jj ,jk) - gde3w_n(ji,jj,jk) ) * r1_e1u(ji,jj)546 & * ( gde3w(ji+1,jj ,jk) - gde3w(ji,jj,jk) ) * r1_e1u(ji,jj) 536 547 zvap = -zcoef0 * ( rhd (ji ,jj+1,jk) + rhd (ji,jj,jk) + 2._wp * znad ) & 537 & * ( gde3w _n(ji ,jj+1,jk) - gde3w_n(ji,jj,jk) ) * r1_e2v(ji,jj)548 & * ( gde3w(ji ,jj+1,jk) - gde3w(ji,jj,jk) ) * r1_e2v(ji,jj) 538 549 ! 539 550 IF( ln_wd_il ) THEN … … 545 556 ! 546 557 ! add to the general momentum trend 547 ua(ji,jj,jk) = ua(ji,jj,jk) + zhpi(ji,jj,jk) + zuap548 va(ji,jj,jk) = va(ji,jj,jk) + zhpj(ji,jj,jk) + zvap558 puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) + zhpi(ji,jj,jk) + zuap 559 pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) + zhpj(ji,jj,jk) + zvap 549 560 END DO 550 561 END DO … … 556 567 557 568 558 SUBROUTINE hpg_isf( kt )569 SUBROUTINE hpg_isf( kt, Kmm, puu, pvv, Krhs ) 559 570 !!--------------------------------------------------------------------- 560 571 !! *** ROUTINE hpg_isf *** … … 568 579 !! zhpi = grav ..... + 1/e1u mi(rhd) di[ grav dep3w ] 569 580 !! zhpj = grav ..... + 1/e2v mj(rhd) dj[ grav dep3w ] 570 !! add it to the general momentum trend ( ua,va).571 !! ua = ua- 1/e1u * zhpi572 !! va = va- 1/e2v * zhpj581 !! add it to the general momentum trend (puu(:,:,:,Krhs),pvv(:,:,:,Krhs)). 582 !! puu(:,:,:,Krhs) = puu(:,:,:,Krhs) - 1/e1u * zhpi 583 !! pvv(:,:,:,Krhs) = pvv(:,:,:,Krhs) - 1/e2v * zhpj 573 584 !! iceload is added and partial cell case are added to the top and bottom 574 585 !! 575 !! ** Action : - Update (ua,va) with the now hydrastatic pressure trend 576 !!---------------------------------------------------------------------- 577 INTEGER, INTENT(in) :: kt ! ocean time-step index 586 !! ** Action : - Update (puu(:,:,:,Krhs),pvv(:,:,:,Krhs)) with the now hydrastatic pressure trend 587 !!---------------------------------------------------------------------- 588 INTEGER , INTENT( in ) :: kt ! ocean time-step index 589 INTEGER , INTENT( in ) :: Kmm, Krhs ! ocean time level indices 590 REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) :: puu, pvv ! ocean velocities and RHS of momentum equation 578 591 !! 579 592 INTEGER :: ji, jj, jk, ikt, iktp1i, iktp1j ! dummy loop indices … … 596 609 DO jj = 1, jpj 597 610 ikt = mikt(ji,jj) 598 zts_top(ji,jj,1) = ts n(ji,jj,ikt,1)599 zts_top(ji,jj,2) = ts n(ji,jj,ikt,2)611 zts_top(ji,jj,1) = ts(ji,jj,ikt,1,Kmm) 612 zts_top(ji,jj,2) = ts(ji,jj,ikt,2,Kmm) 600 613 END DO 601 614 END DO … … 612 625 ! hydrostatic pressure gradient along s-surfaces and ice shelf pressure 613 626 ! we assume ISF is in isostatic equilibrium 614 zhpi(ji,jj,1) = zcoef0 / e1u(ji,jj) * ( 0.5_wp * e3w _n(ji+1,jj,iktp1i) &627 zhpi(ji,jj,1) = zcoef0 / e1u(ji,jj) * ( 0.5_wp * e3w(ji+1,jj,iktp1i,Kmm) & 615 628 & * ( 2._wp * znad + rhd(ji+1,jj,iktp1i) + zrhdtop_oce(ji+1,jj) ) & 616 & - 0.5_wp * e3w _n(ji,jj,ikt) &629 & - 0.5_wp * e3w(ji,jj,ikt,Kmm) & 617 630 & * ( 2._wp * znad + rhd(ji,jj,ikt) + zrhdtop_oce(ji,jj) ) & 618 631 & + ( riceload(ji+1,jj) - riceload(ji,jj)) ) 619 zhpj(ji,jj,1) = zcoef0 / e2v(ji,jj) * ( 0.5_wp * e3w _n(ji,jj+1,iktp1j) &632 zhpj(ji,jj,1) = zcoef0 / e2v(ji,jj) * ( 0.5_wp * e3w(ji,jj+1,iktp1j,Kmm) & 620 633 & * ( 2._wp * znad + rhd(ji,jj+1,iktp1j) + zrhdtop_oce(ji,jj+1) ) & 621 & - 0.5_wp * e3w _n(ji,jj,ikt) &634 & - 0.5_wp * e3w(ji,jj,ikt,Kmm) & 622 635 & * ( 2._wp * znad + rhd(ji,jj,ikt) + zrhdtop_oce(ji,jj) ) & 623 636 & + ( riceload(ji,jj+1) - riceload(ji,jj)) ) 624 637 ! s-coordinate pressure gradient correction (=0 if z coordinate) 625 638 zuap = -zcoef0 * ( rhd (ji+1,jj,1) + rhd (ji,jj,1) + 2._wp * znad ) & 626 & * ( gde3w _n(ji+1,jj,1) - gde3w_n(ji,jj,1) ) * r1_e1u(ji,jj)639 & * ( gde3w(ji+1,jj,1) - gde3w(ji,jj,1) ) * r1_e1u(ji,jj) 627 640 zvap = -zcoef0 * ( rhd (ji,jj+1,1) + rhd (ji,jj,1) + 2._wp * znad ) & 628 & * ( gde3w _n(ji,jj+1,1) - gde3w_n(ji,jj,1) ) * r1_e2v(ji,jj)641 & * ( gde3w(ji,jj+1,1) - gde3w(ji,jj,1) ) * r1_e2v(ji,jj) 629 642 ! add to the general momentum trend 630 ua(ji,jj,1) = ua(ji,jj,1) + (zhpi(ji,jj,1) + zuap) * umask(ji,jj,1)631 va(ji,jj,1) = va(ji,jj,1) + (zhpj(ji,jj,1) + zvap) * vmask(ji,jj,1)643 puu(ji,jj,1,Krhs) = puu(ji,jj,1,Krhs) + (zhpi(ji,jj,1) + zuap) * umask(ji,jj,1) 644 pvv(ji,jj,1,Krhs) = pvv(ji,jj,1,Krhs) + (zhpj(ji,jj,1) + zvap) * vmask(ji,jj,1) 632 645 END DO 633 646 END DO … … 641 654 ! hydrostatic pressure gradient along s-surfaces 642 655 zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) + zcoef0 / e1u(ji,jj) & 643 & * ( e3w _n(ji+1,jj,jk) * ( rhd(ji+1,jj,jk) + rhd(ji+1,jj,jk-1) + 2*znad ) * wmask(ji+1,jj,jk) &644 & - e3w _n(ji ,jj,jk) * ( rhd(ji ,jj,jk) + rhd(ji ,jj,jk-1) + 2*znad ) * wmask(ji ,jj,jk) )656 & * ( e3w(ji+1,jj,jk,Kmm) * ( rhd(ji+1,jj,jk) + rhd(ji+1,jj,jk-1) + 2*znad ) * wmask(ji+1,jj,jk) & 657 & - e3w(ji ,jj,jk,Kmm) * ( rhd(ji ,jj,jk) + rhd(ji ,jj,jk-1) + 2*znad ) * wmask(ji ,jj,jk) ) 645 658 zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1) + zcoef0 / e2v(ji,jj) & 646 & * ( e3w _n(ji,jj+1,jk) * ( rhd(ji,jj+1,jk) + rhd(ji,jj+1,jk-1) + 2*znad ) * wmask(ji,jj+1,jk) &647 & - e3w _n(ji,jj ,jk) * ( rhd(ji,jj, jk) + rhd(ji,jj ,jk-1) + 2*znad ) * wmask(ji,jj ,jk) )659 & * ( e3w(ji,jj+1,jk,Kmm) * ( rhd(ji,jj+1,jk) + rhd(ji,jj+1,jk-1) + 2*znad ) * wmask(ji,jj+1,jk) & 660 & - e3w(ji,jj ,jk,Kmm) * ( rhd(ji,jj, jk) + rhd(ji,jj ,jk-1) + 2*znad ) * wmask(ji,jj ,jk) ) 648 661 ! s-coordinate pressure gradient correction 649 662 zuap = -zcoef0 * ( rhd (ji+1,jj ,jk) + rhd (ji,jj,jk) + 2._wp * znad ) & 650 & * ( gde3w _n(ji+1,jj ,jk) - gde3w_n(ji,jj,jk) ) / e1u(ji,jj)663 & * ( gde3w(ji+1,jj ,jk) - gde3w(ji,jj,jk) ) / e1u(ji,jj) 651 664 zvap = -zcoef0 * ( rhd (ji ,jj+1,jk) + rhd (ji,jj,jk) + 2._wp * znad ) & 652 & * ( gde3w _n(ji ,jj+1,jk) - gde3w_n(ji,jj,jk) ) / e2v(ji,jj)665 & * ( gde3w(ji ,jj+1,jk) - gde3w(ji,jj,jk) ) / e2v(ji,jj) 653 666 ! add to the general momentum trend 654 ua(ji,jj,jk) = ua(ji,jj,jk) + (zhpi(ji,jj,jk) + zuap) * umask(ji,jj,jk)655 va(ji,jj,jk) = va(ji,jj,jk) + (zhpj(ji,jj,jk) + zvap) * vmask(ji,jj,jk)667 puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) + (zhpi(ji,jj,jk) + zuap) * umask(ji,jj,jk) 668 pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) + (zhpj(ji,jj,jk) + zvap) * vmask(ji,jj,jk) 656 669 END DO 657 670 END DO … … 661 674 662 675 663 SUBROUTINE hpg_djc( kt )676 SUBROUTINE hpg_djc( kt, Kmm, puu, pvv, Krhs ) 664 677 !!--------------------------------------------------------------------- 665 678 !! *** ROUTINE hpg_djc *** … … 669 682 !! Reference: Shchepetkin and McWilliams, J. Geophys. Res., 108(C3), 3090, 2003 670 683 !!---------------------------------------------------------------------- 671 INTEGER, INTENT(in) :: kt ! ocean time-step index 684 INTEGER , INTENT( in ) :: kt ! ocean time-step index 685 INTEGER , INTENT( in ) :: Kmm, Krhs ! ocean time level indices 686 REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) :: puu, pvv ! ocean velocities and RHS of momentum equation 672 687 !! 673 688 INTEGER :: ji, jj, jk ! dummy loop indices … … 687 702 DO jj = 2, jpjm1 688 703 DO ji = 2, jpim1 689 ll_tmp1 = MIN( ssh n(ji,jj) , sshn(ji+1,jj) ) > &704 ll_tmp1 = MIN( ssh(ji,jj,Kmm) , ssh(ji+1,jj,Kmm) ) > & 690 705 & MAX( -ht_0(ji,jj) , -ht_0(ji+1,jj) ) .AND. & 691 & MAX( ssh n(ji,jj) + ht_0(ji,jj), sshn(ji+1,jj) + ht_0(ji+1,jj) ) &706 & MAX( ssh(ji,jj,Kmm) + ht_0(ji,jj), ssh(ji+1,jj,Kmm) + ht_0(ji+1,jj) ) & 692 707 & > rn_wdmin1 + rn_wdmin2 693 ll_tmp2 = ( ABS( ssh n(ji,jj) - sshn(ji+1,jj) ) > 1.E-12 ) .AND. ( &694 & MAX( ssh n(ji,jj) , sshn(ji+1,jj) ) > &708 ll_tmp2 = ( ABS( ssh(ji,jj,Kmm) - ssh(ji+1,jj,Kmm) ) > 1.E-12 ) .AND. ( & 709 & MAX( ssh(ji,jj,Kmm) , ssh(ji+1,jj,Kmm) ) > & 695 710 & MAX( -ht_0(ji,jj) , -ht_0(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 ) 696 711 IF(ll_tmp1) THEN 697 712 zcpx(ji,jj) = 1.0_wp 698 713 ELSE IF(ll_tmp2) THEN 699 ! no worries about ssh n(ji+1,jj) - sshn(ji ,jj) = 0, it won't happen ! here700 zcpx(ji,jj) = ABS( (ssh n(ji+1,jj) + ht_0(ji+1,jj) - sshn(ji,jj) - ht_0(ji,jj)) &701 & / (ssh n(ji+1,jj) - sshn(ji ,jj)) )714 ! no worries about ssh(ji+1,jj,Kmm) - ssh(ji ,jj,Kmm) = 0, it won't happen ! here 715 zcpx(ji,jj) = ABS( (ssh(ji+1,jj,Kmm) + ht_0(ji+1,jj) - ssh(ji,jj,Kmm) - ht_0(ji,jj)) & 716 & / (ssh(ji+1,jj,Kmm) - ssh(ji ,jj,Kmm)) ) 702 717 ELSE 703 718 zcpx(ji,jj) = 0._wp 704 719 END IF 705 720 706 ll_tmp1 = MIN( ssh n(ji,jj) , sshn(ji,jj+1) ) > &721 ll_tmp1 = MIN( ssh(ji,jj,Kmm) , ssh(ji,jj+1,Kmm) ) > & 707 722 & MAX( -ht_0(ji,jj) , -ht_0(ji,jj+1) ) .AND. & 708 & MAX( ssh n(ji,jj) + ht_0(ji,jj), sshn(ji,jj+1) + ht_0(ji,jj+1) ) &723 & MAX( ssh(ji,jj,Kmm) + ht_0(ji,jj), ssh(ji,jj+1,Kmm) + ht_0(ji,jj+1) ) & 709 724 & > rn_wdmin1 + rn_wdmin2 710 ll_tmp2 = ( ABS( ssh n(ji,jj) - sshn(ji,jj+1) ) > 1.E-12 ) .AND. ( &711 & MAX( ssh n(ji,jj) , sshn(ji,jj+1) ) > &725 ll_tmp2 = ( ABS( ssh(ji,jj,Kmm) - ssh(ji,jj+1,Kmm) ) > 1.E-12 ) .AND. ( & 726 & MAX( ssh(ji,jj,Kmm) , ssh(ji,jj+1,Kmm) ) > & 712 727 & MAX( -ht_0(ji,jj) , -ht_0(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 ) 713 728 … … 715 730 zcpy(ji,jj) = 1.0_wp 716 731 ELSE IF(ll_tmp2) THEN 717 ! no worries about ssh n(ji,jj+1) - sshn(ji,jj) = 0, it won't happen ! here718 zcpy(ji,jj) = ABS( (ssh n(ji,jj+1) + ht_0(ji,jj+1) - sshn(ji,jj) - ht_0(ji,jj)) &719 & / (ssh n(ji,jj+1) - sshn(ji,jj)) )732 ! no worries about ssh(ji,jj+1,Kmm) - ssh(ji,jj ,Kmm) = 0, it won't happen ! here 733 zcpy(ji,jj) = ABS( (ssh(ji,jj+1,Kmm) + ht_0(ji,jj+1) - ssh(ji,jj,Kmm) - ht_0(ji,jj)) & 734 & / (ssh(ji,jj+1,Kmm) - ssh(ji,jj ,Kmm)) ) 720 735 ELSE 721 736 zcpy(ji,jj) = 0._wp … … 747 762 DO ji = fs_2, fs_jpim1 ! vector opt. 748 763 drhoz(ji,jj,jk) = rhd (ji ,jj ,jk) - rhd (ji,jj,jk-1) 749 dzz (ji,jj,jk) = gde3w _n(ji ,jj ,jk) - gde3w_n(ji,jj,jk-1)764 dzz (ji,jj,jk) = gde3w(ji ,jj ,jk) - gde3w(ji,jj,jk-1) 750 765 drhox(ji,jj,jk) = rhd (ji+1,jj ,jk) - rhd (ji,jj,jk ) 751 dzx (ji,jj,jk) = gde3w _n(ji+1,jj ,jk) - gde3w_n(ji,jj,jk )766 dzx (ji,jj,jk) = gde3w(ji+1,jj ,jk) - gde3w(ji,jj,jk ) 752 767 drhoy(ji,jj,jk) = rhd (ji ,jj+1,jk) - rhd (ji,jj,jk ) 753 dzy (ji,jj,jk) = gde3w _n(ji ,jj+1,jk) - gde3w_n(ji,jj,jk )768 dzy (ji,jj,jk) = gde3w(ji ,jj+1,jk) - gde3w(ji,jj,jk ) 754 769 END DO 755 770 END DO … … 838 853 DO jj = 2, jpjm1 839 854 DO ji = fs_2, fs_jpim1 ! vector opt. 840 rho_k(ji,jj,1) = -grav * ( e3w _n(ji,jj,1) - gde3w_n(ji,jj,1) ) &855 rho_k(ji,jj,1) = -grav * ( e3w(ji,jj,1,Kmm) - gde3w(ji,jj,1) ) & 841 856 & * ( rhd(ji,jj,1) & 842 857 & + 0.5_wp * ( rhd (ji,jj,2) - rhd (ji,jj,1) ) & 843 & * ( e3w _n (ji,jj,1) - gde3w_n(ji,jj,1) ) &844 & / ( gde3w _n(ji,jj,2) - gde3w_n(ji,jj,1) ) )858 & * ( e3w (ji,jj,1,Kmm) - gde3w(ji,jj,1) ) & 859 & / ( gde3w(ji,jj,2) - gde3w(ji,jj,1) ) ) 845 860 END DO 846 861 END DO … … 854 869 855 870 rho_k(ji,jj,jk) = zcoef0 * ( rhd (ji,jj,jk) + rhd (ji,jj,jk-1) ) & 856 & * ( gde3w _n(ji,jj,jk) - gde3w_n(ji,jj,jk-1) ) &871 & * ( gde3w(ji,jj,jk) - gde3w(ji,jj,jk-1) ) & 857 872 & - grav * z1_10 * ( & 858 873 & ( drhow (ji,jj,jk) - drhow (ji,jj,jk-1) ) & 859 & * ( gde3w _n(ji,jj,jk) - gde3w_n(ji,jj,jk-1) - z1_12 * ( dzw (ji,jj,jk) + dzw (ji,jj,jk-1) ) ) &874 & * ( gde3w(ji,jj,jk) - gde3w(ji,jj,jk-1) - z1_12 * ( dzw (ji,jj,jk) + dzw (ji,jj,jk-1) ) ) & 860 875 & - ( dzw (ji,jj,jk) - dzw (ji,jj,jk-1) ) & 861 876 & * ( rhd (ji,jj,jk) - rhd (ji,jj,jk-1) - z1_12 * ( drhow(ji,jj,jk) + drhow(ji,jj,jk-1) ) ) & … … 863 878 864 879 rho_i(ji,jj,jk) = zcoef0 * ( rhd (ji+1,jj,jk) + rhd (ji,jj,jk) ) & 865 & * ( gde3w _n(ji+1,jj,jk) - gde3w_n(ji,jj,jk) ) &880 & * ( gde3w(ji+1,jj,jk) - gde3w(ji,jj,jk) ) & 866 881 & - grav* z1_10 * ( & 867 882 & ( drhou (ji+1,jj,jk) - drhou (ji,jj,jk) ) & 868 & * ( gde3w _n(ji+1,jj,jk) - gde3w_n(ji,jj,jk) - z1_12 * ( dzu (ji+1,jj,jk) + dzu (ji,jj,jk) ) ) &883 & * ( gde3w(ji+1,jj,jk) - gde3w(ji,jj,jk) - z1_12 * ( dzu (ji+1,jj,jk) + dzu (ji,jj,jk) ) ) & 869 884 & - ( dzu (ji+1,jj,jk) - dzu (ji,jj,jk) ) & 870 885 & * ( rhd (ji+1,jj,jk) - rhd (ji,jj,jk) - z1_12 * ( drhou(ji+1,jj,jk) + drhou(ji,jj,jk) ) ) & … … 872 887 873 888 rho_j(ji,jj,jk) = zcoef0 * ( rhd (ji,jj+1,jk) + rhd (ji,jj,jk) ) & 874 & * ( gde3w _n(ji,jj+1,jk) - gde3w_n(ji,jj,jk) ) &889 & * ( gde3w(ji,jj+1,jk) - gde3w(ji,jj,jk) ) & 875 890 & - grav* z1_10 * ( & 876 891 & ( drhov (ji,jj+1,jk) - drhov (ji,jj,jk) ) & 877 & * ( gde3w _n(ji,jj+1,jk) - gde3w_n(ji,jj,jk) - z1_12 * ( dzv (ji,jj+1,jk) + dzv (ji,jj,jk) ) ) &892 & * ( gde3w(ji,jj+1,jk) - gde3w(ji,jj,jk) - z1_12 * ( dzv (ji,jj+1,jk) + dzv (ji,jj,jk) ) ) & 878 893 & - ( dzv (ji,jj+1,jk) - dzv (ji,jj,jk) ) & 879 894 & * ( rhd (ji,jj+1,jk) - rhd (ji,jj,jk) - z1_12 * ( drhov(ji,jj+1,jk) + drhov(ji,jj,jk) ) ) & … … 897 912 ENDIF 898 913 ! add to the general momentum trend 899 ua(ji,jj,1) = ua(ji,jj,1) + zhpi(ji,jj,1)900 va(ji,jj,1) = va(ji,jj,1) + zhpj(ji,jj,1)914 puu(ji,jj,1,Krhs) = puu(ji,jj,1,Krhs) + zhpi(ji,jj,1) 915 pvv(ji,jj,1,Krhs) = pvv(ji,jj,1,Krhs) + zhpj(ji,jj,1) 901 916 END DO 902 917 END DO … … 920 935 ENDIF 921 936 ! add to the general momentum trend 922 ua(ji,jj,jk) = ua(ji,jj,jk) + zhpi(ji,jj,jk)923 va(ji,jj,jk) = va(ji,jj,jk) + zhpj(ji,jj,jk)937 puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) + zhpi(ji,jj,jk) 938 pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) + zhpj(ji,jj,jk) 924 939 END DO 925 940 END DO … … 931 946 932 947 933 SUBROUTINE hpg_prj( kt )948 SUBROUTINE hpg_prj( kt, Kmm, puu, pvv, Krhs ) 934 949 !!--------------------------------------------------------------------- 935 950 !! *** ROUTINE hpg_prj *** … … 940 955 !! all vertical coordinate systems 941 956 !! 942 !! ** Action : - Update ( ua,va) with the now hydrastatic pressure trend957 !! ** Action : - Update (puu(:,:,:,Krhs),pvv(:,:,:,Krhs)) with the now hydrastatic pressure trend 943 958 !!---------------------------------------------------------------------- 944 959 INTEGER, PARAMETER :: polynomial_type = 1 ! 1: cubic spline, 2: linear 945 INTEGER, INTENT(in) :: kt ! ocean time-step index 960 INTEGER , INTENT( in ) :: kt ! ocean time-step index 961 INTEGER , INTENT( in ) :: Kmm, Krhs ! ocean time level indices 962 REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) :: puu, pvv ! ocean velocities and RHS of momentum equation 946 963 !! 947 964 INTEGER :: ji, jj, jk, jkk ! dummy loop indices … … 975 992 DO jj = 2, jpjm1 976 993 DO ji = 2, jpim1 977 ll_tmp1 = MIN( ssh n(ji,jj) , sshn(ji+1,jj) ) > &994 ll_tmp1 = MIN( ssh(ji,jj,Kmm) , ssh(ji+1,jj,Kmm) ) > & 978 995 & MAX( -ht_0(ji,jj) , -ht_0(ji+1,jj) ) .AND. & 979 & MAX( ssh n(ji,jj) + ht_0(ji,jj), sshn(ji+1,jj) + ht_0(ji+1,jj) ) &996 & MAX( ssh(ji,jj,Kmm) + ht_0(ji,jj), ssh(ji+1,jj,Kmm) + ht_0(ji+1,jj) ) & 980 997 & > rn_wdmin1 + rn_wdmin2 981 ll_tmp2 = ( ABS( ssh n(ji,jj) - sshn(ji+1,jj) ) > 1.E-12 ) .AND. ( &982 & MAX( ssh n(ji,jj) , sshn(ji+1,jj) ) > &998 ll_tmp2 = ( ABS( ssh(ji,jj,Kmm) - ssh(ji+1,jj,Kmm) ) > 1.E-12 ) .AND. ( & 999 & MAX( ssh(ji,jj,Kmm) , ssh(ji+1,jj,Kmm) ) > & 983 1000 & MAX( -ht_0(ji,jj) , -ht_0(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 ) 984 1001 … … 986 1003 zcpx(ji,jj) = 1.0_wp 987 1004 ELSE IF(ll_tmp2) THEN 988 ! no worries about ssh n(ji+1,jj) - sshn(ji ,jj) = 0, it won't happen ! here989 zcpx(ji,jj) = ABS( (ssh n(ji+1,jj) + ht_0(ji+1,jj) - sshn(ji,jj) - ht_0(ji,jj)) &990 & / (ssh n(ji+1,jj) - sshn(ji ,jj)) )1005 ! no worries about ssh(ji+1,jj,Kmm) - ssh(ji ,jj,Kmm) = 0, it won't happen ! here 1006 zcpx(ji,jj) = ABS( (ssh(ji+1,jj,Kmm) + ht_0(ji+1,jj) - ssh(ji,jj,Kmm) - ht_0(ji,jj)) & 1007 & / (ssh(ji+1,jj,Kmm) - ssh(ji ,jj,Kmm)) ) 991 1008 992 1009 zcpx(ji,jj) = max(min( zcpx(ji,jj) , 1.0_wp),0.0_wp) … … 995 1012 END IF 996 1013 997 ll_tmp1 = MIN( ssh n(ji,jj) , sshn(ji,jj+1) ) > &1014 ll_tmp1 = MIN( ssh(ji,jj,Kmm) , ssh(ji,jj+1,Kmm) ) > & 998 1015 & MAX( -ht_0(ji,jj) , -ht_0(ji,jj+1) ) .AND. & 999 & MAX( ssh n(ji,jj) + ht_0(ji,jj), sshn(ji,jj+1) + ht_0(ji,jj+1) ) &1016 & MAX( ssh(ji,jj,Kmm) + ht_0(ji,jj), ssh(ji,jj+1,Kmm) + ht_0(ji,jj+1) ) & 1000 1017 & > rn_wdmin1 + rn_wdmin2 1001 ll_tmp2 = ( ABS( ssh n(ji,jj) - sshn(ji,jj+1) ) > 1.E-12 ) .AND. ( &1002 & MAX( ssh n(ji,jj) , sshn(ji,jj+1) ) > &1018 ll_tmp2 = ( ABS( ssh(ji,jj,Kmm) - ssh(ji,jj+1,Kmm) ) > 1.E-12 ) .AND. ( & 1019 & MAX( ssh(ji,jj,Kmm) , ssh(ji,jj+1,Kmm) ) > & 1003 1020 & MAX( -ht_0(ji,jj) , -ht_0(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 ) 1004 1021 … … 1006 1023 zcpy(ji,jj) = 1.0_wp 1007 1024 ELSE IF(ll_tmp2) THEN 1008 ! no worries about ssh n(ji,jj+1) - sshn(ji,jj) = 0, it won't happen ! here1009 zcpy(ji,jj) = ABS( (ssh n(ji,jj+1) + ht_0(ji,jj+1) - sshn(ji,jj) - ht_0(ji,jj)) &1010 & / (ssh n(ji,jj+1) - sshn(ji,jj)) )1025 ! no worries about ssh(ji,jj+1,Kmm) - ssh(ji,jj ,Kmm) = 0, it won't happen ! here 1026 zcpy(ji,jj) = ABS( (ssh(ji,jj+1,Kmm) + ht_0(ji,jj+1) - ssh(ji,jj,Kmm) - ht_0(ji,jj)) & 1027 & / (ssh(ji,jj+1,Kmm) - ssh(ji,jj ,Kmm)) ) 1011 1028 zcpy(ji,jj) = max(min( zcpy(ji,jj) , 1.0_wp),0.0_wp) 1012 1029 … … 1031 1048 ELSEIF( jk < jpkm1 ) THEN 1032 1049 DO jkk = jk+1, jpk 1033 zrhh(ji,jj,jkk) = interp1(gde3w _n(ji,jj,jkk ), gde3w_n(ji,jj,jkk-1), &1034 & gde3w _n(ji,jj,jkk-2), rhd (ji,jj,jkk-1), rhd(ji,jj,jkk-2))1050 zrhh(ji,jj,jkk) = interp1(gde3w(ji,jj,jkk ), gde3w(ji,jj,jkk-1), & 1051 & gde3w(ji,jj,jkk-2), rhd (ji,jj,jkk-1), rhd(ji,jj,jkk-2)) 1035 1052 END DO 1036 1053 ENDIF … … 1041 1058 DO jj = 1, jpj 1042 1059 DO ji = 1, jpi 1043 zdept(ji,jj,1) = 0.5_wp * e3w _n(ji,jj,1) - sshn(ji,jj) * znad1060 zdept(ji,jj,1) = 0.5_wp * e3w(ji,jj,1,Kmm) - ssh(ji,jj,Kmm) * znad 1044 1061 END DO 1045 1062 END DO … … 1048 1065 DO jj = 1, jpj 1049 1066 DO ji = 1, jpi 1050 zdept(ji,jj,jk) = zdept(ji,jj,jk-1) + e3w _n(ji,jj,jk)1067 zdept(ji,jj,jk) = zdept(ji,jj,jk-1) + e3w(ji,jj,jk,Kmm) 1051 1068 END DO 1052 1069 END DO … … 1065 1082 DO ji = 2, jpi 1066 1083 zrhdt1 = zrhh(ji,jj,1) - interp3( zdept(ji,jj,1), asp(ji,jj,1), bsp(ji,jj,1), & 1067 & csp(ji,jj,1), dsp(ji,jj,1) ) * 0.25_wp * e3w _n(ji,jj,1)1084 & csp(ji,jj,1), dsp(ji,jj,1) ) * 0.25_wp * e3w(ji,jj,1,Kmm) 1068 1085 1069 1086 ! assuming linear profile across the top half surface layer 1070 zhpi(ji,jj,1) = 0.5_wp * e3w _n(ji,jj,1) * zrhdt11087 zhpi(ji,jj,1) = 0.5_wp * e3w(ji,jj,1,Kmm) * zrhdt1 1071 1088 END DO 1072 1089 END DO … … 1090 1107 DO ji = 2, jpim1 1091 1108 !!gm BUG ? if it is ssh at u- & v-point then it should be: 1092 ! zsshu_n(ji,jj) = (e1e2t(ji,jj) * ssh n(ji,jj) + e1e2t(ji+1,jj) * sshn(ji+1,jj)) * &1109 ! zsshu_n(ji,jj) = (e1e2t(ji,jj) * ssh(ji,jj,Kmm) + e1e2t(ji+1,jj) * ssh(ji+1,jj,Kmm)) * & 1093 1110 ! & r1_e1e2u(ji,jj) * umask(ji,jj,1) * 0.5_wp 1094 ! zsshv_n(ji,jj) = (e1e2t(ji,jj) * ssh n(ji,jj) + e1e2t(ji,jj+1) * sshn(ji,jj+1)) * &1111 ! zsshv_n(ji,jj) = (e1e2t(ji,jj) * ssh(ji,jj,Kmm) + e1e2t(ji,jj+1) * ssh(ji,jj+1,Kmm)) * & 1095 1112 ! & r1_e1e2v(ji,jj) * vmask(ji,jj,1) * 0.5_wp 1096 1113 !!gm not this: 1097 zsshu_n(ji,jj) = (e1e2u(ji,jj) * ssh n(ji,jj) + e1e2u(ji+1, jj) * sshn(ji+1,jj)) * &1114 zsshu_n(ji,jj) = (e1e2u(ji,jj) * ssh(ji,jj,Kmm) + e1e2u(ji+1, jj) * ssh(ji+1,jj,Kmm)) * & 1098 1115 & r1_e1e2u(ji,jj) * umask(ji,jj,1) * 0.5_wp 1099 zsshv_n(ji,jj) = (e1e2v(ji,jj) * ssh n(ji,jj) + e1e2v(ji+1, jj) * sshn(ji,jj+1)) * &1116 zsshv_n(ji,jj) = (e1e2v(ji,jj) * ssh(ji,jj,Kmm) + e1e2v(ji+1, jj) * ssh(ji,jj+1,Kmm)) * & 1100 1117 & r1_e1e2v(ji,jj) * vmask(ji,jj,1) * 0.5_wp 1101 1118 END DO … … 1106 1123 DO jj = 2, jpjm1 1107 1124 DO ji = 2, jpim1 1108 zu(ji,jj,1) = - ( e3u _n(ji,jj,1) - zsshu_n(ji,jj) * znad)1109 zv(ji,jj,1) = - ( e3v _n(ji,jj,1) - zsshv_n(ji,jj) * znad)1125 zu(ji,jj,1) = - ( e3u(ji,jj,1,Kmm) - zsshu_n(ji,jj) * znad) 1126 zv(ji,jj,1) = - ( e3v(ji,jj,1,Kmm) - zsshv_n(ji,jj) * znad) 1110 1127 END DO 1111 1128 END DO … … 1114 1131 DO jj = 2, jpjm1 1115 1132 DO ji = 2, jpim1 1116 zu(ji,jj,jk) = zu(ji,jj,jk-1) - e3u _n(ji,jj,jk)1117 zv(ji,jj,jk) = zv(ji,jj,jk-1) - e3v _n(ji,jj,jk)1133 zu(ji,jj,jk) = zu(ji,jj,jk-1) - e3u(ji,jj,jk,Kmm) 1134 zv(ji,jj,jk) = zv(ji,jj,jk-1) - e3v(ji,jj,jk,Kmm) 1118 1135 END DO 1119 1136 END DO … … 1123 1140 DO jj = 2, jpjm1 1124 1141 DO ji = 2, jpim1 1125 zu(ji,jj,jk) = zu(ji,jj,jk) + 0.5_wp * e3u _n(ji,jj,jk)1126 zv(ji,jj,jk) = zv(ji,jj,jk) + 0.5_wp * e3v _n(ji,jj,jk)1142 zu(ji,jj,jk) = zu(ji,jj,jk) + 0.5_wp * e3u(ji,jj,jk,Kmm) 1143 zv(ji,jj,jk) = zv(ji,jj,jk) + 0.5_wp * e3v(ji,jj,jk,Kmm) 1127 1144 END DO 1128 1145 END DO … … 1176 1193 DO WHILE ( -zdept(jid,jj,jk1) < zuijk ) 1177 1194 IF( jk1 == 1 ) THEN 1178 zdeps = zdept(jid,jj,1) + MIN(zuijk, ssh n(jid,jj)*znad)1195 zdeps = zdept(jid,jj,1) + MIN(zuijk, ssh(jid,jj,Kmm)*znad) 1179 1196 zrhdt1 = zrhh(jid,jj,1) - interp3(zdept(jid,jj,1), asp(jid,jj,1), & 1180 1197 bsp(jid,jj,1), csp(jid,jj,1), & … … 1196 1213 IF( .NOT.ln_linssh ) THEN 1197 1214 zdpdx2 = zcoef0 * r1_e1u(ji,jj) * & 1198 & ( REAL(jis-jid, wp) * (zpwes + zpwed) + (ssh n(ji+1,jj)-sshn(ji,jj)) )1215 & ( REAL(jis-jid, wp) * (zpwes + zpwed) + (ssh(ji+1,jj,Kmm)-ssh(ji,jj,Kmm)) ) 1199 1216 ELSE 1200 1217 zdpdx2 = zcoef0 * r1_e1u(ji,jj) * REAL(jis-jid, wp) * (zpwes + zpwed) … … 1204 1221 zdpdx2 = zdpdx2 * zcpx(ji,jj) * wdrampu(ji,jj) 1205 1222 ENDIF 1206 ua(ji,jj,jk) = ua(ji,jj,jk) + (zdpdx1 + zdpdx2) * umask(ji,jj,jk)1223 puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) + (zdpdx1 + zdpdx2) * umask(ji,jj,jk) 1207 1224 ENDIF 1208 1225 … … 1234 1251 DO WHILE ( -zdept(ji,jjd,jk1) < zvijk ) 1235 1252 IF( jk1 == 1 ) THEN 1236 zdeps = zdept(ji,jjd,1) + MIN(zvijk, ssh n(ji,jjd)*znad)1253 zdeps = zdept(ji,jjd,1) + MIN(zvijk, ssh(ji,jjd,Kmm)*znad) 1237 1254 zrhdt1 = zrhh(ji,jjd,1) - interp3(zdept(ji,jjd,1), asp(ji,jjd,1), & 1238 1255 bsp(ji,jjd,1), csp(ji,jjd,1), & … … 1255 1272 IF( .NOT.ln_linssh ) THEN 1256 1273 zdpdy2 = zcoef0 * r1_e2v(ji,jj) * & 1257 ( REAL(jjs-jjd, wp) * (zpnss + zpnsd) + (ssh n(ji,jj+1)-sshn(ji,jj)) )1274 ( REAL(jjs-jjd, wp) * (zpnss + zpnsd) + (ssh(ji,jj+1,Kmm)-ssh(ji,jj,Kmm)) ) 1258 1275 ELSE 1259 1276 zdpdy2 = zcoef0 * r1_e2v(ji,jj) * REAL(jjs-jjd, wp) * (zpnss + zpnsd ) … … 1264 1281 ENDIF 1265 1282 1266 va(ji,jj,jk) = va(ji,jj,jk) + (zdpdy1 + zdpdy2) * vmask(ji,jj,jk)1283 pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) + (zdpdy1 + zdpdy2) * vmask(ji,jj,jk) 1267 1284 ENDIF 1268 1285 !
Note: See TracChangeset
for help on using the changeset viewer.