- Timestamp:
- 2016-07-19T10:38:35+02:00 (8 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/NERC/dev_r5549_BDY_ZEROGRAD/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90
r5224 r6808 33 33 USE sbc_oce ! surface variable (only for the flag with ice shelf) 34 34 USE dom_oce ! ocean space and time domain 35 USE wet_dry ! wetting and drying 35 36 USE phycst ! physical constants 36 37 USE trd_oce ! trends: ocean variables 37 38 USE trddyn ! trend manager: dynamics 39 !jc USE zpshde ! partial step: hor. derivative (zps_hde routine) 38 40 ! 39 41 USE in_out_manager ! I/O manager … … 44 46 USE wrk_nemo ! Memory Allocation 45 47 USE timing ! Timing 48 USE iom 46 49 47 50 IMPLICIT NONE … … 51 54 PUBLIC dyn_hpg_init ! routine called by opa module 52 55 53 ! !!* Namelist namdyn_hpg : hydrostatic pressure gradient 54 LOGICAL , PUBLIC :: ln_hpg_zco !: z-coordinate - full steps 55 LOGICAL , PUBLIC :: ln_hpg_zps !: z-coordinate - partial steps (interpolation) 56 LOGICAL , PUBLIC :: ln_hpg_sco !: s-coordinate (standard jacobian formulation) 57 LOGICAL , PUBLIC :: ln_hpg_djc !: s-coordinate (Density Jacobian with Cubic polynomial) 58 LOGICAL , PUBLIC :: ln_hpg_prj !: s-coordinate (Pressure Jacobian scheme) 59 LOGICAL , PUBLIC :: ln_hpg_isf !: s-coordinate similar to sco modify for isf 60 LOGICAL , PUBLIC :: ln_dynhpg_imp !: semi-implicite hpg flag 56 ! !!* Namelist namdyn_hpg : hydrostatic pressure gradient 57 LOGICAL , PUBLIC :: ln_hpg_zco !: z-coordinate - full steps 58 LOGICAL , PUBLIC :: ln_hpg_zps !: z-coordinate - partial steps (interpolation) 59 LOGICAL , PUBLIC :: ln_hpg_sco !: s-coordinate (standard jacobian formulation) 60 LOGICAL , PUBLIC :: ln_hpg_djc !: s-coordinate (Density Jacobian with Cubic polynomial) 61 LOGICAL , PUBLIC :: ln_hpg_prj !: s-coordinate (Pressure Jacobian scheme) 62 LOGICAL , PUBLIC :: ln_hpg_isf !: s-coordinate similar to sco modify for isf 61 63 62 64 INTEGER , PUBLIC :: nhpg = 0 ! = 0 to 7, type of pressure gradient scheme used ! (deduced from ln_hpg_... flags) (PUBLIC for TAM) 63 65 64 66 !! * Substitutions 65 # include "domzgr_substitute.h90"66 67 # include "vectopt_loop_substitute.h90" 67 68 !!---------------------------------------------------------------------- … … 131 132 INTEGER :: ios ! Local integer output status for namelist read 132 133 !! 134 INTEGER :: ji, jj, jk, ikt ! dummy loop indices ISF 135 REAL(wp) :: znad 136 REAL(wp), POINTER, DIMENSION(:,:,:) :: ztstop, zrhd ! hypothesys on isf density 137 REAL(wp), POINTER, DIMENSION(:,:) :: zrhdtop_isf ! density at bottom of ISF 138 REAL(wp), POINTER, DIMENSION(:,:) :: ziceload ! density at bottom of ISF 139 !! 133 140 NAMELIST/namdyn_hpg/ ln_hpg_zco, ln_hpg_zps, ln_hpg_sco, & 134 & ln_hpg_djc, ln_hpg_prj, ln_hpg_isf , ln_dynhpg_imp141 & ln_hpg_djc, ln_hpg_prj, ln_hpg_isf 135 142 !!---------------------------------------------------------------------- 136 143 ! 137 144 REWIND( numnam_ref ) ! Namelist namdyn_hpg in reference namelist : Hydrostatic pressure gradient 138 145 READ ( numnam_ref, namdyn_hpg, IOSTAT = ios, ERR = 901) 139 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdyn_hpg in reference namelist', lwp )140 146 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdyn_hpg in reference namelist', lwp ) 147 ! 141 148 REWIND( numnam_cfg ) ! Namelist namdyn_hpg in configuration namelist : Hydrostatic pressure gradient 142 149 READ ( numnam_cfg, namdyn_hpg, IOSTAT = ios, ERR = 902 ) 143 902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdyn_hpg in configuration namelist', lwp )150 902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdyn_hpg in configuration namelist', lwp ) 144 151 IF(lwm) WRITE ( numond, namdyn_hpg ) 145 152 ! … … 155 162 WRITE(numout,*) ' s-coord. (Density Jacobian: Cubic polynomial) ln_hpg_djc = ', ln_hpg_djc 156 163 WRITE(numout,*) ' s-coord. (Pressure Jacobian: Cubic polynomial) ln_hpg_prj = ', ln_hpg_prj 157 WRITE(numout,*) ' time stepping: centered (F) or semi-implicit (T) ln_dynhpg_imp = ', ln_dynhpg_imp158 164 ENDIF 159 165 ! … … 163 169 & either ln_hpg_sco or ln_hpg_prj instead') 164 170 ! 165 IF( lk_vvl .AND. .NOT. (ln_hpg_sco.OR.ln_hpg_prj.OR.ln_hpg_isf) )&166 & CALL ctl_stop('dyn_hpg_init : variable volume key_vvl requires:&167 & the standard jacobian formulation hpg_sco or&168 & the pressure jacobian formulation hpg_prj')171 IF( .NOT.ln_linssh .AND. .NOT.(ln_hpg_sco.OR.ln_hpg_prj.OR.ln_hpg_isf) ) & 172 & CALL ctl_stop('dyn_hpg_init : non-linear free surface requires either ', & 173 & ' the standard jacobian formulation hpg_sco or ' , & 174 & ' the pressure jacobian formulation hpg_prj' ) 169 175 170 176 IF( ln_hpg_isf .AND. .NOT. ln_isfcav ) & … … 191 197 IF( ioptio /= 1 ) CALL ctl_stop( 'NO or several hydrostatic pressure gradient options used' ) 192 198 ! 193 ! initialisation of ice load 194 riceload(:,:)=0.0 199 ! initialisation of ice shelf load 200 IF ( .NOT. ln_isfcav ) riceload(:,:)=0.0 201 IF ( ln_isfcav ) THEN 202 CALL wrk_alloc( jpi,jpj, 2, ztstop) 203 CALL wrk_alloc( jpi,jpj,jpk, zrhd ) 204 CALL wrk_alloc( jpi,jpj, zrhdtop_isf, ziceload) 205 ! 206 IF(lwp) WRITE(numout,*) 207 IF(lwp) WRITE(numout,*) 'dyn:hpg_isf : hydrostatic pressure gradient trend for ice shelf' 208 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~' 209 210 ! To use density and not density anomaly 211 znad=1._wp 212 213 ! assume water displaced by the ice shelf is at T=-1.9 and S=34.4 (rude) 214 ztstop(:,:,1)=-1.9_wp ; ztstop(:,:,2)=34.4_wp 215 216 ! compute density of the water displaced by the ice shelf 217 DO jk = 1, jpk 218 CALL eos(ztstop(:,:,:),gdept_n(:,:,jk),zrhd(:,:,jk)) 219 END DO 220 221 ! compute rhd at the ice/oce interface (ice shelf side) 222 CALL eos(ztstop,risfdep,zrhdtop_isf) 223 224 ! Surface value + ice shelf gradient 225 ! compute pressure due to ice shelf load (used to compute hpgi/j for all the level from 1 to miku/v) 226 ! divided by 2 later 227 ziceload = 0._wp 228 DO jj = 1, jpj 229 DO ji = 1, jpi 230 ikt=mikt(ji,jj) 231 ziceload(ji,jj) = ziceload(ji,jj) + (znad + zrhd(ji,jj,1) ) * e3w_n(ji,jj,1) * (1._wp - tmask(ji,jj,1)) 232 DO jk=2,ikt-1 233 ziceload(ji,jj) = ziceload(ji,jj) + (2._wp * znad + zrhd(ji,jj,jk-1) + zrhd(ji,jj,jk)) * e3w_n(ji,jj,jk) & 234 & * (1._wp - tmask(ji,jj,jk)) 235 END DO 236 IF (ikt >= 2) ziceload(ji,jj) = ziceload(ji,jj) + (2._wp * znad + zrhdtop_isf(ji,jj) + zrhd(ji,jj,ikt-1)) & 237 & * ( risfdep(ji,jj) - gdept_1d(ikt-1) ) 238 END DO 239 END DO 240 riceload(:,:)=ziceload(:,:) ! need to be saved for diaar5 241 242 CALL wrk_dealloc( jpi,jpj, 2, ztstop) 243 CALL wrk_dealloc( jpi,jpj,jpk, zrhd ) 244 CALL wrk_dealloc( jpi,jpj, zrhdtop_isf, ziceload) 245 END IF 195 246 ! 196 247 END SUBROUTINE dyn_hpg_init … … 214 265 !!---------------------------------------------------------------------- 215 266 INTEGER, INTENT(in) :: kt ! ocean time-step index 216 ! !267 ! 217 268 INTEGER :: ji, jj, jk ! dummy loop indices 218 269 REAL(wp) :: zcoef0, zcoef1 ! temporary scalars … … 220 271 !!---------------------------------------------------------------------- 221 272 ! 222 CALL wrk_alloc( jpi,jpj,jpk, zhpi, zhpj )273 CALL wrk_alloc( jpi,jpj,jpk, zhpi, zhpj ) 223 274 ! 224 275 IF( kt == nit000 ) THEN … … 233 284 DO jj = 2, jpjm1 234 285 DO ji = fs_2, fs_jpim1 ! vector opt. 235 zcoef1 = zcoef0 * fse3w(ji,jj,1)286 zcoef1 = zcoef0 * e3w_n(ji,jj,1) 236 287 ! hydrostatic pressure gradient 237 zhpi(ji,jj,1) = zcoef1 * ( rhd(ji+1,jj,1) - rhd(ji,jj,1) ) /e1u(ji,jj)238 zhpj(ji,jj,1) = zcoef1 * ( rhd(ji,jj+1,1) - rhd(ji,jj,1) ) /e2v(ji,jj)288 zhpi(ji,jj,1) = zcoef1 * ( rhd(ji+1,jj,1) - rhd(ji,jj,1) ) * r1_e1u(ji,jj) 289 zhpj(ji,jj,1) = zcoef1 * ( rhd(ji,jj+1,1) - rhd(ji,jj,1) ) * r1_e2v(ji,jj) 239 290 ! add to the general momentum trend 240 291 ua(ji,jj,1) = ua(ji,jj,1) + zhpi(ji,jj,1) … … 248 299 DO jj = 2, jpjm1 249 300 DO ji = fs_2, fs_jpim1 ! vector opt. 250 zcoef1 = zcoef0 * fse3w(ji,jj,jk)301 zcoef1 = zcoef0 * e3w_n(ji,jj,jk) 251 302 ! hydrostatic pressure gradient 252 303 zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) & 253 & + zcoef1 * ( ( rhd(ji+1,jj,jk)+rhd(ji+1,jj,jk-1) ) &254 & - ( rhd(ji ,jj,jk)+rhd(ji ,jj,jk-1) ) ) /e1u(ji,jj)304 & + zcoef1 * ( ( rhd(ji+1,jj,jk)+rhd(ji+1,jj,jk-1) ) & 305 & - ( rhd(ji ,jj,jk)+rhd(ji ,jj,jk-1) ) ) * r1_e1u(ji,jj) 255 306 256 307 zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1) & 257 & + zcoef1 * ( ( rhd(ji,jj+1,jk)+rhd(ji,jj+1,jk-1) ) &258 & - ( rhd(ji,jj, jk)+rhd(ji,jj ,jk-1) ) ) /e2v(ji,jj)308 & + zcoef1 * ( ( rhd(ji,jj+1,jk)+rhd(ji,jj+1,jk-1) ) & 309 & - ( rhd(ji,jj, jk)+rhd(ji,jj ,jk-1) ) ) * r1_e2v(ji,jj) 259 310 ! add to the general momentum trend 260 311 ua(ji,jj,jk) = ua(ji,jj,jk) + zhpi(ji,jj,jk) … … 264 315 END DO 265 316 ! 266 CALL wrk_dealloc( jpi,jpj,jpk, zhpi, zhpj )317 CALL wrk_dealloc( jpi,jpj,jpk, zhpi, zhpj ) 267 318 ! 268 319 END SUBROUTINE hpg_zco … … 285 336 !!---------------------------------------------------------------------- 286 337 ! 287 CALL wrk_alloc( jpi,jpj,jpk, zhpi, zhpj )338 CALL wrk_alloc( jpi,jpj,jpk, zhpi, zhpj ) 288 339 ! 289 340 IF( kt == nit000 ) THEN … … 293 344 ENDIF 294 345 346 ! Partial steps: bottom before horizontal gradient of t, s, rd at the last ocean level 347 !jc CALL zps_hde ( kt, jpts, tsn, gtsu, gtsv, rhd, gru , grv ) 295 348 296 349 ! Local constant initialization … … 300 353 DO jj = 2, jpjm1 301 354 DO ji = fs_2, fs_jpim1 ! vector opt. 302 zcoef1 = zcoef0 * fse3w(ji,jj,1)355 zcoef1 = zcoef0 * e3w_n(ji,jj,1) 303 356 ! hydrostatic pressure gradient 304 zhpi(ji,jj,1) = zcoef1 * ( rhd(ji+1,jj ,1) - rhd(ji,jj,1) ) /e1u(ji,jj)305 zhpj(ji,jj,1) = zcoef1 * ( rhd(ji ,jj+1,1) - rhd(ji,jj,1) ) /e2v(ji,jj)357 zhpi(ji,jj,1) = zcoef1 * ( rhd(ji+1,jj ,1) - rhd(ji,jj,1) ) * r1_e1u(ji,jj) 358 zhpj(ji,jj,1) = zcoef1 * ( rhd(ji ,jj+1,1) - rhd(ji,jj,1) ) * r1_e2v(ji,jj) 306 359 ! add to the general momentum trend 307 360 ua(ji,jj,1) = ua(ji,jj,1) + zhpi(ji,jj,1) … … 309 362 END DO 310 363 END DO 311 312 364 313 365 ! interior value (2=<jk=<jpkm1) … … 315 367 DO jj = 2, jpjm1 316 368 DO ji = fs_2, fs_jpim1 ! vector opt. 317 zcoef1 = zcoef0 * fse3w(ji,jj,jk)369 zcoef1 = zcoef0 * e3w_n(ji,jj,jk) 318 370 ! hydrostatic pressure gradient 319 371 zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) & 320 372 & + zcoef1 * ( ( rhd(ji+1,jj,jk) + rhd(ji+1,jj,jk-1) ) & 321 & - ( rhd(ji ,jj,jk) + rhd(ji ,jj,jk-1) ) ) /e1u(ji,jj)373 & - ( rhd(ji ,jj,jk) + rhd(ji ,jj,jk-1) ) ) * r1_e1u(ji,jj) 322 374 323 375 zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1) & 324 376 & + zcoef1 * ( ( rhd(ji,jj+1,jk) + rhd(ji,jj+1,jk-1) ) & 325 & - ( rhd(ji,jj, jk) + rhd(ji,jj ,jk-1) ) ) /e2v(ji,jj)377 & - ( rhd(ji,jj, jk) + rhd(ji,jj ,jk-1) ) ) * r1_e2v(ji,jj) 326 378 ! add to the general momentum trend 327 379 ua(ji,jj,jk) = ua(ji,jj,jk) + zhpi(ji,jj,jk) … … 330 382 END DO 331 383 END DO 332 333 384 334 385 ! partial steps correction at the last level (use gru & grv computed in zpshde.F90) … … 337 388 iku = mbku(ji,jj) 338 389 ikv = mbkv(ji,jj) 339 zcoef2 = zcoef0 * MIN( fse3w(ji,jj,iku), fse3w(ji+1,jj ,iku) )340 zcoef3 = zcoef0 * MIN( fse3w(ji,jj,ikv), fse3w(ji ,jj+1,ikv) )390 zcoef2 = zcoef0 * MIN( e3w_n(ji,jj,iku), e3w_n(ji+1,jj ,iku) ) 391 zcoef3 = zcoef0 * MIN( e3w_n(ji,jj,ikv), e3w_n(ji ,jj+1,ikv) ) 341 392 IF( iku > 1 ) THEN ! on i-direction (level 2 or more) 342 393 ua (ji,jj,iku) = ua(ji,jj,iku) - zhpi(ji,jj,iku) ! subtract old value 343 394 zhpi(ji,jj,iku) = zhpi(ji,jj,iku-1) & ! compute the new one 344 & + zcoef2 * ( rhd(ji+1,jj,iku-1) - rhd(ji,jj,iku-1) + gru(ji,jj) ) /e1u(ji,jj)395 & + zcoef2 * ( rhd(ji+1,jj,iku-1) - rhd(ji,jj,iku-1) + gru(ji,jj) ) * r1_e1u(ji,jj) 345 396 ua (ji,jj,iku) = ua(ji,jj,iku) + zhpi(ji,jj,iku) ! add the new one to the general momentum trend 346 397 ENDIF … … 348 399 va (ji,jj,ikv) = va(ji,jj,ikv) - zhpj(ji,jj,ikv) ! subtract old value 349 400 zhpj(ji,jj,ikv) = zhpj(ji,jj,ikv-1) & ! compute the new one 350 & + zcoef3 * ( rhd(ji,jj+1,ikv-1) - rhd(ji,jj,ikv-1) + grv(ji,jj) ) /e2v(ji,jj)401 & + zcoef3 * ( rhd(ji,jj+1,ikv-1) - rhd(ji,jj,ikv-1) + grv(ji,jj) ) * r1_e2v(ji,jj) 351 402 va (ji,jj,ikv) = va(ji,jj,ikv) + zhpj(ji,jj,ikv) ! add the new one to the general momentum trend 352 403 ENDIF … … 354 405 END DO 355 406 ! 356 CALL wrk_dealloc( jpi,jpj,jpk, zhpi, zhpj )407 CALL wrk_dealloc( jpi,jpj,jpk, zhpi, zhpj ) 357 408 ! 358 409 END SUBROUTINE hpg_zps 410 359 411 360 412 SUBROUTINE hpg_sco( kt ) … … 378 430 INTEGER, INTENT(in) :: kt ! ocean time-step index 379 431 !! 380 INTEGER :: ji, jj, jk ! dummy loop indices 381 REAL(wp) :: zcoef0, zuap, zvap, znad ! temporary scalars 432 INTEGER :: ji, jj, jk, jii, jjj ! dummy loop indices 433 REAL(wp) :: zcoef0, zuap, zvap, znad, ztmp ! temporary scalars 434 LOGICAL :: ll_tmp1, ll_tmp2, ll_tmp3 ! local logical variables 382 435 REAL(wp), POINTER, DIMENSION(:,:,:) :: zhpi, zhpj 436 REAL(wp), POINTER, DIMENSION(:,:) :: zcpx, zcpy !W/D pressure filter 383 437 !!---------------------------------------------------------------------- 384 438 ! 385 439 CALL wrk_alloc( jpi,jpj,jpk, zhpi, zhpj ) 440 IF(ln_wd) CALL wrk_alloc( jpi,jpj, zcpx, zcpy ) 386 441 ! 387 442 IF( kt == nit000 ) THEN … … 390 445 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ s-coordinate case, OPA original scheme used' 391 446 ENDIF 392 393 ! Local constant initialization 447 ! 394 448 zcoef0 = - grav * 0.5_wp 395 ! To use density and not density anomaly 396 IF ( lk_vvl ) THEN ; znad = 1._wp ! Variable volume 397 ELSE ; znad = 0._wp ! Fixed volume 449 IF ( ln_linssh ) THEN ; znad = 0._wp ! Fixed volume: density anomaly 450 ELSE ; znad = 1._wp ! Variable volume: density 398 451 ENDIF 452 ! 453 IF(ln_wd) THEN 454 DO jj = 2, jpjm1 455 DO ji = 2, jpim1 456 ll_tmp1 = MIN(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj)) 457 ll_tmp2 = MAX(sshn(ji,jj) + bathy(ji,jj), sshn(ji+1,jj) + bathy(ji+1,jj)) > rn_wdmin1 + rn_wdmin2 458 ll_tmp3 = MAX(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj)) + & 459 & rn_wdmin1 + rn_wdmin2 460 461 IF(ll_tmp1.AND.ll_tmp2) THEN 462 zcpx(ji,jj) = 1.0_wp 463 wduflt(ji,jj) = 1.0_wp 464 ELSE IF(ll_tmp3) THEN 465 ! no worries about sshn(ji+1,jj)-sshn(ji,jj) = 0, it won't happen ! here 466 zcpx(ji,jj) = ABS((sshn(ji+1,jj) + bathy(ji+1,jj) - sshn(ji,jj) - bathy(ji,jj)) / & 467 & (sshn(ji+1,jj) - sshn(ji,jj))) 468 wduflt(ji,jj) = 1.0_wp 469 ELSE 470 zcpx(ji,jj) = 0._wp 471 wduflt(ji,jj) = 0.0_wp 472 END IF 473 474 ll_tmp1 = MIN(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1)) 475 ll_tmp2 = MAX(sshn(ji,jj) + bathy(ji,jj), sshn(ji,jj+1) + bathy(ji,jj+1)) > rn_wdmin1 + rn_wdmin2 476 ll_tmp3 = MAX(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1)) + & 477 & rn_wdmin1 + rn_wdmin2 478 479 IF(ll_tmp1.AND.ll_tmp2) THEN 480 zcpy(ji,jj) = 1.0_wp 481 wdvflt(ji,jj) = 1.0_wp 482 ELSE IF(ll_tmp3) THEN 483 ! no worries about sshn(ji,jj+1)-sshn(ji,jj) = 0, it won't happen ! here 484 zcpy(ji,jj) = ABS((sshn(ji,jj+1) + bathy(ji,jj+1) - sshn(ji,jj) - bathy(ji,jj)) / & 485 & (sshn(ji,jj+1) - sshn(ji,jj))) 486 wdvflt(ji,jj) = 1.0_wp 487 ELSE 488 zcpy(ji,jj) = 0._wp 489 wdvflt(ji,jj) = 0.0_wp 490 END IF 491 END DO 492 END DO 493 CALL lbc_lnk( zcpx, 'U', 1._wp ) ; CALL lbc_lnk( zcpy, 'V', 1._wp ) 494 ENDIF 495 399 496 400 497 ! Surface value … … 402 499 DO ji = fs_2, fs_jpim1 ! vector opt. 403 500 ! hydrostatic pressure gradient along s-surfaces 404 zhpi(ji,jj,1) = zcoef0 / e1u(ji,jj) * ( fse3w(ji+1,jj ,1) * ( znad + rhd(ji+1,jj ,1) )&405 & - fse3w(ji ,jj ,1) * ( znad + rhd(ji ,jj ,1) ))406 zhpj(ji,jj,1) = zcoef0 / e2v(ji,jj) * ( fse3w(ji ,jj+1,1) * ( znad + rhd(ji ,jj+1,1) )&407 & - fse3w(ji ,jj ,1) * ( znad + rhd(ji ,jj ,1) ))501 zhpi(ji,jj,1) = zcoef0 * ( e3w_n(ji+1,jj ,1) * ( znad + rhd(ji+1,jj ,1) ) & 502 & - e3w_n(ji ,jj ,1) * ( znad + rhd(ji ,jj ,1) ) ) * r1_e1u(ji,jj) 503 zhpj(ji,jj,1) = zcoef0 * ( e3w_n(ji ,jj+1,1) * ( znad + rhd(ji ,jj+1,1) ) & 504 & - e3w_n(ji ,jj ,1) * ( znad + rhd(ji ,jj ,1) ) ) * r1_e2v(ji,jj) 408 505 ! s-coordinate pressure gradient correction 409 zuap = -zcoef0 * ( rhd (ji+1,jj,1) + rhd (ji,jj,1) + 2._wp * znad ) & 410 & * ( fsde3w(ji+1,jj,1) - fsde3w(ji,jj,1) ) / e1u(ji,jj) 411 zvap = -zcoef0 * ( rhd (ji,jj+1,1) + rhd (ji,jj,1) + 2._wp * znad ) & 412 & * ( fsde3w(ji,jj+1,1) - fsde3w(ji,jj,1) ) / e2v(ji,jj) 506 zuap = -zcoef0 * ( rhd (ji+1,jj,1) + rhd (ji,jj,1) + 2._wp * znad ) & 507 & * ( gde3w_n(ji+1,jj,1) - gde3w_n(ji,jj,1) ) * r1_e1u(ji,jj) 508 zvap = -zcoef0 * ( rhd (ji,jj+1,1) + rhd (ji,jj,1) + 2._wp * znad ) & 509 & * ( gde3w_n(ji,jj+1,1) - gde3w_n(ji,jj,1) ) * r1_e2v(ji,jj) 510 511 512 IF(ln_wd) THEN 513 514 zhpi(ji,jj,1) = zhpi(ji,jj,1) * zcpx(ji,jj) 515 zhpj(ji,jj,1) = zhpj(ji,jj,1) * zcpy(ji,jj) 516 zuap = zuap * zcpx(ji,jj) 517 zvap = zvap * zcpy(ji,jj) 518 ENDIF 519 413 520 ! add to the general momentum trend 414 521 ua(ji,jj,1) = ua(ji,jj,1) + zhpi(ji,jj,1) + zuap … … 422 529 DO ji = fs_2, fs_jpim1 ! vector opt. 423 530 ! hydrostatic pressure gradient along s-surfaces 424 zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) + zcoef0 /e1u(ji,jj) &425 & * ( fse3w(ji+1,jj,jk) * ( rhd(ji+1,jj,jk) + rhd(ji+1,jj,jk-1) + 2*znad ) &426 & - fse3w(ji ,jj,jk) * ( rhd(ji ,jj,jk) + rhd(ji ,jj,jk-1) + 2*znad ) )427 zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1) + zcoef0 /e2v(ji,jj) &428 & * ( fse3w(ji,jj+1,jk) * ( rhd(ji,jj+1,jk) + rhd(ji,jj+1,jk-1) + 2*znad ) &429 & - fse3w(ji,jj ,jk) * ( rhd(ji,jj, jk) + rhd(ji,jj ,jk-1) + 2*znad ) )531 zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) + zcoef0 * r1_e1u(ji,jj) & 532 & * ( e3w_n(ji+1,jj,jk) * ( rhd(ji+1,jj,jk) + rhd(ji+1,jj,jk-1) + 2*znad ) & 533 & - e3w_n(ji ,jj,jk) * ( rhd(ji ,jj,jk) + rhd(ji ,jj,jk-1) + 2*znad ) ) 534 zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1) + zcoef0 * r1_e2v(ji,jj) & 535 & * ( e3w_n(ji,jj+1,jk) * ( rhd(ji,jj+1,jk) + rhd(ji,jj+1,jk-1) + 2*znad ) & 536 & - e3w_n(ji,jj ,jk) * ( rhd(ji,jj, jk) + rhd(ji,jj ,jk-1) + 2*znad ) ) 430 537 ! s-coordinate pressure gradient correction 431 zuap = -zcoef0 * ( rhd (ji+1,jj ,jk) + rhd (ji,jj,jk) + 2._wp * znad ) & 432 & * ( fsde3w(ji+1,jj ,jk) - fsde3w(ji,jj,jk) ) / e1u(ji,jj) 433 zvap = -zcoef0 * ( rhd (ji ,jj+1,jk) + rhd (ji,jj,jk) + 2._wp * znad ) & 434 & * ( fsde3w(ji ,jj+1,jk) - fsde3w(ji,jj,jk) ) / e2v(ji,jj) 538 zuap = -zcoef0 * ( rhd (ji+1,jj ,jk) + rhd (ji,jj,jk) + 2._wp * znad ) & 539 & * ( gde3w_n(ji+1,jj ,jk) - gde3w_n(ji,jj,jk) ) * r1_e1u(ji,jj) 540 zvap = -zcoef0 * ( rhd (ji ,jj+1,jk) + rhd (ji,jj,jk) + 2._wp * znad ) & 541 & * ( gde3w_n(ji ,jj+1,jk) - gde3w_n(ji,jj,jk) ) * r1_e2v(ji,jj) 542 543 IF(ln_wd) THEN 544 zhpi(ji,jj,jk) = zhpi(ji,jj,jk) * zcpx(ji,jj) 545 zhpj(ji,jj,jk) = zhpj(ji,jj,jk) * zcpy(ji,jj) 546 zuap = zuap * zcpx(ji,jj) 547 zvap = zvap * zcpy(ji,jj) 548 ENDIF 549 435 550 ! add to the general momentum trend 436 551 ua(ji,jj,jk) = ua(ji,jj,jk) + zhpi(ji,jj,jk) + zuap … … 441 556 ! 442 557 CALL wrk_dealloc( jpi,jpj,jpk, zhpi, zhpj ) 558 IF(ln_wd) CALL wrk_dealloc( jpi,jpj, zcpx, zcpy ) 443 559 ! 444 560 END SUBROUTINE hpg_sco 561 445 562 446 563 SUBROUTINE hpg_isf( kt ) 447 564 !!--------------------------------------------------------------------- 448 !! *** ROUTINE hpg_ sco***565 !! *** ROUTINE hpg_isf *** 449 566 !! 450 567 !! ** Method : s-coordinate case. Jacobian scheme. … … 465 582 INTEGER, INTENT(in) :: kt ! ocean time-step index 466 583 !! 467 INTEGER :: ji, jj, jk, ik u, ikv, ikt, iktp1i, iktp1j! dummy loop indices468 REAL(wp) :: zcoef0, zuap, zvap, znad , ze3wu, ze3wv, zuapint, zvapint, zhpjint, zhpiint, zdzwt, zdzwtjp1, zdzwtip1! temporary scalars469 REAL(wp), POINTER, DIMENSION(:,:,:) :: zhpi, zhpj , zrhd584 INTEGER :: ji, jj, jk, ikt, iktp1i, iktp1j ! dummy loop indices 585 REAL(wp) :: zcoef0, zuap, zvap, znad ! temporary scalars 586 REAL(wp), POINTER, DIMENSION(:,:,:) :: zhpi, zhpj 470 587 REAL(wp), POINTER, DIMENSION(:,:,:) :: ztstop 471 REAL(wp), POINTER, DIMENSION(:,:) :: ze3w, zp, zrhdtop_isf, zrhdtop_oce, ziceload, zdept, zpshpi, zpshpj 472 !!---------------------------------------------------------------------- 473 ! 474 CALL wrk_alloc( jpi,jpj, 2, ztstop) 475 CALL wrk_alloc( jpi,jpj,jpk, zhpi, zhpj, zrhd) 476 CALL wrk_alloc( jpi,jpj, ze3w, zp, zrhdtop_isf, zrhdtop_oce, ziceload, zdept, zpshpi, zpshpj) 477 ! 478 IF( kt == nit000 ) THEN 479 IF(lwp) WRITE(numout,*) 480 IF(lwp) WRITE(numout,*) 'dyn:hpg_isf : hydrostatic pressure gradient trend for ice shelf' 481 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ s-coordinate case, OPA original scheme used' 482 ENDIF 483 588 REAL(wp), POINTER, DIMENSION(:,:) :: zrhdtop_oce 589 !!---------------------------------------------------------------------- 590 ! 591 CALL wrk_alloc( jpi,jpj, 2, ztstop) 592 CALL wrk_alloc( jpi,jpj,jpk, zhpi, zhpj) 593 CALL wrk_alloc( jpi,jpj, zrhdtop_oce ) 594 ! 484 595 ! Local constant initialization 485 596 zcoef0 = - grav * 0.5_wp 597 486 598 ! To use density and not density anomaly 487 ! IF ( lk_vvl ) THEN ; znad = 1._wp ! Variable volume488 ! ELSE ; znad = 0._wp ! Fixed volume489 ! ENDIF490 599 znad=1._wp 600 491 601 ! iniitialised to 0. zhpi zhpi 492 602 zhpi(:,:,:)=0._wp ; zhpj(:,:,:)=0._wp 493 603 494 !==================================================================================495 !=====Compute iceload and contribution of the half first wet layer =================496 !===================================================================================497 498 ! assume water displaced by the ice shelf is at T=-1.9 and S=34.4 (rude)499 ztstop(:,:,1)=-1.9_wp ; ztstop(:,:,2)=34.4_wp500 501 ! compute density of the water displaced by the ice shelf502 zrhd = rhd ! save rhd503 DO jk = 1, jpk504 zdept(:,:)=gdept_1d(jk)505 CALL eos(ztstop(:,:,:),zdept(:,:),rhd(:,:,jk))506 END DO507 WHERE ( tmask(:,:,:) == 1._wp)508 rhd(:,:,:) = zrhd(:,:,:) ! replace wet cell by the saved rhd509 END WHERE510 511 ! compute rhd at the ice/oce interface (ice shelf side)512 CALL eos(ztstop,risfdep,zrhdtop_isf)513 514 604 ! compute rhd at the ice/oce interface (ocean side) 605 ! usefull to reduce residual current in the test case ISOMIP with no melting 515 606 DO ji=1,jpi 516 607 DO jj=1,jpj … … 520 611 END DO 521 612 END DO 522 CALL eos(ztstop,risfdep,zrhdtop_oce) 523 ! 524 ! Surface value + ice shelf gradient 525 ! compute pressure due to ice shelf load (used to compute hpgi/j for all the level from 1 to miku/v) 526 ziceload = 0._wp 527 DO jj = 1, jpj 528 DO ji = 1, jpi ! vector opt. 529 ikt=mikt(ji,jj) 530 ziceload(ji,jj) = ziceload(ji,jj) + (znad + rhd(ji,jj,1) ) * fse3w(ji,jj,1) * (1._wp - tmask(ji,jj,1)) 531 DO jk=2,ikt-1 532 ziceload(ji,jj) = ziceload(ji,jj) + (2._wp * znad + rhd(ji,jj,jk-1) + rhd(ji,jj,jk)) * fse3w(ji,jj,jk) & 533 & * (1._wp - tmask(ji,jj,jk)) 534 END DO 535 IF (ikt .GE. 2) ziceload(ji,jj) = ziceload(ji,jj) + (2._wp * znad + zrhdtop_isf(ji,jj) + rhd(ji,jj,ikt-1)) & 536 & * ( risfdep(ji,jj) - gdept_1d(ikt-1) ) 537 END DO 538 END DO 539 riceload(:,:) = 0.0_wp ; riceload(:,:)=ziceload(:,:) ! need to be saved for diaar5 540 ! compute zp from z=0 to first T wet point (correction due to zps not yet applied) 613 CALL eos( ztstop, risfdep, zrhdtop_oce ) 614 615 !================================================================================== 616 !===== Compute surface value ===================================================== 617 !================================================================================== 541 618 DO jj = 2, jpjm1 542 619 DO ji = fs_2, fs_jpim1 ! vector opt. 543 ikt=mikt(ji,jj) ; iktp1i=mikt(ji+1,jj); iktp1j=mikt(ji,jj+1) 620 ikt = mikt(ji,jj) 621 iktp1i = mikt(ji+1,jj) 622 iktp1j = mikt(ji,jj+1) 544 623 ! hydrostatic pressure gradient along s-surfaces and ice shelf pressure 545 624 ! we assume ISF is in isostatic equilibrium 546 zhpi(ji,jj,1) = zcoef0 / e1u(ji,jj) * ( 0.5_wp * fse3w(ji+1,jj,iktp1i) &547 & * ( 2._wp * znad + rhd(ji+1,jj ,iktp1i) + zrhdtop_oce(ji+1,jj) ) &548 & - 0.5_wp * fse3w(ji ,jj ,ikt )&549 & * ( 2._wp * znad + rhd(ji ,jj ,ikt ) + zrhdtop_oce(ji ,jj ) )&550 & + ( ziceload(ji+1,jj) - ziceload(ji,jj)))551 zhpj(ji,jj,1) = zcoef0 / e2v(ji,jj) * ( 0.5_wp * fse3w(ji,jj+1,iktp1j) &552 & * ( 2._wp * znad + rhd(ji ,jj+1,iktp1j) + zrhdtop_oce(ji,jj+1) ) &553 & - 0.5_wp * fse3w(ji ,jj ,ikt )&554 & * ( 2._wp * znad + rhd(ji ,jj ,ikt ) + zrhdtop_oce(ji ,jj ) )&555 & + ( ziceload(ji,jj+1) - ziceload(ji,jj) ))625 zhpi(ji,jj,1) = zcoef0 / e1u(ji,jj) * ( 0.5_wp * e3w_n(ji+1,jj,iktp1i) & 626 & * ( 2._wp * znad + rhd(ji+1,jj,iktp1i) + zrhdtop_oce(ji+1,jj) ) & 627 & - 0.5_wp * e3w_n(ji,jj,ikt) & 628 & * ( 2._wp * znad + rhd(ji,jj,ikt) + zrhdtop_oce(ji,jj) ) & 629 & + ( riceload(ji+1,jj) - riceload(ji,jj)) ) 630 zhpj(ji,jj,1) = zcoef0 / e2v(ji,jj) * ( 0.5_wp * e3w_n(ji,jj+1,iktp1j) & 631 & * ( 2._wp * znad + rhd(ji,jj+1,iktp1j) + zrhdtop_oce(ji,jj+1) ) & 632 & - 0.5_wp * e3w_n(ji,jj,ikt) & 633 & * ( 2._wp * znad + rhd(ji,jj,ikt) + zrhdtop_oce(ji,jj) ) & 634 & + ( riceload(ji,jj+1) - riceload(ji,jj)) ) 556 635 ! s-coordinate pressure gradient correction (=0 if z coordinate) 557 zuap = -zcoef0 * ( rhd (ji+1,jj,1) + rhd(ji,jj,1) + 2._wp * znad ) &558 & * ( fsde3w(ji+1,jj,1) - fsde3w(ji,jj,1) ) /e1u(ji,jj)559 zvap = -zcoef0 * ( rhd (ji,jj+1,1) + rhd(ji,jj,1) + 2._wp * znad ) &560 & * ( fsde3w(ji,jj+1,1) - fsde3w(ji,jj,1) ) /e2v(ji,jj)636 zuap = -zcoef0 * ( rhd (ji+1,jj,1) + rhd (ji,jj,1) + 2._wp * znad ) & 637 & * ( gde3w_n(ji+1,jj,1) - gde3w_n(ji,jj,1) ) * r1_e1u(ji,jj) 638 zvap = -zcoef0 * ( rhd (ji,jj+1,1) + rhd (ji,jj,1) + 2._wp * znad ) & 639 & * ( gde3w_n(ji,jj+1,1) - gde3w_n(ji,jj,1) ) * r1_e2v(ji,jj) 561 640 ! add to the general momentum trend 562 641 ua(ji,jj,1) = ua(ji,jj,1) + (zhpi(ji,jj,1) + zuap) * umask(ji,jj,1) … … 565 644 END DO 566 645 !================================================================================== 567 !===== Compute partial cell contribution for the top cell =========================568 !==================================================================================569 DO jj = 2, jpjm1570 DO ji = fs_2, fs_jpim1 ! vector opt.571 iku = miku(ji,jj) ;572 zpshpi(ji,jj)=0.0_wp ; zpshpj(ji,jj)=0.0_wp573 ze3wu = (gdepw_0(ji+1,jj,iku+1) - gdept_0(ji+1,jj,iku)) - (gdepw_0(ji,jj,iku+1) - gdept_0(ji,jj,iku))574 ! u direction575 IF ( iku .GT. 1 ) THEN576 ! case iku577 zhpi(ji,jj,iku) = zcoef0 / e1u(ji,jj) * ze3wu &578 & * ( rhd (ji+1,jj,iku) + rhd (ji,jj,iku) &579 & + SIGN(1._wp,ze3wu) * grui(ji,jj) + 2._wp * znad )580 ! corrective term ( = 0 if z coordinate )581 zuap = -zcoef0 * ( arui(ji,jj) + 2._wp * znad ) * gzui(ji,jj) / e1u(ji,jj)582 ! zhpi will be added in interior loop583 ua(ji,jj,iku) = ua(ji,jj,iku) + zuap584 ! in case of 2 cell water column, need to save the pressure gradient to compute the bottom pressure585 IF (mbku(ji,jj) == iku + 1) zpshpi(ji,jj) = zhpi(ji,jj,iku)586 587 ! case iku + 1 (remove the zphi term added in the interior loop and compute the one corrected for zps)588 zhpiint = zcoef0 / e1u(ji,jj) &589 & * ( fse3w(ji+1,jj ,iku+1) * ( (rhd(ji+1,jj,iku+1) + znad) &590 & + (rhd(ji+1,jj,iku ) + znad) ) * tmask(ji+1,jj,iku) &591 & - fse3w(ji ,jj ,iku+1) * ( (rhd(ji ,jj,iku+1) + znad) &592 & + (rhd(ji ,jj,iku ) + znad) ) * tmask(ji ,jj,iku) )593 zhpi(ji,jj,iku+1) = zcoef0 / e1u(ji,jj) * ge3rui(ji,jj) - zhpiint594 END IF595 596 ! v direction597 ikv = mikv(ji,jj)598 ze3wv = (gdepw_0(ji,jj+1,ikv+1) - gdept_0(ji,jj+1,ikv)) - (gdepw_0(ji,jj,ikv+1) - gdept_0(ji,jj,ikv))599 IF ( ikv .GT. 1 ) THEN600 ! case ikv601 zhpj(ji,jj,ikv) = zcoef0 / e2v(ji,jj) * ze3wv &602 & * ( rhd(ji,jj+1,ikv) + rhd (ji,jj,ikv) &603 & + SIGN(1._wp,ze3wv) * grvi(ji,jj) + 2._wp * znad )604 ! corrective term ( = 0 if z coordinate )605 zvap = -zcoef0 * ( arvi(ji,jj) + 2._wp * znad ) * gzvi(ji,jj) / e2v(ji,jj)606 ! zhpi will be added in interior loop607 va(ji,jj,ikv) = va(ji,jj,ikv) + zvap608 ! in case of 2 cell water column, need to save the pressure gradient to compute the bottom pressure609 IF (mbkv(ji,jj) == ikv + 1) zpshpj(ji,jj) = zhpj(ji,jj,ikv)610 611 ! case ikv + 1 (remove the zphj term added in the interior loop and compute the one corrected for zps)612 zhpjint = zcoef0 / e2v(ji,jj) &613 & * ( fse3w(ji ,jj+1,ikv+1) * ( (rhd(ji,jj+1,ikv+1) + znad) &614 & + (rhd(ji,jj+1,ikv ) + znad) ) * tmask(ji,jj+1,ikv) &615 & - fse3w(ji ,jj ,ikv+1) * ( (rhd(ji,jj ,ikv+1) + znad) &616 & + (rhd(ji,jj ,ikv ) + znad) ) * tmask(ji,jj ,ikv) )617 zhpj(ji,jj,ikv+1) = zcoef0 / e2v(ji,jj) * ge3rvi(ji,jj) - zhpjint618 END IF619 END DO620 END DO621 622 !==================================================================================623 646 !===== Compute interior value ===================================================== 624 647 !================================================================================== 625 626 DO jj = 2, jpjm1 627 DO ji = fs_2, fs_jpim1 ! vector opt. 628 iku=miku(ji,jj); ikv=mikv(ji,jj) 629 DO jk = 2, jpkm1 648 ! interior value (2=<jk=<jpkm1) 649 DO jk = 2, jpkm1 650 DO jj = 2, jpjm1 651 DO ji = fs_2, fs_jpim1 ! vector opt. 630 652 ! hydrostatic pressure gradient along s-surfaces 631 ! zhpi is masked for the first wet cell (contribution already done in the upper bloc) 632 zhpi(ji,jj,jk) = zhpi(ji,jj,jk) + zhpi(ji,jj,jk-1) & 633 & + zcoef0 / e1u(ji,jj) & 634 & * ( fse3w(ji+1,jj ,jk) * ( (rhd(ji+1,jj,jk ) + znad) & 635 & + (rhd(ji+1,jj,jk-1) + znad) ) * tmask(ji+1,jj,jk-1) & 636 & - fse3w(ji ,jj ,jk) * ( (rhd(ji ,jj,jk ) + znad) & 637 & + (rhd(ji ,jj,jk-1) + znad) ) * tmask(ji ,jj,jk-1) ) 653 zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) + zcoef0 / e1u(ji,jj) & 654 & * ( e3w_n(ji+1,jj,jk) * ( rhd(ji+1,jj,jk) + rhd(ji+1,jj,jk-1) + 2*znad ) * wmask(ji+1,jj,jk) & 655 & - e3w_n(ji ,jj,jk) * ( rhd(ji ,jj,jk) + rhd(ji ,jj,jk-1) + 2*znad ) * wmask(ji ,jj,jk) ) 656 zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1) + zcoef0 / e2v(ji,jj) & 657 & * ( e3w_n(ji,jj+1,jk) * ( rhd(ji,jj+1,jk) + rhd(ji,jj+1,jk-1) + 2*znad ) * wmask(ji,jj+1,jk) & 658 & - e3w_n(ji,jj ,jk) * ( rhd(ji,jj, jk) + rhd(ji,jj ,jk-1) + 2*znad ) * wmask(ji,jj ,jk) ) 638 659 ! s-coordinate pressure gradient correction 639 ! corrective term, we mask this term for the first wet level beneath the ice shelf (contribution done in the upper bloc) 640 zuap = - zcoef0 * ( rhd (ji+1,jj ,jk) + rhd (ji,jj,jk) + 2._wp * znad ) & 641 & * ( fsde3w(ji+1,jj ,jk) - fsde3w(ji,jj,jk) ) / e1u(ji,jj) * umask(ji,jj,jk-1) 642 ua(ji,jj,jk) = ua(ji,jj,jk) + ( zhpi(ji,jj,jk) + zuap) * umask(ji,jj,jk) 643 644 ! hydrostatic pressure gradient along s-surfaces 645 ! zhpi is masked for the first wet cell (contribution already done in the upper bloc) 646 zhpj(ji,jj,jk) = zhpj(ji,jj,jk) + zhpj(ji,jj,jk-1) & 647 & + zcoef0 / e2v(ji,jj) & 648 & * ( fse3w(ji ,jj+1,jk) * ( (rhd(ji,jj+1,jk ) + znad) & 649 & + (rhd(ji,jj+1,jk-1) + znad) ) * tmask(ji,jj+1,jk-1) & 650 & - fse3w(ji ,jj ,jk) * ( (rhd(ji,jj ,jk ) + znad) & 651 & + (rhd(ji,jj ,jk-1) + znad) ) * tmask(ji,jj ,jk-1) ) 652 ! s-coordinate pressure gradient correction 653 ! corrective term, we mask this term for the first wet level beneath the ice shelf (contribution done in the upper bloc) 654 zvap = - zcoef0 * ( rhd (ji ,jj+1,jk) + rhd (ji,jj,jk) + 2._wp * znad ) & 655 & * ( fsde3w(ji ,jj+1,jk) - fsde3w(ji,jj,jk) ) / e2v(ji,jj) * vmask(ji,jj,jk-1) 660 zuap = -zcoef0 * ( rhd (ji+1,jj ,jk) + rhd (ji,jj,jk) + 2._wp * znad ) & 661 & * ( gde3w_n(ji+1,jj ,jk) - gde3w_n(ji,jj,jk) ) / e1u(ji,jj) 662 zvap = -zcoef0 * ( rhd (ji ,jj+1,jk) + rhd (ji,jj,jk) + 2._wp * znad ) & 663 & * ( gde3w_n(ji ,jj+1,jk) - gde3w_n(ji,jj,jk) ) / e2v(ji,jj) 656 664 ! add to the general momentum trend 657 va(ji,jj,jk) = va(ji,jj,jk) + ( zhpj(ji,jj,jk) + zvap ) * vmask(ji,jj,jk) 665 ua(ji,jj,jk) = ua(ji,jj,jk) + (zhpi(ji,jj,jk) + zuap) * umask(ji,jj,jk) 666 va(ji,jj,jk) = va(ji,jj,jk) + (zhpj(ji,jj,jk) + zvap) * vmask(ji,jj,jk) 658 667 END DO 659 668 END DO 660 669 END DO 661 662 !================================================================================== 663 !===== Compute bottom cell contribution (partial cell) ============================ 664 !================================================================================== 665 666 DO jj = 2, jpjm1 667 DO ji = 2, jpim1 668 iku = mbku(ji,jj) 669 ikv = mbkv(ji,jj) 670 671 IF (iku .GT. 1) THEN 672 ! remove old value (interior case) 673 zuap = -zcoef0 * ( rhd (ji+1,jj ,iku) + rhd (ji,jj,iku) + 2._wp * znad ) & 674 & * ( fsde3w(ji+1,jj ,iku) - fsde3w(ji,jj,iku) ) / e1u(ji,jj) 675 ua(ji,jj,iku) = ua(ji,jj,iku) - zhpi(ji,jj,iku) - zuap 676 ! put new value 677 ! -zpshpi to avoid double contribution of the partial step in the top layer 678 zuap = -zcoef0 * ( aru(ji,jj) + 2._wp * znad ) * gzu(ji,jj) / e1u(ji,jj) 679 zhpi(ji,jj,iku) = zhpi(ji,jj,iku-1) + zcoef0 / e1u(ji,jj) * ge3ru(ji,jj) - zpshpi(ji,jj) 680 ua(ji,jj,iku) = ua(ji,jj,iku) + zhpi(ji,jj,iku) + zuap 681 END IF 682 ! v direction 683 IF (ikv .GT. 1) THEN 684 ! remove old value (interior case) 685 zvap = -zcoef0 * ( rhd (ji ,jj+1,ikv) + rhd (ji,jj,ikv) + 2._wp * znad ) & 686 & * ( fsde3w(ji ,jj+1,ikv) - fsde3w(ji,jj,ikv) ) / e2v(ji,jj) 687 va(ji,jj,ikv) = va(ji,jj,ikv) - zhpj(ji,jj,ikv) - zvap 688 ! put new value 689 ! -zpshpj to avoid double contribution of the partial step in the top layer 690 zvap = -zcoef0 * ( arv(ji,jj) + 2._wp * znad ) * gzv(ji,jj) / e2v(ji,jj) 691 zhpj(ji,jj,ikv) = zhpj(ji,jj,ikv-1) + zcoef0 / e2v(ji,jj) * ge3rv(ji,jj) - zpshpj(ji,jj) 692 va(ji,jj,ikv) = va(ji,jj,ikv) + zhpj(ji,jj,ikv) + zvap 693 END IF 694 END DO 695 END DO 696 697 ! set back to original density value into the ice shelf cell (maybe useless because it is masked) 698 rhd = zrhd 699 ! 700 CALL wrk_dealloc( jpi,jpj,2, ztstop) 701 CALL wrk_dealloc( jpi,jpj,jpk, zhpi, zhpj, zrhd) 702 CALL wrk_dealloc( jpi,jpj, ze3w, zp, zrhdtop_isf, zrhdtop_oce, ziceload, zdept, zpshpi, zpshpj) 670 ! 671 CALL wrk_dealloc( jpi,jpj,2 , ztstop) 672 CALL wrk_dealloc( jpi,jpj,jpk, zhpi, zhpj) 673 CALL wrk_dealloc( jpi,jpj , zrhdtop_oce ) 703 674 ! 704 675 END SUBROUTINE hpg_isf … … 719 690 REAL(wp) :: z1_10, cffu, cffx ! " " 720 691 REAL(wp) :: z1_12, cffv, cffy ! " " 692 LOGICAL :: ll_tmp1, ll_tmp2 ! local logical variables 721 693 REAL(wp), POINTER, DIMENSION(:,:,:) :: zhpi, zhpj 722 694 REAL(wp), POINTER, DIMENSION(:,:,:) :: dzx, dzy, dzz, dzu, dzv, dzw 723 695 REAL(wp), POINTER, DIMENSION(:,:,:) :: drhox, drhoy, drhoz, drhou, drhov, drhow 724 696 REAL(wp), POINTER, DIMENSION(:,:,:) :: rho_i, rho_j, rho_k 697 REAL(wp), POINTER, DIMENSION(:,:) :: zcpx, zcpy !W/D pressure filter 725 698 !!---------------------------------------------------------------------- 726 699 ! … … 728 701 CALL wrk_alloc( jpi, jpj, jpk, drhox, drhoy, drhoz, drhou, drhov, drhow ) 729 702 CALL wrk_alloc( jpi, jpj, jpk, rho_i, rho_j, rho_k, zhpi, zhpj ) 730 ! 703 IF(ln_wd) CALL wrk_alloc( jpi,jpj, zcpx, zcpy ) 704 ! 705 ! 706 IF(ln_wd) THEN 707 DO jj = 2, jpjm1 708 DO ji = 2, jpim1 709 ll_tmp1 = MIN(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj)) & 710 & .and. MAX(sshn(ji,jj) + bathy(ji,jj), sshn(ji+1,jj) + bathy(ji+1,jj)) & 711 & > rn_wdmin1 + rn_wdmin2 712 ll_tmp2 = MAX(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj)) +& 713 & rn_wdmin1 + rn_wdmin2 714 715 IF(ll_tmp1) THEN 716 zcpx(ji,jj) = 1.0_wp 717 ELSE IF(ll_tmp2) THEN 718 ! no worries about sshn(ji+1,jj)-sshn(ji,jj) = 0, it won't happen ! here 719 zcpx(ji,jj) = ABS((sshn(ji+1,jj) + bathy(ji+1,jj) - sshn(ji,jj) - bathy(ji,jj)) /& 720 & (sshn(ji+1,jj) - sshn(ji,jj))) 721 ELSE 722 zcpx(ji,jj) = 0._wp 723 END IF 724 725 ll_tmp1 = MIN(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1)) & 726 & .and. MAX(sshn(ji,jj) + bathy(ji,jj), sshn(ji,jj+1) + bathy(ji,jj+1)) & 727 & > rn_wdmin1 + rn_wdmin2 728 ll_tmp2 = MAX(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1)) +& 729 & rn_wdmin1 + rn_wdmin2 730 731 IF(ll_tmp1) THEN 732 zcpy(ji,jj) = 1.0_wp 733 ELSE IF(ll_tmp2) THEN 734 ! no worries about sshn(ji,jj+1)-sshn(ji,jj) = 0, it won't happen ! here 735 zcpy(ji,jj) = ABS((sshn(ji,jj+1) + bathy(ji,jj+1) - sshn(ji,jj) - bathy(ji,jj)) /& 736 & (sshn(ji,jj+1) - sshn(ji,jj))) 737 ELSE 738 zcpy(ji,jj) = 0._wp 739 END IF 740 END DO 741 END DO 742 CALL lbc_lnk( zcpx, 'U', 1._wp ) ; CALL lbc_lnk( zcpy, 'V', 1._wp ) 743 ENDIF 744 731 745 732 746 IF( kt == nit000 ) THEN … … 750 764 DO jj = 2, jpjm1 751 765 DO ji = fs_2, fs_jpim1 ! vector opt. 752 drhoz(ji,jj,jk) = rhd (ji ,jj ,jk) - rhd(ji,jj,jk-1)753 dzz (ji,jj,jk) = fsde3w(ji ,jj ,jk) - fsde3w(ji,jj,jk-1)754 drhox(ji,jj,jk) = rhd (ji+1,jj ,jk) - rhd(ji,jj,jk )755 dzx (ji,jj,jk) = fsde3w(ji+1,jj ,jk) - fsde3w(ji,jj,jk )756 drhoy(ji,jj,jk) = rhd (ji ,jj+1,jk) - rhd(ji,jj,jk )757 dzy (ji,jj,jk) = fsde3w(ji ,jj+1,jk) - fsde3w(ji,jj,jk )766 drhoz(ji,jj,jk) = rhd (ji ,jj ,jk) - rhd (ji,jj,jk-1) 767 dzz (ji,jj,jk) = gde3w_n(ji ,jj ,jk) - gde3w_n(ji,jj,jk-1) 768 drhox(ji,jj,jk) = rhd (ji+1,jj ,jk) - rhd (ji,jj,jk ) 769 dzx (ji,jj,jk) = gde3w_n(ji+1,jj ,jk) - gde3w_n(ji,jj,jk ) 770 drhoy(ji,jj,jk) = rhd (ji ,jj+1,jk) - rhd (ji,jj,jk ) 771 dzy (ji,jj,jk) = gde3w_n(ji ,jj+1,jk) - gde3w_n(ji,jj,jk ) 758 772 END DO 759 773 END DO … … 837 851 !------------------------------------------------------------- 838 852 839 !!bug gm : e3w- de3w = 0.5*e3w .... and de3w(2)-de3w(1)=e3w(2) .... to be verified840 ! true if de3w is really defined as the sum of the e3w scale factors as, it seems to me, it should be853 !!bug gm : e3w-gde3w = 0.5*e3w .... and gde3w(2)-gde3w(1)=e3w(2) .... to be verified 854 ! true if gde3w is really defined as the sum of the e3w scale factors as, it seems to me, it should be 841 855 842 856 DO jj = 2, jpjm1 843 857 DO ji = fs_2, fs_jpim1 ! vector opt. 844 rho_k(ji,jj,1) = -grav * ( fse3w(ji,jj,1) - fsde3w(ji,jj,1) ) &845 & * ( rhd(ji,jj,1) &846 & + 0.5_wp * ( rhd (ji,jj,2) - rhd(ji,jj,1) )&847 & * ( fse3w (ji,jj,1) - fsde3w(ji,jj,1) )&848 & / ( fsde3w(ji,jj,2) - fsde3w(ji,jj,1) ) )858 rho_k(ji,jj,1) = -grav * ( e3w_n(ji,jj,1) - gde3w_n(ji,jj,1) ) & 859 & * ( rhd(ji,jj,1) & 860 & + 0.5_wp * ( rhd (ji,jj,2) - rhd (ji,jj,1) ) & 861 & * ( e3w_n (ji,jj,1) - gde3w_n(ji,jj,1) ) & 862 & / ( gde3w_n(ji,jj,2) - gde3w_n(ji,jj,1) ) ) 849 863 END DO 850 864 END DO … … 857 871 DO ji = fs_2, fs_jpim1 ! vector opt. 858 872 859 rho_k(ji,jj,jk) = zcoef0 * ( rhd (ji,jj,jk) + rhd(ji,jj,jk-1) ) &860 & * ( fsde3w(ji,jj,jk) - fsde3w(ji,jj,jk-1) ) &861 & - grav * z1_10 * ( 862 & ( drhow (ji,jj,jk) - drhow(ji,jj,jk-1) ) &863 & * ( fsde3w(ji,jj,jk) - fsde3w(ji,jj,jk-1) - z1_12 * ( dzw (ji,jj,jk) + dzw (ji,jj,jk-1) ) ) &864 & - ( dzw (ji,jj,jk) - dzw(ji,jj,jk-1) ) &865 & * ( rhd (ji,jj,jk) - rhd(ji,jj,jk-1) - z1_12 * ( drhow(ji,jj,jk) + drhow(ji,jj,jk-1) ) ) &873 rho_k(ji,jj,jk) = zcoef0 * ( rhd (ji,jj,jk) + rhd (ji,jj,jk-1) ) & 874 & * ( gde3w_n(ji,jj,jk) - gde3w_n(ji,jj,jk-1) ) & 875 & - grav * z1_10 * ( & 876 & ( drhow (ji,jj,jk) - drhow (ji,jj,jk-1) ) & 877 & * ( gde3w_n(ji,jj,jk) - gde3w_n(ji,jj,jk-1) - z1_12 * ( dzw (ji,jj,jk) + dzw (ji,jj,jk-1) ) ) & 878 & - ( dzw (ji,jj,jk) - dzw (ji,jj,jk-1) ) & 879 & * ( rhd (ji,jj,jk) - rhd (ji,jj,jk-1) - z1_12 * ( drhow(ji,jj,jk) + drhow(ji,jj,jk-1) ) ) & 866 880 & ) 867 881 868 rho_i(ji,jj,jk) = zcoef0 * ( rhd (ji+1,jj,jk) + rhd(ji,jj,jk) ) &869 & * ( fsde3w(ji+1,jj,jk) - fsde3w(ji,jj,jk) ) &870 & - grav* z1_10 * ( 871 & ( drhou (ji+1,jj,jk) - drhou(ji,jj,jk) ) &872 & * ( fsde3w(ji+1,jj,jk) - fsde3w(ji,jj,jk) - z1_12 * ( dzu (ji+1,jj,jk) + dzu (ji,jj,jk) ) ) &873 & - ( dzu (ji+1,jj,jk) - dzu(ji,jj,jk) ) &874 & * ( rhd (ji+1,jj,jk) - rhd(ji,jj,jk) - z1_12 * ( drhou(ji+1,jj,jk) + drhou(ji,jj,jk) ) ) &882 rho_i(ji,jj,jk) = zcoef0 * ( rhd (ji+1,jj,jk) + rhd (ji,jj,jk) ) & 883 & * ( gde3w_n(ji+1,jj,jk) - gde3w_n(ji,jj,jk) ) & 884 & - grav* z1_10 * ( & 885 & ( drhou (ji+1,jj,jk) - drhou (ji,jj,jk) ) & 886 & * ( gde3w_n(ji+1,jj,jk) - gde3w_n(ji,jj,jk) - z1_12 * ( dzu (ji+1,jj,jk) + dzu (ji,jj,jk) ) ) & 887 & - ( dzu (ji+1,jj,jk) - dzu (ji,jj,jk) ) & 888 & * ( rhd (ji+1,jj,jk) - rhd (ji,jj,jk) - z1_12 * ( drhou(ji+1,jj,jk) + drhou(ji,jj,jk) ) ) & 875 889 & ) 876 890 877 rho_j(ji,jj,jk) = zcoef0 * ( rhd (ji,jj+1,jk) + rhd (ji,jj,jk) )&878 & * ( fsde3w(ji,jj+1,jk) - fsde3w(ji,jj,jk) ) &879 & - grav* z1_10 * ( 880 & ( drhov (ji,jj+1,jk) - drhov(ji,jj,jk) ) &881 & * ( fsde3w(ji,jj+1,jk) - fsde3w(ji,jj,jk) - z1_12 * ( dzv (ji,jj+1,jk) + dzv (ji,jj,jk) ) ) &882 & - ( dzv (ji,jj+1,jk) - dzv(ji,jj,jk) ) &883 & * ( rhd (ji,jj+1,jk) - rhd(ji,jj,jk) - z1_12 * ( drhov(ji,jj+1,jk) + drhov(ji,jj,jk) ) ) &891 rho_j(ji,jj,jk) = zcoef0 * ( rhd (ji,jj+1,jk) + rhd (ji,jj,jk) ) & 892 & * ( gde3w_n(ji,jj+1,jk) - gde3w_n(ji,jj,jk) ) & 893 & - grav* z1_10 * ( & 894 & ( drhov (ji,jj+1,jk) - drhov (ji,jj,jk) ) & 895 & * ( gde3w_n(ji,jj+1,jk) - gde3w_n(ji,jj,jk) - z1_12 * ( dzv (ji,jj+1,jk) + dzv (ji,jj,jk) ) ) & 896 & - ( dzv (ji,jj+1,jk) - dzv (ji,jj,jk) ) & 897 & * ( rhd (ji,jj+1,jk) - rhd (ji,jj,jk) - z1_12 * ( drhov(ji,jj+1,jk) + drhov(ji,jj,jk) ) ) & 884 898 & ) 885 899 … … 897 911 DO jj = 2, jpjm1 898 912 DO ji = fs_2, fs_jpim1 ! vector opt. 899 zhpi(ji,jj,1) = ( rho_k(ji+1,jj ,1) - rho_k(ji,jj,1) - rho_i(ji,jj,1) ) / e1u(ji,jj) 900 zhpj(ji,jj,1) = ( rho_k(ji ,jj+1,1) - rho_k(ji,jj,1) - rho_j(ji,jj,1) ) / e2v(ji,jj) 913 zhpi(ji,jj,1) = ( rho_k(ji+1,jj ,1) - rho_k(ji,jj,1) - rho_i(ji,jj,1) ) * r1_e1u(ji,jj) 914 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 IF(ln_wd) THEN 916 zhpi(ji,jj,1) = zhpi(ji,jj,1) * zcpx(ji,jj) 917 zhpj(ji,jj,1) = zhpj(ji,jj,1) * zcpy(ji,jj) 918 ENDIF 901 919 ! add to the general momentum trend 902 920 ua(ji,jj,1) = ua(ji,jj,1) + zhpi(ji,jj,1) … … 914 932 zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) & 915 933 & + ( ( rho_k(ji+1,jj,jk) - rho_k(ji,jj,jk ) ) & 916 & - ( rho_i(ji ,jj,jk) - rho_i(ji,jj,jk-1) ) ) /e1u(ji,jj)934 & - ( rho_i(ji ,jj,jk) - rho_i(ji,jj,jk-1) ) ) * r1_e1u(ji,jj) 917 935 zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1) & 918 936 & + ( ( rho_k(ji,jj+1,jk) - rho_k(ji,jj,jk ) ) & 919 & -( rho_j(ji,jj ,jk) - rho_j(ji,jj,jk-1) ) ) / e2v(ji,jj) 937 & -( rho_j(ji,jj ,jk) - rho_j(ji,jj,jk-1) ) ) * r1_e2v(ji,jj) 938 IF(ln_wd) THEN 939 zhpi(ji,jj,jk) = zhpi(ji,jj,jk) * zcpx(ji,jj) 940 zhpj(ji,jj,jk) = zhpj(ji,jj,jk) * zcpy(ji,jj) 941 ENDIF 920 942 ! add to the general momentum trend 921 943 ua(ji,jj,jk) = ua(ji,jj,jk) + zhpi(ji,jj,jk) … … 928 950 CALL wrk_dealloc( jpi, jpj, jpk, drhox, drhoy, drhoz, drhou, drhov, drhow ) 929 951 CALL wrk_dealloc( jpi, jpj, jpk, rho_i, rho_j, rho_k, zhpi, zhpj ) 952 IF(ln_wd) CALL wrk_dealloc( jpi,jpj, zcpx, zcpy ) 930 953 ! 931 954 END SUBROUTINE hpg_djc … … 947 970 !! 948 971 INTEGER :: ji, jj, jk, jkk ! dummy loop indices 949 REAL(wp) :: zcoef0, znad ! temporaryscalars950 ! !972 REAL(wp) :: zcoef0, znad ! local scalars 973 ! 951 974 !! The local variables for the correction term 952 975 INTEGER :: jk1, jis, jid, jjs, jjd 976 LOGICAL :: ll_tmp1, ll_tmp2 ! local logical variables 953 977 REAL(wp) :: zuijk, zvijk, zpwes, zpwed, zpnss, zpnsd, zdeps 954 978 REAL(wp) :: zrhdt1 … … 957 981 REAL(wp), POINTER, DIMENSION(:,:,:) :: zhpi, zu, zv, fsp, xsp, asp, bsp, csp, dsp 958 982 REAL(wp), POINTER, DIMENSION(:,:) :: zsshu_n, zsshv_n 959 !!---------------------------------------------------------------------- 960 ! 961 CALL wrk_alloc( jpi,jpj,jpk, zhpi, zu, zv, fsp, xsp, asp, bsp, csp, dsp ) 962 CALL wrk_alloc( jpi,jpj,jpk, zdept, zrhh ) 963 CALL wrk_alloc( jpi,jpj, zsshu_n, zsshv_n ) 983 REAL(wp), POINTER, DIMENSION(:,:) :: zcpx, zcpy !W/D pressure filter 984 !!---------------------------------------------------------------------- 985 ! 986 CALL wrk_alloc( jpi,jpj,jpk, zhpi, zu, zv, fsp, xsp, asp, bsp, csp, dsp ) 987 CALL wrk_alloc( jpi,jpj,jpk, zdept, zrhh ) 988 CALL wrk_alloc( jpi,jpj, zsshu_n, zsshv_n ) 989 IF(ln_wd) CALL wrk_alloc( jpi,jpj, zcpx, zcpy ) 964 990 ! 965 991 IF( kt == nit000 ) THEN … … 969 995 ENDIF 970 996 971 !!----------------------------------------------------------------------972 997 ! Local constant initialization 973 998 zcoef0 = - grav 974 znad = 0.0_wp 975 IF( lk_vvl ) znad = 1._wp 999 znad = 1._wp 1000 IF( ln_linssh ) znad = 0._wp 1001 1002 IF(ln_wd) THEN 1003 DO jj = 2, jpjm1 1004 DO ji = 2, jpim1 1005 ll_tmp1 = MIN(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj)) & 1006 & .and. MAX(sshn(ji,jj) + bathy(ji,jj), sshn(ji+1,jj) + bathy(ji+1,jj)) & 1007 & > rn_wdmin1 + rn_wdmin2 1008 ll_tmp2 = MAX(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj)) +& 1009 & rn_wdmin1 + rn_wdmin2 1010 1011 IF(ll_tmp1) THEN 1012 zcpx(ji,jj) = 1.0_wp 1013 ELSE IF(ll_tmp2) THEN 1014 ! no worries about sshn(ji+1,jj)-sshn(ji,jj) = 0, it won't happen ! here 1015 zcpx(ji,jj) = ABS((sshn(ji+1,jj) + bathy(ji+1,jj) - sshn(ji,jj) - bathy(ji,jj)) /& 1016 & (sshn(ji+1,jj) - sshn(ji,jj))) 1017 ELSE 1018 zcpx(ji,jj) = 0._wp 1019 END IF 1020 1021 ll_tmp1 = MIN(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1)) & 1022 & .and. MAX(sshn(ji,jj) + bathy(ji,jj), sshn(ji,jj+1) + bathy(ji,jj+1)) & 1023 & > rn_wdmin1 + rn_wdmin2 1024 ll_tmp2 = MAX(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1)) +& 1025 & rn_wdmin1 + rn_wdmin2 1026 1027 IF(ll_tmp1.OR.ll_tmp2) THEN 1028 zcpy(ji,jj) = 1.0_wp 1029 ELSE IF(ll_tmp2) THEN 1030 ! no worries about sshn(ji,jj+1)-sshn(ji,jj) = 0, it won't happen ! here 1031 zcpy(ji,jj) = ABS((sshn(ji,jj+1) + bathy(ji,jj+1) - sshn(ji,jj) - bathy(ji,jj)) /& 1032 & (sshn(ji,jj+1) - sshn(ji,jj))) 1033 ELSE 1034 zcpy(ji,jj) = 0._wp 1035 END IF 1036 END DO 1037 END DO 1038 CALL lbc_lnk( zcpx, 'U', 1._wp ) ; CALL lbc_lnk( zcpy, 'V', 1._wp ) 1039 ENDIF 976 1040 977 1041 ! Clean 3-D work arrays … … 983 1047 DO ji = 1, jpi 984 1048 jk = mbathy(ji,jj) 985 IF( jk <= 0 ) THEN; zrhh(ji,jj,:) = 0._wp986 ELSE IF(jk == 1) THEN; zrhh(ji,jj,jk+1:jpk) = rhd(ji,jj,jk)987 ELSE IF(jk < jpkm1) THEN1049 IF( jk <= 0 ) THEN ; zrhh(ji,jj, : ) = 0._wp 1050 ELSEIF( jk == 1 ) THEN ; zrhh(ji,jj,jk+1:jpk) = rhd(ji,jj,jk) 1051 ELSEIF( jk < jpkm1 ) THEN 988 1052 DO jkk = jk+1, jpk 989 zrhh(ji,jj,jkk) = interp1( fsde3w(ji,jj,jkk), fsde3w(ji,jj,jkk-1),&990 fsde3w(ji,jj,jkk-2), rhd(ji,jj,jkk-1), rhd(ji,jj,jkk-2))1053 zrhh(ji,jj,jkk) = interp1(gde3w_n(ji,jj,jkk ), gde3w_n(ji,jj,jkk-1), & 1054 & gde3w_n(ji,jj,jkk-2), rhd (ji,jj,jkk-1), rhd(ji,jj,jkk-2)) 991 1055 END DO 992 1056 ENDIF … … 997 1061 DO jj = 1, jpj 998 1062 DO ji = 1, jpi 999 zdept(ji,jj,1) = 0.5_wp * fse3w(ji,jj,1) - sshn(ji,jj) * znad1063 zdept(ji,jj,1) = 0.5_wp * e3w_n(ji,jj,1) - sshn(ji,jj) * znad 1000 1064 END DO 1001 1065 END DO … … 1004 1068 DO jj = 1, jpj 1005 1069 DO ji = 1, jpi 1006 zdept(ji,jj,jk) = zdept(ji,jj,jk-1) + fse3w(ji,jj,jk)1070 zdept(ji,jj,jk) = zdept(ji,jj,jk-1) + e3w_n(ji,jj,jk) 1007 1071 END DO 1008 1072 END DO … … 1015 1079 ! constrained cubic spline interpolation 1016 1080 ! rho(z) = asp + bsp*z + csp*z^2 + dsp*z^3 1017 CALL cspline( fsp,xsp,asp,bsp,csp,dsp,polynomial_type)1081 CALL cspline( fsp, xsp, asp, bsp, csp, dsp, polynomial_type ) 1018 1082 1019 1083 ! Integrate the hydrostatic pressure "zhpi(:,:,:)" at "T(ji,jj,1)" 1020 1084 DO jj = 2, jpj 1021 1085 DO ji = 2, jpi 1022 zrhdt1 = zrhh(ji,jj,1) - interp3(zdept(ji,jj,1),asp(ji,jj,1), & 1023 bsp(ji,jj,1), csp(ji,jj,1), & 1024 dsp(ji,jj,1) ) * 0.25_wp * fse3w(ji,jj,1) 1086 zrhdt1 = zrhh(ji,jj,1) - interp3( zdept(ji,jj,1), asp(ji,jj,1), bsp(ji,jj,1), & 1087 & csp(ji,jj,1), dsp(ji,jj,1) ) * 0.25_wp * e3w_n(ji,jj,1) 1025 1088 1026 1089 ! assuming linear profile across the top half surface layer 1027 zhpi(ji,jj,1) = 0.5_wp * fse3w(ji,jj,1) * zrhdt11090 zhpi(ji,jj,1) = 0.5_wp * e3w_n(ji,jj,1) * zrhdt1 1028 1091 END DO 1029 1092 END DO … … 1033 1096 DO jj = 2, jpj 1034 1097 DO ji = 2, jpi 1035 zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) + &1036 integ_spline(zdept(ji,jj,jk-1), zdept(ji,jj,jk),&1037 asp(ji,jj,jk-1), bsp(ji,jj,jk-1), &1038 csp(ji,jj,jk-1), dsp(ji,jj,jk-1))1098 zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) + & 1099 & integ_spline( zdept(ji,jj,jk-1), zdept(ji,jj,jk), & 1100 & asp (ji,jj,jk-1), bsp (ji,jj,jk-1), & 1101 & csp (ji,jj,jk-1), dsp (ji,jj,jk-1) ) 1039 1102 END DO 1040 1103 END DO … … 1046 1109 DO jj = 2, jpjm1 1047 1110 DO ji = 2, jpim1 1048 zsshu_n(ji,jj) = (e12u(ji,jj) * sshn(ji,jj) + e12u(ji+1, jj) * sshn(ji+1,jj)) * & 1049 & r1_e12u(ji,jj) * umask(ji,jj,1) * 0.5_wp 1050 zsshv_n(ji,jj) = (e12v(ji,jj) * sshn(ji,jj) + e12v(ji+1, jj) * sshn(ji,jj+1)) * & 1051 & r1_e12v(ji,jj) * vmask(ji,jj,1) * 0.5_wp 1111 !!gm BUG ? if it is ssh at u- & v-point then it should be: 1112 ! zsshu_n(ji,jj) = (e1e2t(ji,jj) * sshn(ji,jj) + e1e2t(ji+1,jj) * sshn(ji+1,jj)) * & 1113 ! & r1_e1e2u(ji,jj) * umask(ji,jj,1) * 0.5_wp 1114 ! zsshv_n(ji,jj) = (e1e2t(ji,jj) * sshn(ji,jj) + e1e2t(ji,jj+1) * sshn(ji,jj+1)) * & 1115 ! & r1_e1e2v(ji,jj) * vmask(ji,jj,1) * 0.5_wp 1116 !!gm not this: 1117 zsshu_n(ji,jj) = (e1e2u(ji,jj) * sshn(ji,jj) + e1e2u(ji+1, jj) * sshn(ji+1,jj)) * & 1118 & r1_e1e2u(ji,jj) * umask(ji,jj,1) * 0.5_wp 1119 zsshv_n(ji,jj) = (e1e2v(ji,jj) * sshn(ji,jj) + e1e2v(ji+1, jj) * sshn(ji,jj+1)) * & 1120 & r1_e1e2v(ji,jj) * vmask(ji,jj,1) * 0.5_wp 1052 1121 END DO 1053 1122 END DO 1123 1124 CALL lbc_lnk (zsshu_n, 'U', 1.) 1125 CALL lbc_lnk (zsshv_n, 'V', 1.) 1054 1126 1055 1127 DO jj = 2, jpjm1 1056 1128 DO ji = 2, jpim1 1057 zu(ji,jj,1) = - ( fse3u(ji,jj,1) - zsshu_n(ji,jj) * znad)1058 zv(ji,jj,1) = - ( fse3v(ji,jj,1) - zsshv_n(ji,jj) * znad)1129 zu(ji,jj,1) = - ( e3u_n(ji,jj,1) - zsshu_n(ji,jj) * znad) 1130 zv(ji,jj,1) = - ( e3v_n(ji,jj,1) - zsshv_n(ji,jj) * znad) 1059 1131 END DO 1060 1132 END DO … … 1063 1135 DO jj = 2, jpjm1 1064 1136 DO ji = 2, jpim1 1065 zu(ji,jj,jk) = zu(ji,jj,jk-1) - fse3u(ji,jj,jk)1066 zv(ji,jj,jk) = zv(ji,jj,jk-1) - fse3v(ji,jj,jk)1137 zu(ji,jj,jk) = zu(ji,jj,jk-1) - e3u_n(ji,jj,jk) 1138 zv(ji,jj,jk) = zv(ji,jj,jk-1) - e3v_n(ji,jj,jk) 1067 1139 END DO 1068 1140 END DO … … 1072 1144 DO jj = 2, jpjm1 1073 1145 DO ji = 2, jpim1 1074 zu(ji,jj,jk) = zu(ji,jj,jk) + 0.5_wp * fse3u(ji,jj,jk)1075 zv(ji,jj,jk) = zv(ji,jj,jk) + 0.5_wp * fse3v(ji,jj,jk)1146 zu(ji,jj,jk) = zu(ji,jj,jk) + 0.5_wp * e3u_n(ji,jj,jk) 1147 zv(ji,jj,jk) = zv(ji,jj,jk) + 0.5_wp * e3v_n(ji,jj,jk) 1076 1148 END DO 1077 1149 END DO … … 1081 1153 DO jj = 2, jpjm1 1082 1154 DO ji = 2, jpim1 1083 zu(ji,jj,jk) = min(zu(ji,jj,jk), max(-zdept(ji,jj,jk), -zdept(ji+1,jj,jk)))1084 zu(ji,jj,jk) = max(zu(ji,jj,jk), min(-zdept(ji,jj,jk), -zdept(ji+1,jj,jk)))1085 zv(ji,jj,jk) = min(zv(ji,jj,jk), max(-zdept(ji,jj,jk), -zdept(ji,jj+1,jk)))1086 zv(ji,jj,jk) = max(zv(ji,jj,jk), min(-zdept(ji,jj,jk), -zdept(ji,jj+1,jk)))1155 zu(ji,jj,jk) = MIN( zu(ji,jj,jk) , MAX( -zdept(ji,jj,jk) , -zdept(ji+1,jj,jk) ) ) 1156 zu(ji,jj,jk) = MAX( zu(ji,jj,jk) , MIN( -zdept(ji,jj,jk) , -zdept(ji+1,jj,jk) ) ) 1157 zv(ji,jj,jk) = MIN( zv(ji,jj,jk) , MAX( -zdept(ji,jj,jk) , -zdept(ji,jj+1,jk) ) ) 1158 zv(ji,jj,jk) = MAX( zv(ji,jj,jk) , MIN( -zdept(ji,jj,jk) , -zdept(ji,jj+1,jk) ) ) 1087 1159 END DO 1088 1160 END DO … … 1142 1214 ! update the momentum trends in u direction 1143 1215 1144 zdpdx1 = zcoef0 / e1u(ji,jj) * (zhpi(ji+1,jj,jk) - zhpi(ji,jj,jk))1145 IF( lk_vvl) THEN1146 zdpdx2 = zcoef0 /e1u(ji,jj) * &1147 1216 zdpdx1 = zcoef0 * r1_e1u(ji,jj) * ( zhpi(ji+1,jj,jk) - zhpi(ji,jj,jk) ) 1217 IF( .NOT.ln_linssh ) THEN 1218 zdpdx2 = zcoef0 * r1_e1u(ji,jj) * & 1219 & ( REAL(jis-jid, wp) * (zpwes + zpwed) + (sshn(ji+1,jj)-sshn(ji,jj)) ) 1148 1220 ELSE 1149 zdpdx2 = zcoef0 /e1u(ji,jj) * REAL(jis-jid, wp) * (zpwes + zpwed)1221 zdpdx2 = zcoef0 * r1_e1u(ji,jj) * REAL(jis-jid, wp) * (zpwes + zpwed) 1150 1222 ENDIF 1151 1152 ua(ji,jj,jk) = ua(ji,jj,jk) + (zdpdx1 + zdpdx2) * & 1153 & umask(ji,jj,jk) * tmask(ji,jj,jk) * tmask(ji+1,jj,jk) 1223 IF(ln_wd) THEN 1224 zdpdx1 = zdpdx1 * zcpx(ji,jj) 1225 zdpdx2 = zdpdx2 * zcpx(ji,jj) 1226 ENDIF 1227 ua(ji,jj,jk) = ua(ji,jj,jk) + (zdpdx1 + zdpdx2) * umask(ji,jj,jk) 1154 1228 ENDIF 1155 1229 … … 1199 1273 ! update the momentum trends in v direction 1200 1274 1201 zdpdy1 = zcoef0 / e2v(ji,jj) * (zhpi(ji,jj+1,jk) - zhpi(ji,jj,jk))1202 IF( lk_vvl) THEN1203 zdpdy2 = zcoef0 /e2v(ji,jj) * &1204 1275 zdpdy1 = zcoef0 * r1_e2v(ji,jj) * ( zhpi(ji,jj+1,jk) - zhpi(ji,jj,jk) ) 1276 IF( .NOT.ln_linssh ) THEN 1277 zdpdy2 = zcoef0 * r1_e2v(ji,jj) * & 1278 ( REAL(jjs-jjd, wp) * (zpnss + zpnsd) + (sshn(ji,jj+1)-sshn(ji,jj)) ) 1205 1279 ELSE 1206 zdpdy2 = zcoef0 /e2v(ji,jj) * REAL(jjs-jjd, wp) * (zpnss + zpnsd )1280 zdpdy2 = zcoef0 * r1_e2v(ji,jj) * REAL(jjs-jjd, wp) * (zpnss + zpnsd ) 1207 1281 ENDIF 1208 1209 va(ji,jj,jk) = va(ji,jj,jk) + (zdpdy1 + zdpdy2)*& 1210 & vmask(ji,jj,jk)*tmask(ji,jj,jk)*tmask(ji,jj+1,jk) 1282 IF(ln_wd) THEN 1283 zdpdy1 = zdpdy1 * zcpy(ji,jj) 1284 zdpdy2 = zdpdy2 * zcpy(ji,jj) 1285 ENDIF 1286 1287 va(ji,jj,jk) = va(ji,jj,jk) + (zdpdy1 + zdpdy2) * vmask(ji,jj,jk) 1211 1288 ENDIF 1212 1213 1214 1215 1216 END DO1217 !1218 CALL wrk_dealloc( jpi,jpj,jpk, zhpi, zu, zv, fsp, xsp, asp, bsp, csp, dsp)1219 CALL wrk_dealloc( jpi,jpj, jpk, zdept, zrhh)1220 CALL wrk_dealloc( jpi,jpj, zsshu_n, zsshv_n)1289 ! 1290 END DO 1291 END DO 1292 END DO 1293 ! 1294 CALL wrk_dealloc( jpi,jpj,jpk, zhpi, zu, zv, fsp, xsp, asp, bsp, csp, dsp ) 1295 CALL wrk_dealloc( jpi,jpj,jpk, zdept, zrhh ) 1296 CALL wrk_dealloc( jpi,jpj, zsshu_n, zsshv_n ) 1297 IF(ln_wd) CALL wrk_dealloc( jpi,jpj, zcpx, zcpy ) 1221 1298 ! 1222 1299 END SUBROUTINE hpg_prj 1223 1300 1224 1301 1225 SUBROUTINE cspline( fsp, xsp, asp, bsp, csp, dsp, polynomial_type)1302 SUBROUTINE cspline( fsp, xsp, asp, bsp, csp, dsp, polynomial_type ) 1226 1303 !!---------------------------------------------------------------------- 1227 1304 !! *** ROUTINE cspline *** … … 1233 1310 !! Reference: CJC Kruger, Constrained Cubic Spline Interpoltation 1234 1311 !!---------------------------------------------------------------------- 1235 IMPLICIT NONE 1236 REAL(wp), DIMENSION(:,:,:), INTENT(in) :: fsp, xsp ! value and coordinate 1237 REAL(wp), DIMENSION(:,:,:), INTENT(out) :: asp, bsp, csp, dsp ! coefficients of 1238 ! the interpoated function 1239 INTEGER, INTENT(in) :: polynomial_type ! 1: cubic spline 1240 ! 2: Linear 1312 REAL(wp), DIMENSION(:,:,:), INTENT(in ) :: fsp, xsp ! value and coordinate 1313 REAL(wp), DIMENSION(:,:,:), INTENT( out) :: asp, bsp, csp, dsp ! coefficients of the interpoated function 1314 INTEGER , INTENT(in ) :: polynomial_type ! 1: cubic spline ; 2: Linear 1241 1315 ! 1242 1316 INTEGER :: ji, jj, jk ! dummy loop indices … … 1246 1320 REAL(wp) :: zdf(size(fsp,3)) 1247 1321 !!---------------------------------------------------------------------- 1248 1322 ! 1323 !!gm WHAT !!!!! THIS IS VERY DANGEROUS !!!!! 1249 1324 jpi = size(fsp,1) 1250 1325 jpj = size(fsp,2) 1251 1326 jpkm1 = size(fsp,3) - 1 1252 1253 1327 ! 1254 1328 IF (polynomial_type == 1) THEN ! Constrained Cubic Spline 1255 1329 DO ji = 1, jpi … … 1284 1358 1285 1359 zdf(1) = 1.5_wp * ( fsp(ji,jj,2) - fsp(ji,jj,1) ) / & 1286 & ( xsp(ji,jj,2) - xsp(ji,jj,1) ) - 0.5_wp * zdf(2)1360 & ( xsp(ji,jj,2) - xsp(ji,jj,1) ) - 0.5_wp * zdf(2) 1287 1361 zdf(jpkm1) = 1.5_wp * ( fsp(ji,jj,jpkm1) - fsp(ji,jj,jpkm1-1) ) / & 1288 & ( xsp(ji,jj,jpkm1) - xsp(ji,jj,jpkm1-1) ) - & 1289 & 0.5_wp * zdf(jpkm1 - 1) 1362 & ( xsp(ji,jj,jpkm1) - xsp(ji,jj,jpkm1-1) ) - 0.5_wp * zdf(jpkm1 - 1) 1290 1363 1291 1364 DO jk = 1, jpkm1 - 1 … … 1310 1383 END DO 1311 1384 1312 ELSE IF (polynomial_type == 2) THEN ! Linear1385 ELSEIF ( polynomial_type == 2 ) THEN ! Linear 1313 1386 DO ji = 1, jpi 1314 1387 DO jj = 1, jpj … … 1341 1414 !! extrapolation is also permitted (no value limit) 1342 1415 !!---------------------------------------------------------------------- 1343 IMPLICIT NONE1344 1416 REAL(wp), INTENT(in) :: x, xl, xr, fl, fr 1345 1417 REAL(wp) :: f ! result of the interpolation (extrapolation) 1346 1418 REAL(wp) :: zdeltx 1347 1419 !!---------------------------------------------------------------------- 1348 1420 ! 1349 1421 zdeltx = xr - xl 1350 IF( abs(zdeltx) <= 10._wp * EPSILON(x)) THEN1351 f = 0.5_wp * (fl + fr)1422 IF( abs(zdeltx) <= 10._wp * EPSILON(x) ) THEN 1423 f = 0.5_wp * (fl + fr) 1352 1424 ELSE 1353 f = ( (x - xl ) * fr - ( x - xr ) * fl ) / zdeltx1425 f = ( (x - xl ) * fr - ( x - xr ) * fl ) / zdeltx 1354 1426 ENDIF 1355 1427 ! 1356 1428 END FUNCTION interp1 1357 1429 1358 1430 1359 FUNCTION interp2( x, a, b, c, d) RESULT(f)1431 FUNCTION interp2( x, a, b, c, d ) RESULT(f) 1360 1432 !!---------------------------------------------------------------------- 1361 1433 !! *** ROUTINE interp1 *** … … 1366 1438 !! 1367 1439 !!---------------------------------------------------------------------- 1368 IMPLICIT NONE1369 1440 REAL(wp), INTENT(in) :: x, a, b, c, d 1370 1441 REAL(wp) :: f ! value from the interpolation 1371 1442 !!---------------------------------------------------------------------- 1372 1443 ! 1373 1444 f = a + x* ( b + x * ( c + d * x ) ) 1374 1445 ! 1375 1446 END FUNCTION interp2 1376 1447 1377 1448 1378 FUNCTION interp3( x, a, b, c, d) RESULT(f)1449 FUNCTION interp3( x, a, b, c, d ) RESULT(f) 1379 1450 !!---------------------------------------------------------------------- 1380 1451 !! *** ROUTINE interp1 *** … … 1386 1457 !! 1387 1458 !!---------------------------------------------------------------------- 1388 IMPLICIT NONE1389 1459 REAL(wp), INTENT(in) :: x, a, b, c, d 1390 1460 REAL(wp) :: f ! value from the interpolation 1391 1461 !!---------------------------------------------------------------------- 1392 1462 ! 1393 1463 f = b + x * ( 2._wp * c + 3._wp * d * x) 1394 1464 ! 1395 1465 END FUNCTION interp3 1396 1466 1397 1467 1398 FUNCTION integ_spline( xl, xr, a, b, c, d) RESULT(f)1468 FUNCTION integ_spline( xl, xr, a, b, c, d ) RESULT(f) 1399 1469 !!---------------------------------------------------------------------- 1400 1470 !! *** ROUTINE interp1 *** … … 1405 1475 !! 1406 1476 !!---------------------------------------------------------------------- 1407 IMPLICIT NONE1408 1477 REAL(wp), INTENT(in) :: xl, xr, a, b, c, d 1409 1478 REAL(wp) :: za1, za2, za3 1410 1479 REAL(wp) :: f ! integration result 1411 1480 !!---------------------------------------------------------------------- 1412 1481 ! 1413 1482 za1 = 0.5_wp * b 1414 1483 za2 = c / 3.0_wp 1415 1484 za3 = 0.25_wp * d 1416 1485 ! 1417 1486 f = xr * ( a + xr * ( za1 + xr * ( za2 + za3 * xr ) ) ) - & 1418 1487 & xl * ( a + xl * ( za1 + xl * ( za2 + za3 * xl ) ) ) 1419 1488 ! 1420 1489 END FUNCTION integ_spline 1421 1490
Note: See TracChangeset
for help on using the changeset viewer.
