MODULE dynhpgs !!====================================================================== !! *** MODULE dynhpgs *** !! Ocean dynamics: hydrostatic pressure gradient trend ; science routines !!====================================================================== !! History : OPA ! 1987-09 (P. Andrich, M.-A. Foujols) hpg_zco: Original code !! 5.0 ! 1991-11 (G. Madec) !! 7.0 ! 1996-01 (G. Madec) hpg_sco: Original code for s-coordinates !! 8.0 ! 1997-05 (G. Madec) split dynber into dynkeg and dynhpg !! 8.5 ! 2002-07 (G. Madec) F90: Free form and module !! 8.5 ! 2002-08 (A. Bozec) hpg_zps: Original code !! NEMO 1.0 ! 2005-10 (A. Beckmann, B.W. An) various s-coordinate options !! ! Original code for hpg_ctl, hpg_hel hpg_wdj, hpg_djc, hpg_rot !! - ! 2005-11 (G. Madec) style & small optimisation !! 3.3 ! 2010-10 (C. Ethe, G. Madec) reorganisation of initialisation phase !! 3.4 ! 2011-11 (H. Liu) hpg_prj: Original code for s-coordinates !! ! (A. Coward) suppression of hel, wdj and rot options !! 3.6 ! 2014-11 (P. Mathiot) hpg_isf: original code for ice shelf cavity !!---------------------------------------------------------------------- !!---------------------------------------------------------------------- !! hpg_zco : z-coordinate scheme !! hpg_zps : z-coordinate plus partial steps (interpolation) !! hpg_sco : s-coordinate (standard jacobian formulation) !!---------------------------------------------------------------------- ! USE par_kind USE len_oce USE in_out_manager IMPLICIT NONE !! * Substitutions # include "vectopt_loop_substitute.h90" !!---------------------------------------------------------------------- !! NEMO/OPA 3.3 , NEMO Consortium (2010) !! $Id: dynhpg.F90 9490 2018-04-23 08:44:07Z gm $ !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) !!---------------------------------------------------------------------- CONTAINS SUBROUTINE hpg_zco( grav, r1_e1u, r1_e2v, e3w_n, rhd, ua, va, zhpi, zhpj) !!--------------------------------------------------------------------- !! *** ROUTINE hpg_zco *** !! !! ** Method : z-coordinate case, levels are horizontal surfaces. !! 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 = grav ..... !! zhpj = grav ..... !! 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 !!---------------------------------------------------------------------- REAL(wp), INTENT(in) :: grav REAL(wp), INTENT(in), DIMENSION(jpi,jpj) :: r1_e1u, r1_e2v REAL(wp), INTENT(in), DIMENSION(jpi,jpj,jpk) :: e3w_n, rhd REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) :: ua, va REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) :: zhpi, zhpj REAL(wp) :: et !!---------------------------------------------------------------------- ! INTEGER :: ji, jj, jk ! dummy loop indices REAL(wp) :: zcoef0, zcoef1 ! temporary scalars !!---------------------------------------------------------------------- ! ! suppressed by MJB 27/10/2018 - should be moved to a control routine ! IF( kt == nit000 ) THEN ! IF(lwp) WRITE(numout,*) ! IF(lwp) WRITE(numout,*) 'dyn:hpg_zco : hydrostatic pressure gradient trend' ! IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ z-coordinate case ' ! ENDIF et = TIMER() zcoef0 = - grav * 0.5_wp ! Local constant initialization ! Surface value !$OMP PARALLEL DO PRIVATE(zcoef1) DO jj = 2, jpjm1 DO ji = fs_2, fs_jpim1 ! vector opt. zcoef1 = zcoef0 * e3w_n(ji,jj,1) ! hydrostatic pressure gradient zhpi(ji,jj,1) = zcoef1 * ( rhd(ji+1,jj,1) - rhd(ji,jj,1) ) * r1_e1u(ji,jj) zhpj(ji,jj,1) = zcoef1 * ( rhd(ji,jj+1,1) - rhd(ji,jj,1) ) * r1_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) END DO END DO ! ! interior value (2= 1 ) THEN ! on i-direction (level 2 or more) ua (ji,jj,iku) = ua(ji,jj,iku) - zhpi(ji,jj,iku) ! subtract old value zhpi(ji,jj,iku) = zhpi(ji,jj,iku-1) & ! compute the new one & + zcoef2 * ( rhd(ji+1,jj,iku-1) - rhd(ji,jj,iku-1) + gru(ji,jj) ) * r1_e1u(ji,jj) ua (ji,jj,iku) = ua(ji,jj,iku) + zhpi(ji,jj,iku) ! add the new one to the general momentum trend ENDIF IF( ikv > 1 ) THEN ! on j-direction (level 2 or more) va (ji,jj,ikv) = va(ji,jj,ikv) - zhpj(ji,jj,ikv) ! subtract old value zhpj(ji,jj,ikv) = zhpj(ji,jj,ikv-1) & ! compute the new one & + zcoef3 * ( rhd(ji,jj+1,ikv-1) - rhd(ji,jj,ikv-1) + grv(ji,jj) ) * r1_e2v(ji,jj) va (ji,jj,ikv) = va(ji,jj,ikv) + zhpj(ji,jj,ikv) ! add the new one to the general momentum trend ENDIF END DO END DO ! !hpg_zps_time = hpg_zps_time + (TIMER() - et) ! Moved timer up call tree END SUBROUTINE hpg_zps SUBROUTINE hpg_sco( grav, r1_e1u, r1_e2v, e3w_n, rhd, gde3w_n, ln_linssh, ua, va, zhpi, zhpj) !!--------------------------------------------------------------------- !! *** ROUTINE hpg_sco *** !! !! ** Method : s-coordinate case. Jacobian scheme. !! 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 (ln_sco): a corrective term is added !! to the horizontal pressure gradient : !! zhpi = grav ..... + 1/e1u mi(rhd) di[ grav dep3w ] !! zhpj = grav ..... + 1/e2v mj(rhd) dj[ grav 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 !!---------------------------------------------------------------------- REAL(wp), INTENT(in):: grav REAL(wp), INTENT(in), DIMENSION(jpi,jpj) :: r1_e1u, r1_e2v REAL(wp), INTENT(in), DIMENSION(jpi,jpj,jpk) :: e3w_n, rhd, gde3w_n LOGICAL, INTENT(in):: ln_linssh REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) :: ua, va REAL(wp), INTENT(INOUT), DIMENSION(jpi,jpj,jpk) :: zhpi, zhpj !!---------------------------------------------------------------------- !! INTEGER :: ji, jj, jk, jii, jjj ! dummy loop indices REAL(wp) :: zcoef0, zuap, zvap, znad, ztmp ! temporary scalars LOGICAL :: ll_tmp1, ll_tmp2 ! local logical variables REAL(wp), DIMENSION(:,:), ALLOCATABLE :: zcpx, zcpy !W/D pressure filter REAL(wp) :: et !!---------------------------------------------------------------------- ! ! suppressed by MJB 27/10/2018 - should be moved to a control routine ! IF( kt == nit000 ) THEN ! IF(lwp) WRITE(numout,*) ! IF(lwp) WRITE(numout,*) 'dyn:hpg_sco : hydrostatic pressure gradient trend' ! IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ s-coordinate case, OPA original scheme used' ! ENDIF ! et = TIMER() zcoef0 = - grav * 0.5_wp IF ( ln_linssh ) THEN ; znad = 0._wp ! Fixed volume: density anomaly ELSE ; znad = 1._wp ! Variable volume: density ENDIF ! ! Surface value !$OMP PARALLEL !$OMP DO PRIVATE(zuap,zvap) DO jj = 2, jpjm1 DO ji = fs_2, fs_jpim1 ! vector opt. ! hydrostatic pressure gradient along s-surfaces zhpi(ji,jj,1) = zcoef0 * ( e3w_n(ji+1,jj ,1) * ( znad + rhd(ji+1,jj ,1) ) & & - e3w_n(ji ,jj ,1) * ( znad + rhd(ji ,jj ,1) ) ) * r1_e1u(ji,jj) zhpj(ji,jj,1) = zcoef0 * ( e3w_n(ji ,jj+1,1) * ( znad + rhd(ji ,jj+1,1) ) & & - e3w_n(ji ,jj ,1) * ( znad + rhd(ji ,jj ,1) ) ) * r1_e2v(ji,jj) ! s-coordinate pressure gradient correction zuap = -zcoef0 * ( rhd (ji+1,jj,1) + rhd (ji,jj,1) + 2._wp * znad ) & & * ( gde3w_n(ji+1,jj,1) - gde3w_n(ji,jj,1) ) * r1_e1u(ji,jj) zvap = -zcoef0 * ( rhd (ji,jj+1,1) + rhd (ji,jj,1) + 2._wp * znad ) & & * ( gde3w_n(ji,jj+1,1) - gde3w_n(ji,jj,1) ) * r1_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 END DO END DO ! interior value (2=