Changeset 3108 for branches/2011/dev_NOC_UKMO_MERGE/NEMOGCM/NEMO
- Timestamp:
- 2011-11-15T12:55:23+01:00 (13 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2011/dev_NOC_UKMO_MERGE/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90
r3102 r3108 14 14 !! - ! 2005-11 (G. Madec) style & small optimisation 15 15 !! 3.3 ! 2010-10 (C. Ethe, G. Madec) reorganisation of initialisation phase 16 !! 3.4 ! 2011-11 (A. Coward, H. Liu) introduction of prj scheme; 17 !! ! suppression of hel, wdj and rot options 16 18 !!---------------------------------------------------------------------- 17 19 … … 23 25 !! hpg_zps : z-coordinate plus partial steps (interpolation) 24 26 !! hpg_sco : s-coordinate (standard jacobian formulation) 25 !! hpg_hel : s-coordinate (helsinki modification)26 !! hpg_wdj : s-coordinate (weighted density jacobian)27 27 !! hpg_djc : s-coordinate (Density Jacobian with Cubic polynomial) 28 !! hpg_rot : s-coordinate (ROTated axes scheme)29 28 !! hpg_prj : s-coordinate (Pressure Jacobian with Cubic polynomial) 30 29 !!---------------------------------------------------------------------- … … 49 48 LOGICAL , PUBLIC :: ln_hpg_zps = .FALSE. !: z-coordinate - partial steps (interpolation) 50 49 LOGICAL , PUBLIC :: ln_hpg_sco = .FALSE. !: s-coordinate (standard jacobian formulation) 51 LOGICAL , PUBLIC :: ln_hpg_hel = .FALSE. !: s-coordinate (helsinki modification)52 LOGICAL , PUBLIC :: ln_hpg_wdj = .FALSE. !: s-coordinate (weighted density jacobian)53 50 LOGICAL , PUBLIC :: ln_hpg_djc = .FALSE. !: s-coordinate (Density Jacobian with Cubic polynomial) 54 LOGICAL , PUBLIC :: ln_hpg_rot = .FALSE. !: s-coordinate (ROTated axes scheme)55 51 LOGICAL , PUBLIC :: ln_hpg_prj = .FALSE. !: s-coordinate (Pressure Jacobian scheme) 56 REAL(wp), PUBLIC :: rn_gamma = 0._wp !: weighting coefficient57 52 LOGICAL , PUBLIC :: ln_dynhpg_imp = .FALSE. !: semi-implicite hpg flag 58 53 … … 98 93 CASE ( 1 ) ; CALL hpg_zps ( kt ) ! z-coordinate plus partial steps (interpolation) 99 94 CASE ( 2 ) ; CALL hpg_sco ( kt ) ! s-coordinate (standard jacobian formulation) 100 CASE ( 3 ) ; CALL hpg_hel ( kt ) ! s-coordinate (helsinki modification) 101 CASE ( 4 ) ; CALL hpg_wdj ( kt ) ! s-coordinate (weighted density jacobian) 102 CASE ( 5 ) ; CALL hpg_djc ( kt ) ! s-coordinate (Density Jacobian with Cubic polynomial) 103 CASE ( 6 ) ; CALL hpg_rot ( kt ) ! s-coordinate (ROTated axes scheme) 104 CASE ( 7 ) ; CALL hpg_prj ( kt ) ! s-coordinate (Pressure Jacobian scheme) 95 CASE ( 3 ) ; CALL hpg_djc ( kt ) ! s-coordinate (Density Jacobian with Cubic polynomial) 96 CASE ( 4 ) ; CALL hpg_prj ( kt ) ! s-coordinate (Pressure Jacobian scheme) 105 97 END SELECT 106 98 ! … … 131 123 INTEGER :: ioptio = 0 ! temporary integer 132 124 !! 133 NAMELIST/namdyn_hpg/ ln_hpg_zco, ln_hpg_zps, ln_hpg_sco, ln_hpg_hel, & 134 & ln_hpg_wdj, ln_hpg_djc, ln_hpg_rot, ln_hpg_prj, & 135 & rn_gamma , ln_dynhpg_imp 125 NAMELIST/namdyn_hpg/ ln_hpg_zco, ln_hpg_zps, ln_hpg_sco, & 126 & ln_hpg_djc, ln_hpg_prj, ln_dynhpg_imp 136 127 !!---------------------------------------------------------------------- 137 128 ! … … 147 138 WRITE(numout,*) ' z-coord. - partial steps (interpolation) ln_hpg_zps = ', ln_hpg_zps 148 139 WRITE(numout,*) ' s-coord. (standard jacobian formulation) ln_hpg_sco = ', ln_hpg_sco 149 WRITE(numout,*) ' s-coord. (helsinki modification) ln_hpg_hel = ', ln_hpg_hel150 WRITE(numout,*) ' s-coord. (weighted density jacobian) ln_hpg_wdj = ', ln_hpg_wdj151 140 WRITE(numout,*) ' s-coord. (Density Jacobian: Cubic polynomial) ln_hpg_djc = ', ln_hpg_djc 152 WRITE(numout,*) ' s-coord. (ROTated axes scheme) ln_hpg_rot = ', ln_hpg_rot153 141 WRITE(numout,*) ' s-coord. (Pressure Jacobian: Cubic polynomial) ln_hpg_prj = ', ln_hpg_prj 154 WRITE(numout,*) ' weighting coeff. (wdj scheme) rn_gamma = ', rn_gamma155 142 WRITE(numout,*) ' time stepping: centered (F) or semi-implicit (T) ln_dynhpg_imp = ', ln_dynhpg_imp 156 143 ENDIF … … 165 152 IF( ln_hpg_zps ) nhpg = 1 166 153 IF( ln_hpg_sco ) nhpg = 2 167 IF( ln_hpg_hel ) nhpg = 3 168 IF( ln_hpg_wdj ) nhpg = 4 169 IF( ln_hpg_djc ) nhpg = 5 170 IF( ln_hpg_rot ) nhpg = 6 171 IF( ln_hpg_prj ) nhpg = 7 154 IF( ln_hpg_djc ) nhpg = 3 155 IF( ln_hpg_prj ) nhpg = 4 172 156 ! 173 157 ! ! Consistency check … … 176 160 IF( ln_hpg_zps ) ioptio = ioptio + 1 177 161 IF( ln_hpg_sco ) ioptio = ioptio + 1 178 IF( ln_hpg_hel ) ioptio = ioptio + 1179 IF( ln_hpg_wdj ) ioptio = ioptio + 1180 162 IF( ln_hpg_djc ) ioptio = ioptio + 1 181 IF( ln_hpg_rot ) ioptio = ioptio + 1182 163 IF( ln_hpg_prj ) ioptio = ioptio + 1 183 164 IF( ioptio /= 1 ) CALL ctl_stop( 'NO or several hydrostatic pressure gradient options used' ) … … 428 409 END SUBROUTINE hpg_sco 429 410 430 431 SUBROUTINE hpg_hel( kt )432 !!---------------------------------------------------------------------433 !! *** ROUTINE hpg_hel ***434 !!435 !! ** Method : s-coordinate case.436 !! The now hydrostatic pressure gradient at a given level437 !! jk is computed by taking the vertical integral of the in-situ438 !! density gradient along the model level from the suface to that439 !! level. s-coordinates (ln_sco): a corrective term is added440 !! to the horizontal pressure gradient :441 !! zhpi = grav ..... + 1/e1u mi(rhd) di[ grav dep3w ]442 !! zhpj = grav ..... + 1/e2v mj(rhd) dj[ grav dep3w ]443 !! add it to the general momentum trend (ua,va).444 !! ua = ua - 1/e1u * zhpi445 !! va = va - 1/e2v * zhpj446 !!447 !! ** Action : - Update (ua,va) with the now hydrastatic pressure trend448 !! - Save the trend (l_trddyn=T)449 !!----------------------------------------------------------------------450 USE oce, ONLY: zhpi => ta , zhpj => sa ! (ta,sa) used as 3D workspace451 !!452 INTEGER, INTENT(in) :: kt ! ocean time-step index453 !!454 INTEGER :: ji, jj, jk ! dummy loop indices455 REAL(wp) :: zcoef0, zuap, zvap ! temporary scalars456 !!----------------------------------------------------------------------457 458 IF( kt == nit000 ) THEN459 IF(lwp) WRITE(numout,*)460 IF(lwp) WRITE(numout,*) 'dyn:hpg_hel : hydrostatic pressure gradient trend'461 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ s-coordinate case, helsinki modified scheme'462 ENDIF463 464 ! Local constant initialization465 zcoef0 = - grav * 0.5_wp466 467 ! Surface value468 DO jj = 2, jpjm1469 DO ji = fs_2, fs_jpim1 ! vector opt.470 ! hydrostatic pressure gradient along s-surfaces471 zhpi(ji,jj,1) = zcoef0 / e1u(ji,jj) * ( fse3t(ji+1,jj ,1) * rhd(ji+1,jj ,1) &472 & - fse3t(ji ,jj ,1) * rhd(ji ,jj ,1) )473 zhpj(ji,jj,1) = zcoef0 / e2v(ji,jj) * ( fse3t(ji ,jj+1,1) * rhd(ji ,jj+1,1) &474 & - fse3t(ji ,jj ,1) * rhd(ji ,jj ,1) )475 ! s-coordinate pressure gradient correction476 zuap = -zcoef0 * ( rhd (ji+1,jj,1) + rhd (ji,jj,1) ) &477 & * ( fsdept(ji+1,jj,1) - fsdept(ji,jj,1) ) / e1u(ji,jj)478 zvap = -zcoef0 * ( rhd (ji,jj+1,1) + rhd (ji,jj,1) ) &479 & * ( fsdept(ji,jj+1,1) - fsdept(ji,jj,1) ) / e2v(ji,jj)480 ! add to the general momentum trend481 ua(ji,jj,1) = ua(ji,jj,1) + zhpi(ji,jj,1) + zuap482 va(ji,jj,1) = va(ji,jj,1) + zhpj(ji,jj,1) + zvap483 END DO484 END DO485 !486 ! interior value (2=<jk=<jpkm1)487 DO jk = 2, jpkm1488 DO jj = 2, jpjm1489 DO ji = fs_2, fs_jpim1 ! vector opt.490 ! hydrostatic pressure gradient along s-surfaces491 zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) &492 & + zcoef0 / e1u(ji,jj) * ( fse3t(ji+1,jj,jk ) * rhd(ji+1,jj,jk) &493 & -fse3t(ji ,jj,jk ) * rhd(ji ,jj,jk) ) &494 & + zcoef0 / e1u(ji,jj) * ( fse3t(ji+1,jj,jk-1) * rhd(ji+1,jj,jk-1) &495 & -fse3t(ji ,jj,jk-1) * rhd(ji ,jj,jk-1) )496 zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1) &497 & + zcoef0 / e2v(ji,jj) * ( fse3t(ji,jj+1,jk ) * rhd(ji,jj+1,jk) &498 & -fse3t(ji,jj ,jk ) * rhd(ji,jj, jk) ) &499 & + zcoef0 / e2v(ji,jj) * ( fse3t(ji,jj+1,jk-1) * rhd(ji,jj+1,jk-1) &500 & -fse3t(ji,jj ,jk-1) * rhd(ji,jj, jk-1) )501 ! s-coordinate pressure gradient correction502 zuap = - zcoef0 * ( rhd (ji+1,jj,jk) + rhd (ji,jj,jk) ) &503 & * ( fsdept(ji+1,jj,jk) - fsdept(ji,jj,jk) ) / e1u(ji,jj)504 zvap = - zcoef0 * ( rhd (ji,jj+1,jk) + rhd (ji,jj,jk) ) &505 & * ( fsdept(ji,jj+1,jk) - fsdept(ji,jj,jk) ) / e2v(ji,jj)506 ! add to the general momentum trend507 ua(ji,jj,jk) = ua(ji,jj,jk) + zhpi(ji,jj,jk) + zuap508 va(ji,jj,jk) = va(ji,jj,jk) + zhpj(ji,jj,jk) + zvap509 END DO510 END DO511 END DO512 !513 END SUBROUTINE hpg_hel514 515 516 SUBROUTINE hpg_wdj( kt )517 !!---------------------------------------------------------------------518 !! *** ROUTINE hpg_wdj ***519 !!520 !! ** Method : Weighted Density Jacobian (wdj) scheme (song 1998)521 !! The weighting coefficients from the namelist parameter rn_gamma522 !! (alpha=0.5-rn_gamma ; beta=1-alpha=0.5+rn_gamma523 !!524 !! Reference : Song, Mon. Wea. Rev., 126, 3213-3230, 1998.525 !!----------------------------------------------------------------------526 USE oce, ONLY: zhpi => ta , zhpj => sa ! (ta,sa) used as 3D workspace527 !!528 INTEGER, INTENT(in) :: kt ! ocean time-step index529 !!530 INTEGER :: ji, jj, jk ! dummy loop indices531 REAL(wp) :: zcoef0, zuap, zvap ! temporary scalars532 REAL(wp) :: zalph , zbeta ! " "533 !!----------------------------------------------------------------------534 535 IF( kt == nit000 ) THEN536 IF(lwp) WRITE(numout,*)537 IF(lwp) WRITE(numout,*) 'dyn:hpg_wdj : hydrostatic pressure gradient trend'538 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ Weighted Density Jacobian'539 ENDIF540 541 ! Local constant initialization542 zcoef0 = - grav * 0.5_wp543 zalph = 0.5_wp - rn_gamma ! weighting coefficients (alpha=0.5-rn_gamma544 zbeta = 0.5_wp + rn_gamma ! (beta =1-alpha=0.5+rn_gamma545 546 ! Surface value (no ponderation)547 DO jj = 2, jpjm1548 DO ji = fs_2, fs_jpim1 ! vector opt.549 ! hydrostatic pressure gradient along s-surfaces550 zhpi(ji,jj,1) = zcoef0 / e1u(ji,jj) * ( fse3w(ji+1,jj ,1) * rhd(ji+1,jj ,1) &551 & - fse3w(ji ,jj ,1) * rhd(ji ,jj ,1) )552 zhpj(ji,jj,1) = zcoef0 / e2v(ji,jj) * ( fse3w(ji ,jj+1,1) * rhd(ji ,jj+1,1) &553 & - fse3w(ji ,jj ,1) * rhd(ji, jj ,1) )554 ! s-coordinate pressure gradient correction555 zuap = -zcoef0 * ( rhd (ji+1,jj,1) + rhd (ji,jj,1) ) &556 & * ( fsde3w(ji+1,jj,1) - fsde3w(ji,jj,1) ) / e1u(ji,jj)557 zvap = -zcoef0 * ( rhd (ji,jj+1,1) + rhd (ji,jj,1) ) &558 & * ( fsde3w(ji,jj+1,1) - fsde3w(ji,jj,1) ) / e2v(ji,jj)559 ! add to the general momentum trend560 ua(ji,jj,1) = ua(ji,jj,1) + zhpi(ji,jj,1) + zuap561 va(ji,jj,1) = va(ji,jj,1) + zhpj(ji,jj,1) + zvap562 END DO563 END DO564 565 ! Interior value (2=<jk=<jpkm1) (weighted with zalph & zbeta)566 DO jk = 2, jpkm1567 DO jj = 2, jpjm1568 DO ji = fs_2, fs_jpim1 ! vector opt.569 zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) + zcoef0 / e1u(ji,jj) &570 & * ( ( fsde3w(ji+1,jj,jk ) + fsde3w(ji,jj,jk ) &571 & - fsde3w(ji+1,jj,jk-1) - fsde3w(ji,jj,jk-1) ) &572 & * ( zalph * ( rhd (ji+1,jj,jk-1) - rhd (ji,jj,jk-1) ) &573 & + zbeta * ( rhd (ji+1,jj,jk ) - rhd (ji,jj,jk ) ) ) &574 & - ( rhd (ji+1,jj,jk ) + rhd (ji,jj,jk ) &575 & - rhd (ji+1,jj,jk-1) - rhd (ji,jj,jk-1) ) &576 & * ( zalph * ( fsde3w(ji+1,jj,jk-1) - fsde3w(ji,jj,jk-1) ) &577 & + zbeta * ( fsde3w(ji+1,jj,jk ) - fsde3w(ji,jj,jk ) ) ) )578 zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1) + zcoef0 / e2v(ji,jj) &579 & * ( ( fsde3w(ji,jj+1,jk ) + fsde3w(ji,jj,jk ) &580 & - fsde3w(ji,jj+1,jk-1) - fsde3w(ji,jj,jk-1) ) &581 & * ( zalph * ( rhd (ji,jj+1,jk-1) - rhd (ji,jj,jk-1) ) &582 & + zbeta * ( rhd (ji,jj+1,jk ) - rhd (ji,jj,jk ) ) ) &583 & - ( rhd (ji,jj+1,jk ) + rhd (ji,jj,jk ) &584 & - rhd (ji,jj+1,jk-1) - rhd (ji,jj,jk-1) ) &585 & * ( zalph * ( fsde3w(ji,jj+1,jk-1) - fsde3w(ji,jj,jk-1) ) &586 & + zbeta * ( fsde3w(ji,jj+1,jk ) - fsde3w(ji,jj,jk ) ) ) )587 ! add to the general momentum trend588 ua(ji,jj,jk) = ua(ji,jj,jk) + zhpi(ji,jj,jk)589 va(ji,jj,jk) = va(ji,jj,jk) + zhpj(ji,jj,jk)590 END DO591 END DO592 END DO593 !594 END SUBROUTINE hpg_wdj595 596 597 411 SUBROUTINE hpg_djc( kt ) 598 412 !!--------------------------------------------------------------------- … … 826 640 827 641 828 SUBROUTINE hpg_rot( kt )829 !!---------------------------------------------------------------------830 !! *** ROUTINE hpg_rot ***831 !!832 !! ** Method : rotated axes scheme (Thiem and Berntsen 2005)833 !!834 !! Reference: Thiem & Berntsen, Ocean Modelling, In press, 2005.835 !!----------------------------------------------------------------------836 USE wrk_nemo, ONLY: wrk_in_use, wrk_not_released837 USE oce , ONLY: zhpi => ta , zhpj => sa ! (ta,sa) used as 3D workspace838 USE wrk_nemo, ONLY: zdistr => wrk_2d_1 , zsina => wrk_2d_2 , zcosa => wrk_2d_3839 USE wrk_nemo, ONLY: zhpiorg => wrk_3d_1 , zhpirot => wrk_3d_2840 USE wrk_nemo, ONLY: zhpitra => wrk_3d_3 , zhpine => wrk_3d_4841 USE wrk_nemo, ONLY: zhpjorg => wrk_3d_5 , zhpjrot => wrk_3d_6842 USE wrk_nemo, ONLY: zhpjtra => wrk_3d_7 , zhpjne => wrk_3d_8843 !!844 INTEGER, INTENT(in) :: kt ! ocean time-step index845 !!846 INTEGER :: ji, jj, jk ! dummy loop indices847 REAL(wp) :: zforg, zcoef0, zuap, zmskd1, zmskd1m ! temporary scalar848 REAL(wp) :: zfrot , zvap, zmskd2, zmskd2m ! " "849 !!----------------------------------------------------------------------850 851 IF( wrk_in_use(2, 1,2,3) .OR. &852 wrk_in_use(3, 1,2,3,4,5,6,7,8) ) THEN853 CALL ctl_stop('dyn:hpg_rot: requested workspace arrays unavailable') ; RETURN854 ENDIF855 856 IF( kt == nit000 ) THEN857 IF(lwp) WRITE(numout,*)858 IF(lwp) WRITE(numout,*) 'dyn:hpg_rot : hydrostatic pressure gradient trend'859 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ s-coordinate case, rotated axes scheme used'860 ENDIF861 862 ! -------------------------------863 ! Local constant initialization864 ! -------------------------------865 zcoef0 = - grav * 0.5_wp866 zforg = 0.95_wp867 zfrot = 1._wp - zforg868 869 ! inverse of the distance between 2 diagonal T-points (defined at F-point) (here zcoef0/distance)870 zdistr(:,:) = zcoef0 / SQRT( e1f(:,:)*e1f(:,:) + e2f(:,:)*e1f(:,:) )871 872 ! sinus and cosinus of diagonal angle at F-point873 zsina(:,:) = ATAN2( e2f(:,:), e1f(:,:) )874 zcosa(:,:) = COS( zsina(:,:) )875 zsina(:,:) = SIN( zsina(:,:) )876 877 ! ---------------878 ! Surface value879 ! ---------------880 ! compute and add to the general trend the pressure gradients along the axes881 DO jj = 2, jpjm1882 DO ji = fs_2, fs_jpim1 ! vector opt.883 ! hydrostatic pressure gradient along s-surfaces884 zhpiorg(ji,jj,1) = zcoef0 / e1u(ji,jj) * ( fse3t(ji+1,jj,1) * rhd(ji+1,jj,1) &885 & - fse3t(ji ,jj,1) * rhd(ji ,jj,1) )886 zhpjorg(ji,jj,1) = zcoef0 / e2v(ji,jj) * ( fse3t(ji,jj+1,1) * rhd(ji,jj+1,1) &887 & - fse3t(ji,jj ,1) * rhd(ji,jj ,1) )888 ! s-coordinate pressure gradient correction889 zuap = -zcoef0 * ( rhd (ji+1,jj ,1) + rhd (ji,jj,1) ) &890 & * ( fsdept(ji+1,jj ,1) - fsdept(ji,jj,1) ) / e1u(ji,jj)891 zvap = -zcoef0 * ( rhd (ji ,jj+1,1) + rhd (ji,jj,1) ) &892 & * ( fsdept(ji ,jj+1,1) - fsdept(ji,jj,1) ) / e2v(ji,jj)893 ! add to the general momentum trend894 ua(ji,jj,1) = ua(ji,jj,1) + zforg * ( zhpiorg(ji,jj,1) + zuap )895 va(ji,jj,1) = va(ji,jj,1) + zforg * ( zhpjorg(ji,jj,1) + zvap )896 END DO897 END DO898 899 ! compute the pressure gradients in the diagonal directions900 DO jj = 1, jpjm1901 DO ji = 1, fs_jpim1 ! vector opt.902 zmskd1 = tmask(ji+1,jj+1,1) * tmask(ji ,jj,1) ! mask in the 1st diagnonal903 zmskd2 = tmask(ji ,jj+1,1) * tmask(ji+1,jj,1) ! mask in the 2nd diagnonal904 ! hydrostatic pressure gradient along s-surfaces905 zhpitra(ji,jj,1) = zdistr(ji,jj) * zmskd1 * ( fse3t(ji+1,jj+1,1) * rhd(ji+1,jj+1,1) &906 & - fse3t(ji ,jj ,1) * rhd(ji ,jj ,1) )907 zhpjtra(ji,jj,1) = zdistr(ji,jj) * zmskd2 * ( fse3t(ji ,jj+1,1) * rhd(ji ,jj+1,1) &908 & - fse3t(ji+1,jj ,1) * rhd(ji+1,jj ,1) )909 ! s-coordinate pressure gradient correction910 zuap = -zdistr(ji,jj) * zmskd1 * ( rhd (ji+1,jj+1,1) + rhd (ji ,jj,1) ) &911 & * ( fsdept(ji+1,jj+1,1) - fsdept(ji ,jj,1) )912 zvap = -zdistr(ji,jj) * zmskd2 * ( rhd (ji ,jj+1,1) + rhd (ji+1,jj,1) ) &913 & * ( fsdept(ji ,jj+1,1) - fsdept(ji+1,jj,1) )914 ! back rotation915 zhpine(ji,jj,1) = zcosa(ji,jj) * ( zhpitra(ji,jj,1) + zuap ) &916 & - zsina(ji,jj) * ( zhpjtra(ji,jj,1) + zvap )917 zhpjne(ji,jj,1) = zsina(ji,jj) * ( zhpitra(ji,jj,1) + zuap ) &918 & + zcosa(ji,jj) * ( zhpjtra(ji,jj,1) + zvap )919 END DO920 END DO921 922 ! interpolate and add to the general trend the diagonal gradient923 DO jj = 2, jpjm1924 DO ji = fs_2, fs_jpim1 ! vector opt.925 ! averaging926 zhpirot(ji,jj,1) = 0.5 * ( zhpine(ji,jj,1) + zhpine(ji ,jj-1,1) )927 zhpjrot(ji,jj,1) = 0.5 * ( zhpjne(ji,jj,1) + zhpjne(ji-1,jj ,1) )928 ! add to the general momentum trend929 ua(ji,jj,1) = ua(ji,jj,1) + zfrot * zhpirot(ji,jj,1)930 va(ji,jj,1) = va(ji,jj,1) + zfrot * zhpjrot(ji,jj,1)931 END DO932 END DO933 934 ! -----------------935 ! 2. interior value (2=<jk=<jpkm1)936 ! -----------------937 ! compute and add to the general trend the pressure gradients along the axes938 DO jk = 2, jpkm1939 DO jj = 2, jpjm1940 DO ji = fs_2, fs_jpim1 ! vector opt.941 ! hydrostatic pressure gradient along s-surfaces942 zhpiorg(ji,jj,jk) = zhpiorg(ji,jj,jk-1) &943 & + zcoef0 / e1u(ji,jj) * ( fse3t(ji+1,jj,jk ) * rhd(ji+1,jj,jk ) &944 & - fse3t(ji ,jj,jk ) * rhd(ji ,jj,jk ) &945 & + fse3t(ji+1,jj,jk-1) * rhd(ji+1,jj,jk-1) &946 & - fse3t(ji ,jj,jk-1) * rhd(ji ,jj,jk-1) )947 zhpjorg(ji,jj,jk) = zhpjorg(ji,jj,jk-1) &948 & + zcoef0 / e2v(ji,jj) * ( fse3t(ji,jj+1,jk ) * rhd(ji,jj+1,jk ) &949 & - fse3t(ji,jj ,jk ) * rhd(ji,jj, jk ) &950 & + fse3t(ji,jj+1,jk-1) * rhd(ji,jj+1,jk-1) &951 & - fse3t(ji,jj ,jk-1) * rhd(ji,jj, jk-1) )952 ! s-coordinate pressure gradient correction953 zuap = - zcoef0 * ( rhd (ji+1,jj ,jk) + rhd (ji,jj,jk) ) &954 & * ( fsdept(ji+1,jj ,jk) - fsdept(ji,jj,jk) ) / e1u(ji,jj)955 zvap = - zcoef0 * ( rhd (ji ,jj+1,jk) + rhd (ji,jj,jk) ) &956 & * ( fsdept(ji ,jj+1,jk) - fsdept(ji,jj,jk) ) / e2v(ji,jj)957 ! add to the general momentum trend958 ua(ji,jj,jk) = ua(ji,jj,jk) + zforg*( zhpiorg(ji,jj,jk) + zuap )959 va(ji,jj,jk) = va(ji,jj,jk) + zforg*( zhpjorg(ji,jj,jk) + zvap )960 END DO961 END DO962 END DO963 964 ! compute the pressure gradients in the diagonal directions965 DO jk = 2, jpkm1966 DO jj = 1, jpjm1967 DO ji = 1, fs_jpim1 ! vector opt.968 zmskd1 = tmask(ji+1,jj+1,jk ) * tmask(ji ,jj,jk ) ! level jk mask in the 1st diagnonal969 zmskd1m = tmask(ji+1,jj+1,jk-1) * tmask(ji ,jj,jk-1) ! level jk-1 " "970 zmskd2 = tmask(ji ,jj+1,jk ) * tmask(ji+1,jj,jk ) ! level jk mask in the 2nd diagnonal971 zmskd2m = tmask(ji ,jj+1,jk-1) * tmask(ji+1,jj,jk-1) ! level jk-1 " "972 ! hydrostatic pressure gradient along s-surfaces973 zhpitra(ji,jj,jk) = zhpitra(ji,jj,jk-1) &974 & + zdistr(ji,jj) * zmskd1 * ( fse3t(ji+1,jj+1,jk ) * rhd(ji+1,jj+1,jk) &975 & -fse3t(ji ,jj ,jk ) * rhd(ji ,jj ,jk) ) &976 & + zdistr(ji,jj) * zmskd1m * ( fse3t(ji+1,jj+1,jk-1) * rhd(ji+1,jj+1,jk-1) &977 & -fse3t(ji ,jj ,jk-1) * rhd(ji ,jj ,jk-1) )978 zhpjtra(ji,jj,jk) = zhpjtra(ji,jj,jk-1) &979 & + zdistr(ji,jj) * zmskd2 * ( fse3t(ji ,jj+1,jk ) * rhd(ji ,jj+1,jk) &980 & -fse3t(ji+1,jj ,jk ) * rhd(ji+1,jj, jk) ) &981 & + zdistr(ji,jj) * zmskd2m * ( fse3t(ji ,jj+1,jk-1) * rhd(ji ,jj+1,jk-1) &982 & -fse3t(ji+1,jj ,jk-1) * rhd(ji+1,jj, jk-1) )983 ! s-coordinate pressure gradient correction984 zuap = - zdistr(ji,jj) * zmskd1 * ( rhd (ji+1,jj+1,jk) + rhd (ji ,jj,jk) ) &985 & * ( fsdept(ji+1,jj+1,jk) - fsdept(ji ,jj,jk) )986 zvap = - zdistr(ji,jj) * zmskd2 * ( rhd (ji ,jj+1,jk) + rhd (ji+1,jj,jk) ) &987 & * ( fsdept(ji ,jj+1,jk) - fsdept(ji+1,jj,jk) )988 ! back rotation989 zhpine(ji,jj,jk) = zcosa(ji,jj) * ( zhpitra(ji,jj,jk) + zuap ) &990 & - zsina(ji,jj) * ( zhpjtra(ji,jj,jk) + zvap )991 zhpjne(ji,jj,jk) = zsina(ji,jj) * ( zhpitra(ji,jj,jk) + zuap ) &992 & + zcosa(ji,jj) * ( zhpjtra(ji,jj,jk) + zvap )993 END DO994 END DO995 END DO996 997 ! interpolate and add to the general trend998 DO jk = 2, jpkm1999 DO jj = 2, jpjm11000 DO ji = fs_2, fs_jpim1 ! vector opt.1001 ! averaging1002 zhpirot(ji,jj,jk) = 0.5 * ( zhpine(ji,jj,jk) + zhpine(ji ,jj-1,jk) )1003 zhpjrot(ji,jj,jk) = 0.5 * ( zhpjne(ji,jj,jk) + zhpjne(ji-1,jj ,jk) )1004 ! add to the general momentum trend1005 ua(ji,jj,jk) = ua(ji,jj,jk) + zfrot * zhpirot(ji,jj,jk)1006 va(ji,jj,jk) = va(ji,jj,jk) + zfrot * zhpjrot(ji,jj,jk)1007 END DO1008 END DO1009 END DO1010 !1011 IF( wrk_not_released(2, 1,2,3) .OR. &1012 wrk_not_released(3, 1,2,3,4,5,6,7,8) ) CALL ctl_stop('dyn:hpg_rot: failed to release workspace arrays')1013 !1014 END SUBROUTINE hpg_rot1015 1016 1017 642 SUBROUTINE hpg_prj( kt ) 1018 643 !!---------------------------------------------------------------------
Note: See TracChangeset
for help on using the changeset viewer.
