MODULE dynhpg !!====================================================================== !! *** MODULE dynhpg *** !! Ocean dynamics: hydrostatic pressure gradient trend !!====================================================================== !!---------------------------------------------------------------------- !! dyn_hpg : update the momentum trend with the horizontal !! gradient of the hydrostatic pressure !! !! default case : use of 3D work arrays (vector opt. available) !! key_s_coord : s-coordinate !! key_partial_steps : z-coordinate with partial steps !! default key : z-coordinate !!---------------------------------------------------------------------- !! * Modules used USE oce ! ocean dynamics and tracers USE dom_oce ! ocean space and time domain USE phycst ! physical constants USE in_out_manager ! I/O manager USE trddyn_oce ! dynamics trends diagnostics variables IMPLICIT NONE PRIVATE !! * Accessibility PUBLIC dyn_hpg ! routine called by step.F90 #if defined key_autotasking !!---------------------------------------------------------------------- !! 'key_autotasking' : j-k-i loop (j-slab) !!---------------------------------------------------------------------- LOGICAL, PUBLIC :: l_dyn_hpg_tsk = .TRUE. ! ??? LOGICAL, PUBLIC :: l_dyn_hpg = .FALSE. ! ??? #else !!---------------------------------------------------------------------- !! default case : k-j-i loop (vector opt.) !!---------------------------------------------------------------------- LOGICAL, PUBLIC :: l_dyn_hpg_tsk = .FALSE. ! ??? LOGICAL, PUBLIC :: l_dyn_hpg = .TRUE. ! ??? #endif !! * Substitutions # include "domzgr_substitute.h90" # include "vectopt_loop_substitute.h90" !!---------------------------------------------------------------------- !! OPA 9.0 , LODYC-IPSL (2003) !!---------------------------------------------------------------------- CONTAINS #if defined key_s_coord !!---------------------------------------------------------------------- !! 'key_s_coord' : s-coordinate !!---------------------------------------------------------------------- SUBROUTINE dyn_hpg( kt ) !!--------------------------------------------------------------------- !! *** ROUTINE dyn_hpg *** !! !! ** Purpose : Compute the now momentum trend due to the hor. gradient !! of the hydrostatic pressure. Add it to the general momentum trend. !! !! ** Method : The now hydrostatic pressure gradient at a given level !! jk is computed by taking the vertical integral of the in-situ !! density gradient along the model level from the suface to that !! level. s-coordinates ('key_s_coord'): a corrective term is added !! to the horizontal pressure gradient : !! zhpi = g ..... + 1/e1u mi(rhd) di[ g dep3w ] !! zhpj = g ..... + 1/e2v mj(rhd) dj[ g dep3w ] !! add it to the general momentum trend (ua,va). !! ua = ua - 1/e1u * zhpi !! va = va - 1/e2v * zhpj !! !! ** Action : - Update (ua,va) with the now hydrastatic pressure trend !! - Save the trend in (utrd,vtrd) ('key_trddyn') !! !! History : !! 1.0 ! 87-09 (P. Andrich, m.-a. Foujols) Original code !! ! 91-11 (G. Madec) !! ! 96-01 (G. Madec) s-coordinates !! ! 97-05 (G. Madec) split dynber into dynkeg and dynhpg !! 8.5 ! 02-08 (G. Madec) F90: Free form and module, vector opt. !!---------------------------------------------------------------------- !! * modules used USE oce, ONLY : zhpi => ta, & ! use ta as 3D workspace & zhpj => sa ! use sa as 3D workspace !! * Arguments INTEGER, INTENT( in ) :: kt ! ocean time-step index !! * Local declarations INTEGER :: ji, jj, jk ! dummy loop indices REAL(wp) :: & zcoef0, zcoef1, zuap, zvap ! temporary scalars !!---------------------------------------------------------------------- IF( kt == nit000 ) THEN IF(lwp) WRITE(numout,*) IF(lwp) WRITE(numout,*) 'dyn_hpg : hydrostatic pressure gradient trend' IF(lwp) WRITE(numout,*) '~~~~~~~ s-coordinate case, vector opt. case' ENDIF ! 0. Local constant initialization ! -------------------------------- zcoef0 = -g * 0.5 zuap = 0.e0 zvap = 0.e0 ! 1. Surface value ! ---------------- DO jj = 2, jpjm1 DO ji = fs_2, fs_jpim1 ! vector opt. ! hydrostatic pressure gradient along s-surfaces zhpi(ji,jj,1) = zcoef0 / e1u(ji,jj) & * ( fse3w(ji+1,jj,1) * rhd(ji+1,jj,1) - fse3w(ji,jj,1) * rhd(ji,jj,1) ) zhpj(ji,jj,1) = zcoef0 / e2v(ji,jj) & * ( fse3w(ji,jj+1,1) * rhd(ji,jj+1,1) - fse3w(ji,jj,1) * rhd(ji,jj,1) ) ! s-coordinate pressure gradient correction zuap = -zcoef0 * ( rhd(ji+1,jj,1) + rhd(ji,jj,1) ) & * ( fsde3w(ji+1,jj,1) - fsde3w(ji,jj,1) ) / e1u(ji,jj) zvap = -zcoef0 * ( rhd(ji,jj+1,1) + rhd(ji,jj,1) ) & * ( fsde3w(ji,jj+1,1) - fsde3w(ji,jj,1) ) / e2v(ji,jj) ! add to the general momentum trend ua(ji,jj,1) = ua(ji,jj,1) + zhpi(ji,jj,1) + zuap va(ji,jj,1) = va(ji,jj,1) + zhpj(ji,jj,1) + zvap # if defined key_trddyn ! save the trend for diagnostics utrd(ji,jj,1,1) = zhpi(ji,jj,1) + zuap vtrd(ji,jj,1,1) = zhpj(ji,jj,1) + zvap # endif END DO END DO ! 2. interior value (2= ta, & ! use ta as 3D workspace & zhpj => sa ! use sa as 3D workspace !! * Arguments INTEGER, INTENT( in ) :: kt ! ocean time-step index !! * local declarations INTEGER :: ji, jj, jk ! dummy loop indices INTEGER :: iku, ikv ! temporary integers REAL(wp) :: & zcoef0, zcoef1, zuap, & ! temporary scalars zcoef2, zcoef3, zvap ! " " !!---------------------------------------------------------------------- IF( kt == nit000 ) THEN IF(lwp) WRITE(numout,*) IF(lwp) WRITE(numout,*) 'dyn_hpg : hydrostatic pressure gradient trend' IF(lwp) WRITE(numout,*) '~~~~~~~ z-coordinate with partial steps' IF(lwp) WRITE(numout,*) ' vector optimization, no autotasking' ENDIF ! 0. Local constant initialization ! -------------------------------- zcoef0 = -g * 0.5 zuap = 0.e0 zvap = 0.e0 ! 1. Surface value ! ---------------- zcoef1 = zcoef0 * fse3w(1,1,1) DO jj = 2, jpjm1 DO ji = fs_2, fs_jpim1 ! vector opt. ! hydrostatic pressure gradient zhpi(ji,jj,1) = zcoef1 * ( rhd(ji+1,jj,1) - rhd(ji,jj,1) ) / e1u(ji,jj) zhpj(ji,jj,1) = zcoef1 * ( rhd(ji,jj+1,1) - rhd(ji,jj,1) ) / e2v(ji,jj) ! add to the general momentum trend ua(ji,jj,1) = ua(ji,jj,1) + zhpi(ji,jj,1) va(ji,jj,1) = va(ji,jj,1) + zhpj(ji,jj,1) # if defined key_trddyn ! save the momentum trends for diagnostics utrd(ji,jj,1,1) = zhpi(ji,jj,1) vtrd(ji,jj,1,1) = zhpj(ji,jj,1) # endif END DO END DO ! 2. interior value (2= 2 ) THEN ! subtract old value ua(ji,jj,iku) = ua(ji,jj,iku) - zhpi(ji,jj,iku) ! compute the new one zhpi (ji,jj,iku) = zhpi(ji,jj,iku-1) & + zcoef2 * ( rhd(ji+1,jj,iku-1) - rhd(ji,jj,iku-1) + gru(ji,jj) ) / e1u(ji,jj) ! add the new one to the general momentum trend ua(ji,jj,iku) = ua(ji,jj,iku) + zhpi(ji,jj,iku) # if defined key_trddyn ! save the momentum trends for diagnostics utrd(ji,jj,iku,1) = zhpi(ji,jj,iku) # endif ENDIF ! on j-direction IF ( ikv > 2 ) THEN ! subtract old value va(ji,jj,ikv) = va(ji,jj,ikv) - zhpj(ji,jj,ikv) ! compute the new one zhpj (ji,jj,ikv) = zhpj(ji,jj,ikv-1) & + zcoef3 * ( rhd(ji,jj+1,ikv-1) - rhd(ji,jj,ikv-1) + grv(ji,jj) ) / e2v(ji,jj) ! add the new one to the general momentum trend va(ji,jj,ikv) = va(ji,jj,ikv) + zhpj(ji,jj,ikv) # if defined key_trddyn ! save the momentum trends for diagnostics vtrd(ji,jj,ikv,1) = zhpj(ji,jj,ikv) # endif ENDIF # if ! defined key_vectopt_loop END DO # endif END DO IF( l_ctl .AND. lwp ) THEN ! print sum trends (used for debugging) zuap = SUM( ua(2:jpim1,2:jpjm1,1:jpkm1) * umask(2:jpim1,2:jpjm1,1:jpkm1) ) zvap = SUM( va(2:jpim1,2:jpjm1,1:jpkm1) * vmask(2:jpim1,2:jpjm1,1:jpkm1) ) WRITE(numout,*) ' hpg - Ua: ', zuap-u_ctl, ' Va: ', zvap-v_ctl u_ctl = zuap ; v_ctl = zvap ENDIF END SUBROUTINE dyn_hpg #else !!--------------------------------------------------------------------- !! Default case : z-coordinate !!--------------------------------------------------------------------- SUBROUTINE dyn_hpg( kt ) !!--------------------------------------------------------------------- !! *** ROUTINE dyn_hpg *** !! !! ** Purpose : Compute the now momentum trend due to the horizontal !! gradient of the hydrostatic pressure. Add it to the general !! momentum trend. !! !! ** Method : The now hydrostatic pressure gradient at a given level !! jk is computed by taking the vertical integral of the in-situ !! density gradient along the model level from the suface to that !! level: zhpi = g ..... !! zhpj = g ..... !! add it to the general momentum trend (ua,va). !! ua = ua - 1/e1u * zhpi !! va = va - 1/e2v * zhpj !! !! ** Action : - Update (ua,va) with the now hydrastatic pressure trend !! - Save the trend in (utrd,vtrd) ('key_trddyn') !! !! History : !! 1.0 ! 87-09 (P. Andrich, m.-a. Foujols) Original code !! ! 91-11 (G. Madec) !! ! 96-01 (G. Madec) s-coordinates !! ! 97-05 (G. Madec) split dynber into dynkeg and dynhpg !! 8.5 ! 02-07 (G. Madec) F90: Free form and module !!---------------------------------------------------------------------- !! * modules used USE oce, ONLY : zhpi => ta, & ! use ta as 3D workspace & zhpj => sa ! use sa as 3D workspace !! * Arguments INTEGER, INTENT( in ) :: kt ! ocean time-step index !! * local declarations INTEGER :: ji, jj, jk ! dummy loop indices REAL(wp) :: & zcoef0, zcoef1, zuap, zvap ! temporary scalars !!---------------------------------------------------------------------- IF( kt == nit000 ) THEN IF(lwp) WRITE(numout,*) IF(lwp) WRITE(numout,*) 'dyn_hpg : hydrostatic pressure gradient trend' IF(lwp) WRITE(numout,*) '~~~~~~~ z-coordinate case ' ENDIF ! 0. Local constant initialization ! -------------------------------- zcoef0 = -g * 0.5 zuap = 0.e0 zvap = 0.e0 ! 1. Surface value ! ---------------- zcoef1 = zcoef0 * fse3w(1,1,1) DO jj = 2, jpjm1 DO ji = fs_2, fs_jpim1 ! vector opt. ! hydrostatic pressure gradient zhpi(ji,jj,1) = zcoef1 * ( rhd(ji+1,jj,1) - rhd(ji,jj,1) ) / e1u(ji,jj) zhpj(ji,jj,1) = zcoef1 * ( rhd(ji,jj+1,1) - rhd(ji,jj,1) ) / e2v(ji,jj) ! add to the general momentum trend ua(ji,jj,1) = ua(ji,jj,1) + zhpi(ji,jj,1) va(ji,jj,1) = va(ji,jj,1) + zhpj(ji,jj,1) # if defined key_trddyn ! save the momentum trends for diagnostics utrd(ji,jj,1,1) = zhpi(ji,jj,1) vtrd(ji,jj,1,1) = zhpj(ji,jj,1) # endif END DO END DO ! 2. interior value (2=