- Timestamp:
- 2015-02-20T20:11:47+01:00 (9 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2015/dev_r5094_UKMO_ISFCLEAN/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90
r4990 r5098 16 16 !! 3.4 ! 2011-11 (H. Liu) hpg_prj: Original code for s-coordinates 17 17 !! ! (A. Coward) suppression of hel, wdj and rot options 18 !! 3.6 ! 2014-11 (P. Mathiot) hpg_isf: original code for ice shelf cavity 18 19 !!---------------------------------------------------------------------- 19 20 … … 25 26 !! hpg_zps : z-coordinate plus partial steps (interpolation) 26 27 !! hpg_sco : s-coordinate (standard jacobian formulation) 28 !! hpg_isf : s-coordinate (sco formulation) adapted to ice shelf 27 29 !! hpg_djc : s-coordinate (Density Jacobian with Cubic polynomial) 28 30 !! hpg_prj : s-coordinate (Pressure Jacobian with Cubic polynomial) … … 55 57 LOGICAL , PUBLIC :: ln_hpg_djc !: s-coordinate (Density Jacobian with Cubic polynomial) 56 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 57 60 LOGICAL , PUBLIC :: ln_dynhpg_imp !: semi-implicite hpg flag 58 61 … … 97 100 CASE ( 3 ) ; CALL hpg_djc ( kt ) ! s-coordinate (Density Jacobian with Cubic polynomial) 98 101 CASE ( 4 ) ; CALL hpg_prj ( kt ) ! s-coordinate (Pressure Jacobian scheme) 102 CASE ( 5 ) ; CALL hpg_isf ( kt ) ! s-coordinate similar to sco modify for ice shelf 99 103 END SELECT 100 104 ! … … 128 132 !! 129 133 NAMELIST/namdyn_hpg/ ln_hpg_zco, ln_hpg_zps, ln_hpg_sco, & 130 & ln_hpg_djc, ln_hpg_prj, ln_ dynhpg_imp134 & ln_hpg_djc, ln_hpg_prj, ln_hpg_isf, ln_dynhpg_imp 131 135 !!---------------------------------------------------------------------- 132 136 ! … … 148 152 WRITE(numout,*) ' z-coord. - partial steps (interpolation) ln_hpg_zps = ', ln_hpg_zps 149 153 WRITE(numout,*) ' s-coord. (standard jacobian formulation) ln_hpg_sco = ', ln_hpg_sco 154 WRITE(numout,*) ' s-coord. (standard jacobian formulation) for isf ln_hpg_isf = ', ln_hpg_isf 150 155 WRITE(numout,*) ' s-coord. (Density Jacobian: Cubic polynomial) ln_hpg_djc = ', ln_hpg_djc 151 156 WRITE(numout,*) ' s-coord. (Pressure Jacobian: Cubic polynomial) ln_hpg_prj = ', ln_hpg_prj … … 158 163 & either ln_hpg_sco or ln_hpg_prj instead') 159 164 ! 160 IF( lk_vvl .AND. .NOT. (ln_hpg_sco.OR.ln_hpg_prj ) ) &165 IF( lk_vvl .AND. .NOT. (ln_hpg_sco.OR.ln_hpg_prj.OR.ln_hpg_isf) ) & 161 166 & CALL ctl_stop('dyn_hpg_init : variable volume key_vvl requires:& 162 167 & the standard jacobian formulation hpg_sco or & … … 169 174 IF( ln_hpg_djc ) nhpg = 3 170 175 IF( ln_hpg_prj ) nhpg = 4 176 IF( ln_hpg_isf ) nhpg = 5 171 177 ! 172 178 ! ! Consistency check … … 177 183 IF( ln_hpg_djc ) ioptio = ioptio + 1 178 184 IF( ln_hpg_prj ) ioptio = ioptio + 1 185 IF( ln_hpg_isf ) ioptio = ioptio + 1 179 186 IF( ioptio /= 1 ) CALL ctl_stop( 'NO or several hydrostatic pressure gradient options used' ) 180 IF( ( ln_hpg_zco .OR. ln_hpg_zps .OR. ln_hpg_djc .OR. ln_hpg_prj ) .AND. nn_isf .NE. 0) &181 & CALL ctl_stop( 'Only hpg_ scohas been corrected to work with ice shelf cavity.' )187 IF( ( .NOT. ln_hpg_isf ) .AND. ln_isfcav ) & 188 & CALL ctl_stop( 'Only hpg_isf has been corrected to work with ice shelf cavity.' ) 182 189 ! 183 190 END SUBROUTINE dyn_hpg_init … … 345 352 END SUBROUTINE hpg_zps 346 353 347 348 354 SUBROUTINE hpg_sco( kt ) 349 355 !!--------------------------------------------------------------------- … … 366 372 INTEGER, INTENT(in) :: kt ! ocean time-step index 367 373 !! 374 INTEGER :: ji, jj, jk ! dummy loop indices 375 REAL(wp) :: zcoef0, zuap, zvap, znad ! temporary scalars 376 REAL(wp), POINTER, DIMENSION(:,:,:) :: zhpi, zhpj 377 !!---------------------------------------------------------------------- 378 ! 379 CALL wrk_alloc( jpi,jpj,jpk, zhpi, zhpj ) 380 ! 381 IF( kt == nit000 ) THEN 382 IF(lwp) WRITE(numout,*) 383 IF(lwp) WRITE(numout,*) 'dyn:hpg_sco : hydrostatic pressure gradient trend' 384 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ s-coordinate case, OPA original scheme used' 385 ENDIF 386 387 ! Local constant initialization 388 zcoef0 = - grav * 0.5_wp 389 ! To use density and not density anomaly 390 IF ( lk_vvl ) THEN ; znad = 1._wp ! Variable volume 391 ELSE ; znad = 0._wp ! Fixed volume 392 ENDIF 393 394 ! Surface value 395 DO jj = 2, jpjm1 396 DO ji = fs_2, fs_jpim1 ! vector opt. 397 ! hydrostatic pressure gradient along s-surfaces 398 zhpi(ji,jj,1) = zcoef0 / e1u(ji,jj) * ( fse3w(ji+1,jj ,1) * ( znad + rhd(ji+1,jj ,1) ) & 399 & - fse3w(ji ,jj ,1) * ( znad + rhd(ji ,jj ,1) ) ) 400 zhpj(ji,jj,1) = zcoef0 / e2v(ji,jj) * ( fse3w(ji ,jj+1,1) * ( znad + rhd(ji ,jj+1,1) ) & 401 & - fse3w(ji ,jj ,1) * ( znad + rhd(ji ,jj ,1) ) ) 402 ! s-coordinate pressure gradient correction 403 zuap = -zcoef0 * ( rhd (ji+1,jj,1) + rhd (ji,jj,1) + 2._wp * znad ) & 404 & * ( fsde3w(ji+1,jj,1) - fsde3w(ji,jj,1) ) / e1u(ji,jj) 405 zvap = -zcoef0 * ( rhd (ji,jj+1,1) + rhd (ji,jj,1) + 2._wp * znad ) & 406 & * ( fsde3w(ji,jj+1,1) - fsde3w(ji,jj,1) ) / e2v(ji,jj) 407 ! add to the general momentum trend 408 ua(ji,jj,1) = ua(ji,jj,1) + zhpi(ji,jj,1) + zuap 409 va(ji,jj,1) = va(ji,jj,1) + zhpj(ji,jj,1) + zvap 410 END DO 411 END DO 412 413 ! interior value (2=<jk=<jpkm1) 414 DO jk = 2, jpkm1 415 DO jj = 2, jpjm1 416 DO ji = fs_2, fs_jpim1 ! vector opt. 417 ! hydrostatic pressure gradient along s-surfaces 418 zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) + zcoef0 / e1u(ji,jj) & 419 & * ( fse3w(ji+1,jj,jk) * ( rhd(ji+1,jj,jk) + rhd(ji+1,jj,jk-1) + 2*znad ) & 420 & - fse3w(ji ,jj,jk) * ( rhd(ji ,jj,jk) + rhd(ji ,jj,jk-1) + 2*znad ) ) 421 zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1) + zcoef0 / e2v(ji,jj) & 422 & * ( fse3w(ji,jj+1,jk) * ( rhd(ji,jj+1,jk) + rhd(ji,jj+1,jk-1) + 2*znad ) & 423 & - fse3w(ji,jj ,jk) * ( rhd(ji,jj, jk) + rhd(ji,jj ,jk-1) + 2*znad ) ) 424 ! s-coordinate pressure gradient correction 425 zuap = -zcoef0 * ( rhd (ji+1,jj ,jk) + rhd (ji,jj,jk) + 2._wp * znad ) & 426 & * ( fsde3w(ji+1,jj ,jk) - fsde3w(ji,jj,jk) ) / e1u(ji,jj) 427 zvap = -zcoef0 * ( rhd (ji ,jj+1,jk) + rhd (ji,jj,jk) + 2._wp * znad ) & 428 & * ( fsde3w(ji ,jj+1,jk) - fsde3w(ji,jj,jk) ) / e2v(ji,jj) 429 ! add to the general momentum trend 430 ua(ji,jj,jk) = ua(ji,jj,jk) + zhpi(ji,jj,jk) + zuap 431 va(ji,jj,jk) = va(ji,jj,jk) + zhpj(ji,jj,jk) + zvap 432 END DO 433 END DO 434 END DO 435 ! 436 CALL wrk_dealloc( jpi,jpj,jpk, zhpi, zhpj ) 437 ! 438 END SUBROUTINE hpg_sco 439 440 SUBROUTINE hpg_isf( kt ) 441 !!--------------------------------------------------------------------- 442 !! *** ROUTINE hpg_sco *** 443 !! 444 !! ** Method : s-coordinate case. Jacobian scheme. 445 !! The now hydrostatic pressure gradient at a given level, jk, 446 !! is computed by taking the vertical integral of the in-situ 447 !! density gradient along the model level from the suface to that 448 !! level. s-coordinates (ln_sco): a corrective term is added 449 !! to the horizontal pressure gradient : 450 !! zhpi = grav ..... + 1/e1u mi(rhd) di[ grav dep3w ] 451 !! zhpj = grav ..... + 1/e2v mj(rhd) dj[ grav dep3w ] 452 !! add it to the general momentum trend (ua,va). 453 !! ua = ua - 1/e1u * zhpi 454 !! va = va - 1/e2v * zhpj 455 !! iceload is added and partial cell case are added to the top and bottom 456 !! 457 !! ** Action : - Update (ua,va) with the now hydrastatic pressure trend 458 !!---------------------------------------------------------------------- 459 INTEGER, INTENT(in) :: kt ! ocean time-step index 460 !! 368 461 INTEGER :: ji, jj, jk, iku, ikv, ikt, iktp1i, iktp1j ! dummy loop indices 369 462 REAL(wp) :: zcoef0, zuap, zvap, znad, ze3wu, ze3wv, zuapint, zvapint, zhpjint, zhpiint, zdzwt, zdzwtjp1, zdzwtip1 ! temporary scalars … … 379 472 IF( kt == nit000 ) THEN 380 473 IF(lwp) WRITE(numout,*) 381 IF(lwp) WRITE(numout,*) 'dyn:hpg_ sco : hydrostatic pressure gradient trend'474 IF(lwp) WRITE(numout,*) 'dyn:hpg_isf : hydrostatic pressure gradient trend for ice shelf' 382 475 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ s-coordinate case, OPA original scheme used' 383 476 ENDIF … … 565 658 !================================================================================== 566 659 567 # if defined key_vectopt_loop568 jj = 1569 DO ji = jpi+2, jpij-jpi-1 ! vector opt. (forced unrolling)570 # else571 660 DO jj = 2, jpjm1 572 661 DO ji = 2, jpim1 573 # endif574 662 iku = mbku(ji,jj) 575 663 ikv = mbkv(ji,jj) … … 598 686 va(ji,jj,ikv) = va(ji,jj,ikv) + zhpj(ji,jj,ikv) + zvap 599 687 END IF 600 # if ! defined key_vectopt_loop 601 END DO 602 # endif 688 END DO 603 689 END DO 604 690 … … 610 696 CALL wrk_dealloc( jpi,jpj, ze3w, zp, zrhdtop_isf, zrhdtop_oce, ziceload, zdept, zpshpi, zpshpj) 611 697 ! 612 END SUBROUTINE hpg_ sco698 END SUBROUTINE hpg_isf 613 699 614 700
Note: See TracChangeset
for help on using the changeset viewer.
