New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
Changeset 14060 – NEMO

Changeset 14060


Ignore:
Timestamp:
2020-12-03T17:03:42+01:00 (4 years ago)
Author:
ayoung
Message:

Merging ticket #2480 into trunk. (Revised djc scheme and two new namelist variables)

Location:
NEMO/trunk
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • NEMO/trunk/cfgs/SHARED/namelist_ref

    r14056 r14060  
    10131013   ln_hpg_isf  = .false.   !  s-coordinate (sco ) adapted to isf 
    10141014   ln_hpg_djc  = .false.   !  s-coordinate (Density Jacobian with Cubic polynomial) 
     1015      ln_hpg_djc_vnh = .true.  !  hor.  bc type for djc scheme (T=von Neumann, F=linear extrapolation) 
     1016      ln_hpg_djc_vnv = .true.  !  vert. bc type for djc scheme (T=von Neumann, F=linear extrapolation) 
    10151017   ln_hpg_prj  = .false.   !  s-coordinate (Pressure Jacobian scheme) 
    10161018/ 
  • NEMO/trunk/src/OCE/DYN/dynhpg.F90

    r13982 r14060  
    1717   !!                 !           (A. Coward) suppression of hel, wdj and rot options 
    1818   !!            3.6  !  2014-11  (P. Mathiot) hpg_isf: original code for ice shelf cavity 
     19   !!            4.2  !  2020-12  (M. Bell, A. Young) hpg_djc: revised djc scheme 
    1920   !!---------------------------------------------------------------------- 
    2021 
     
    7273   INTEGER, PARAMETER ::   np_isf    =  5   ! s-coordinate similar to sco modify for isf 
    7374   ! 
    74    INTEGER, PUBLIC ::   nhpg         !: type of pressure gradient scheme used ! (deduced from ln_hpg_... flags) (PUBLIC for TAM) 
     75   INTEGER, PUBLIC  ::   nhpg         !: type of pressure gradient scheme used ! (deduced from ln_hpg_... flags) (PUBLIC for TAM) 
     76   ! 
     77   LOGICAL          ::   ln_hpg_djc_vnh, ln_hpg_djc_vnv                 ! flag to specify hpg_djc boundary condition type 
     78   REAL(wp), PUBLIC ::   aco_bc_hor, bco_bc_hor, aco_bc_vrt, bco_bc_vrt !: coefficients for hpg_djc hor and vert boundary conditions 
    7579 
    7680   !! * Substitutions 
     
    156160      !! 
    157161      NAMELIST/namdyn_hpg/ ln_hpg_zco, ln_hpg_zps, ln_hpg_sco,     & 
    158          &                 ln_hpg_djc, ln_hpg_prj, ln_hpg_isf 
     162         &                 ln_hpg_djc, ln_hpg_prj, ln_hpg_isf,     & 
     163         &                 ln_hpg_djc_vnh, ln_hpg_djc_vnv 
    159164      !!---------------------------------------------------------------------- 
    160165      ! 
     
    179184      ENDIF 
    180185      ! 
    181       IF( ln_hpg_djc )   & 
    182          &   CALL ctl_stop('dyn_hpg_init : Density Jacobian: Cubic polynominal method',   & 
    183          &                 '   currently disabled (bugs under investigation).'        ,   & 
    184          &                 '   Please select either  ln_hpg_sco or  ln_hpg_prj instead' ) 
    185          ! 
    186       IF( .NOT.ln_linssh .AND. .NOT.(ln_hpg_sco.OR.ln_hpg_prj.OR.ln_hpg_isf) )          & 
     186      IF( .NOT.ln_linssh .AND. .NOT.(ln_hpg_sco.OR.ln_hpg_prj.OR.ln_hpg_isf.OR.ln_hpg_djc) )          & 
    187187         &   CALL ctl_stop('dyn_hpg_init : non-linear free surface requires either ',   & 
    188188         &                 '   the standard jacobian formulation hpg_sco    or '    ,   & 
     
    220220      ENDIF 
    221221      !                           
     222      IF ( ln_hpg_djc ) THEN 
     223         IF (ln_hpg_djc_vnh) THEN ! Von Neumann boundary condition 
     224           IF(lwp) WRITE(numout,*) '           horizontal bc: von Neumann ' 
     225           aco_bc_hor = 6.0_wp/5.0_wp 
     226           bco_bc_hor = 7.0_wp/15.0_wp 
     227         ELSE ! Linear extrapolation 
     228           IF(lwp) WRITE(numout,*) '           horizontal bc: linear extrapolation' 
     229           aco_bc_hor = 3.0_wp/2.0_wp 
     230           bco_bc_hor = 1.0_wp/2.0_wp 
     231         END IF 
     232         IF (ln_hpg_djc_vnv) THEN ! Von Neumann boundary condition 
     233           IF(lwp) WRITE(numout,*) '           vertical bc: von Neumann ' 
     234           aco_bc_vrt = 6.0_wp/5.0_wp 
     235           bco_bc_vrt = 7.0_wp/15.0_wp 
     236         ELSE ! Linear extrapolation 
     237           IF(lwp) WRITE(numout,*) '           vertical bc: linear extrapolation' 
     238           aco_bc_vrt = 3.0_wp/2.0_wp 
     239           bco_bc_vrt = 1.0_wp/2.0_wp 
     240         END IF 
     241      END IF 
    222242   END SUBROUTINE dyn_hpg_init 
    223243 
     
    631651      !! 
    632652      INTEGER  ::   ji, jj, jk          ! dummy loop indices 
     653      INTEGER  ::   iktb, iktt          ! jk indices at tracer points for top and bottom points  
    633654      REAL(wp) ::   zcoef0, zep, cffw   ! temporary scalars 
    634       REAL(wp) ::   z1_10, cffu, cffx   !    "         " 
    635       REAL(wp) ::   z1_12, cffv, cffy   !    "         " 
     655      REAL(wp) ::   z_grav_10, z1_12 
     656      REAL(wp) ::   cffu, cffx          !    "         " 
     657      REAL(wp) ::   cffv, cffy          !    "         " 
    636658      LOGICAL  ::   ll_tmp1, ll_tmp2    ! local logical variables 
    637659      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zhpi, zhpj 
    638       REAL(wp), DIMENSION(jpi,jpj,jpk) ::   dzx, dzy, dzz, dzu, dzv, dzw 
    639       REAL(wp), DIMENSION(jpi,jpj,jpk) ::   drhox, drhoy, drhoz, drhou, drhov, drhow 
    640       REAL(wp), DIMENSION(jpi,jpj,jpk) ::   rho_i, rho_j, rho_k 
     660      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   rhd_opt 
     661  
     662      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zdzx, zdzy, zdzz                          ! Primitive grid differences ('delta_xyz') 
     663      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zdz_i, zdz_j, zdz_k                       ! Harmonic average of primitive grid differences ('d_xyz') 
     664      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zdrhox, zdrhoy, zdrhoz                    ! Primitive rho differences ('delta_rho') 
     665      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zdrho_i, zdrho_j, zdrho_k                 ! Harmonic average of primitive rho differences ('d_rho') 
     666      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   z_rho_i, z_rho_j, z_rho_k                 ! Face intergrals 
     667      REAL(wp), DIMENSION(jpi,jpj)     ::   zz_dz_i, zz_dz_j, zz_drho_i, zz_drho_j    ! temporary arrays  
    641668      REAL(wp), DIMENSION(:,:), ALLOCATABLE ::   zcpx, zcpy   !W/D pressure filter 
    642669      !!---------------------------------------------------------------------- 
     
    691718      ! Local constant initialization 
    692719      zcoef0 = - grav * 0.5_wp 
    693       z1_10  = 1._wp / 10._wp 
    694       z1_12  = 1._wp / 12._wp 
     720      z_grav_10  = grav / 10._wp 
     721      z1_12  = 1.0_wp / 12._wp 
     722 
     723 
     724!! To be removed once combined with KERNEL08.  KERNEL08 action removes znad from the hpg module, spg is calculated in spg moduel instead.  The final code will just have rhd(:,:,:), no rhd_opt(:,:,:) 
     725      IF( .NOT. ln_linssh ) THEN 
     726         rhd_opt(:,:,:) = rhd(:,:,:) + 1.0_wp ! for vvl option 
     727      ELSE 
     728         rhd_opt(:,:,:) = rhd(:,:,:) 
     729      END IF 
    695730 
    696731      !---------------------------------------------------------------------------------------- 
    697       !  compute and store in provisional arrays elementary vertical and horizontal differences 
     732      !  1. compute and store elementary vertical differences in provisional arrays  
    698733      !---------------------------------------------------------------------------------------- 
    699734 
    700 !!bug gm   Not a true bug, but... dzz=e3w  for dzx, dzy verify what it is really 
    701  
    702       DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 
    703          drhoz(ji,jj,jk) = rhd    (ji  ,jj  ,jk) - rhd    (ji,jj,jk-1) 
    704          dzz  (ji,jj,jk) = gde3w(ji  ,jj  ,jk) - gde3w(ji,jj,jk-1) 
    705          drhox(ji,jj,jk) = rhd    (ji+1,jj  ,jk) - rhd    (ji,jj,jk  ) 
    706          dzx  (ji,jj,jk) = gde3w(ji+1,jj  ,jk) - gde3w(ji,jj,jk  ) 
    707          drhoy(ji,jj,jk) = rhd    (ji  ,jj+1,jk) - rhd    (ji,jj,jk  ) 
    708          dzy  (ji,jj,jk) = gde3w(ji  ,jj+1,jk) - gde3w(ji,jj,jk  ) 
     735!!bug gm   Not a true bug, but... zdzz=e3w  for zdzx, zdzy verify what it is really 
     736 
     737      DO_3D( 1, 1, 1, 1, 2, jpkm1 )   
     738         zdrhoz(ji,jj,jk) =   rhd_opt    (ji  ,jj  ,jk) - rhd_opt    (ji,jj,jk-1) 
     739         zdzz  (ji,jj,jk) = - gde3w(ji  ,jj  ,jk) + gde3w(ji,jj,jk-1) 
    709740      END_3D 
    710741 
    711742      !------------------------------------------------------------------------- 
    712       ! compute harmonic averages using eq. 5.18 
     743      ! 2. compute harmonic averages for vertical differences using eq. 5.18 
    713744      !------------------------------------------------------------------------- 
    714745      zep = 1.e-15 
    715746 
    716 !!bug  gm  drhoz not defined at level 1 and used (jk-1 with jk=2) 
    717 !!bug  gm  idem for drhox, drhoy et ji=jpi and jj=jpj 
    718  
    719       DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 
    720          cffw = 2._wp * drhoz(ji  ,jj  ,jk) * drhoz(ji,jj,jk-1) 
    721  
    722          cffu = 2._wp * drhox(ji+1,jj  ,jk) * drhox(ji,jj,jk  ) 
    723          cffx = 2._wp * dzx  (ji+1,jj  ,jk) * dzx  (ji,jj,jk  ) 
    724  
    725          cffv = 2._wp * drhoy(ji  ,jj+1,jk) * drhoy(ji,jj,jk  ) 
    726          cffy = 2._wp * dzy  (ji  ,jj+1,jk) * dzy  (ji,jj,jk  ) 
    727  
     747!! mb zdrho_k, zdz_k, zdrho_i, zdz_i, zdrho_j, zdz_j re-centred about the point (ji,jj,jk)  
     748      zdrho_k(:,:,:) = 0._wp 
     749      zdz_k  (:,:,:) = 0._wp 
     750 
     751      DO_3D( 1, 1, 1, 1, 2, jpk-2 )  
     752         cffw = 2._wp * zdrhoz(ji  ,jj  ,jk) * zdrhoz(ji,jj,jk+1) 
    728753         IF( cffw > zep) THEN 
    729             drhow(ji,jj,jk) = 2._wp *   drhoz(ji,jj,jk) * drhoz(ji,jj,jk-1)   & 
    730                &                    / ( drhoz(ji,jj,jk) + drhoz(ji,jj,jk-1) ) 
     754            zdrho_k(ji,jj,jk) = cffw / ( zdrhoz(ji,jj,jk) + zdrhoz(ji,jj,jk+1) ) 
     755         ENDIF 
     756         zdz_k(ji,jj,jk) = 2._wp *   zdzz(ji,jj,jk) * zdzz(ji,jj,jk+1)   & 
     757            &                  / ( zdzz(ji,jj,jk) + zdzz(ji,jj,jk+1) ) 
     758      END_3D 
     759 
     760      !---------------------------------------------------------------------------------- 
     761      ! 3. apply boundary conditions at top and bottom using 5.36-5.37 
     762      !---------------------------------------------------------------------------------- 
     763 
     764! mb for sea-ice shelves we will need to re-write this upper boundary condition in the same form as the lower boundary condition 
     765      zdrho_k(:,:,1) = aco_bc_vrt * ( rhd_opt    (:,:,2) - rhd_opt    (:,:,1) ) - bco_bc_vrt * zdrho_k(:,:,2) 
     766      zdz_k  (:,:,1) = aco_bc_vrt * (-gde3w(:,:,2) + gde3w(:,:,1) ) - bco_bc_vrt * zdz_k  (:,:,2) 
     767 
     768      DO_2D( 1, 1, 1, 1 ) 
     769         IF ( mbkt(ji,jj)>1 ) THEN 
     770            iktb = mbkt(ji,jj) 
     771            zdrho_k(ji,jj,iktb) = aco_bc_vrt * (     rhd_opt(ji,jj,iktb) -     rhd_opt(ji,jj,iktb-1) ) - bco_bc_vrt * zdrho_k(ji,jj,iktb-1) 
     772            zdz_k  (ji,jj,iktb) = aco_bc_vrt * (-gde3w(ji,jj,iktb) + gde3w(ji,jj,iktb-1) ) - bco_bc_vrt * zdz_k  (ji,jj,iktb-1)  
     773         END IF 
     774      END_2D 
     775 
     776      !-------------------------------------------------------------- 
     777      ! 4. Compute side face integrals 
     778      !------------------------------------------------------------- 
     779 
     780!! ssh replaces e3w_n ; gde3w is a depth; the formulae involve heights   
     781!! rho_k stores grav * FX / rho_0   
     782 
     783      !-------------------------------------------------------------- 
     784      ! 4. a) Upper half of top-most grid box, compute and store 
     785      !------------------------------------------------------------- 
     786! *** AY note: ssh(ji,jj,Kmm) + gde3w(ji,jj,1) = e3w(ji,jj,1) 
     787      DO_2D( 0, 1, 0, 1) 
     788         z_rho_k(ji,jj,1) =  grav * ( ssh(ji,jj,Kmm) + gde3w(ji,jj,1) )                        &  
     789            &                     * (  rhd_opt(ji,jj,1)                                        & 
     790            &                     + 0.5_wp * (   rhd_opt    (ji,jj,2) - rhd_opt    (ji,jj,1) ) & 
     791            &                              * (   ssh   (ji,jj,Kmm) + gde3w(ji,jj,1) )          & 
     792            &                              / ( - gde3w(ji,jj,2) + gde3w(ji,jj,1) )  ) 
     793      END_2D 
     794 
     795      !-------------------------------------------------------------- 
     796      ! 4. b) Interior faces, compute and store 
     797      !------------------------------------------------------------- 
     798 
     799      DO_3D( 0, 1, 0, 1, 2, jpkm1 ) 
     800         z_rho_k(ji,jj,jk) = zcoef0 * (   rhd_opt    (ji,jj,jk) + rhd_opt    (ji,jj,jk-1) )                                   & 
     801            &                       * ( - gde3w(ji,jj,jk) + gde3w(ji,jj,jk-1) )                                               & 
     802            &                       + z_grav_10 * (                                                                           & 
     803            &     (   zdrho_k  (ji,jj,jk) - zdrho_k  (ji,jj,jk-1) )                                                           & 
     804            &   * ( - gde3w(ji,jj,jk) + gde3w(ji,jj,jk-1) - z1_12 * ( zdz_k  (ji,jj,jk) + zdz_k  (ji,jj,jk-1) ) )             & 
     805            &   - ( zdz_k    (ji,jj,jk) - zdz_k    (ji,jj,jk-1) )                                                             & 
     806            &   * ( rhd_opt    (ji,jj,jk) - rhd_opt    (ji,jj,jk-1) - z1_12 * ( zdrho_k(ji,jj,jk) + zdrho_k(ji,jj,jk-1) ) )   & 
     807            &                             ) 
     808      END_3D 
     809 
     810      !---------------------------------------------------------------------------------------- 
     811      !  5. compute and store elementary horizontal differences in provisional arrays  
     812      !---------------------------------------------------------------------------------------- 
     813 
     814      DO_3D( 1, 0, 1, 0, 1, jpkm1 ) 
     815         zdrhox(ji,jj,jk) =   rhd_opt    (ji+1,jj  ,jk) - rhd_opt    (ji,jj,jk  ) 
     816         zdzx  (ji,jj,jk) = - gde3w(ji+1,jj  ,jk) + gde3w(ji,jj,jk  ) 
     817         zdrhoy(ji,jj,jk) =   rhd_opt    (ji  ,jj+1,jk) - rhd_opt    (ji,jj,jk  ) 
     818         zdzy  (ji,jj,jk) = - gde3w(ji  ,jj+1,jk) + gde3w(ji,jj,jk  ) 
     819      END_3D 
     820 
     821      CALL lbc_lnk_multi( 'dynhpg', zdrhox, 'U', 1., zdzx, 'U', 1., zdrhoy, 'V', 1., zdzy, 'V', 1. )  
     822 
     823      !------------------------------------------------------------------------- 
     824      ! 6. compute harmonic averages using eq. 5.18 
     825      !------------------------------------------------------------------------- 
     826 
     827      DO_3D( 0, 1, 0, 1, 1, jpkm1 ) 
     828         cffu = 2._wp * zdrhox(ji-1,jj  ,jk) * zdrhox(ji,jj,jk  ) 
     829         IF( cffu > zep ) THEN 
     830            zdrho_i(ji,jj,jk) = cffu / ( zdrhox(ji-1,jj,jk) + zdrhox(ji,jj,jk) ) 
    731831         ELSE 
    732             drhow(ji,jj,jk) = 0._wp 
    733          ENDIF 
    734  
    735          dzw(ji,jj,jk) = 2._wp *   dzz(ji,jj,jk) * dzz(ji,jj,jk-1)   & 
    736             &                  / ( dzz(ji,jj,jk) + dzz(ji,jj,jk-1) ) 
    737  
    738          IF( cffu > zep ) THEN 
    739             drhou(ji,jj,jk) = 2._wp *   drhox(ji+1,jj,jk) * drhox(ji,jj,jk)   & 
    740                &                    / ( drhox(ji+1,jj,jk) + drhox(ji,jj,jk) ) 
     832            zdrho_i(ji,jj,jk ) = 0._wp 
     833         ENDIF 
     834 
     835         cffx = 2._wp * zdzx  (ji-1,jj  ,jk) * zdzx  (ji,jj,jk  ) 
     836         IF( cffx > zep ) THEN 
     837            zdz_i(ji,jj,jk) = cffx / ( zdzx(ji-1,jj,jk) + zdzx(ji,jj,jk) ) 
    741838         ELSE 
    742             drhou(ji,jj,jk ) = 0._wp 
    743          ENDIF 
    744  
    745          IF( cffx > zep ) THEN 
    746             dzu(ji,jj,jk) = 2._wp *   dzx(ji+1,jj,jk) * dzx(ji,jj,jk)   & 
    747                &                  / ( dzx(ji+1,jj,jk) + dzx(ji,jj,jk) ) 
     839            zdz_i(ji,jj,jk) = 0._wp 
     840         ENDIF 
     841 
     842         cffv = 2._wp * zdrhoy(ji  ,jj-1,jk) * zdrhoy(ji,jj,jk  ) 
     843         IF( cffv > zep ) THEN 
     844            zdrho_j(ji,jj,jk) = cffv / ( zdrhoy(ji,jj-1,jk) + zdrhoy(ji,jj,jk) ) 
    748845         ELSE 
    749             dzu(ji,jj,jk) = 0._wp 
    750          ENDIF 
    751  
    752          IF( cffv > zep ) THEN 
    753             drhov(ji,jj,jk) = 2._wp *   drhoy(ji,jj+1,jk) * drhoy(ji,jj,jk)   & 
    754                &                    / ( drhoy(ji,jj+1,jk) + drhoy(ji,jj,jk) ) 
     846            zdrho_j(ji,jj,jk) = 0._wp 
     847         ENDIF 
     848 
     849         cffy = 2._wp * zdzy  (ji  ,jj-1,jk) * zdzy  (ji,jj,jk  ) 
     850         IF( cffy > zep ) THEN 
     851            zdz_j(ji,jj,jk) = cffy / ( zdzy(ji,jj-1,jk) + zdzy(ji,jj,jk) ) 
    755852         ELSE 
    756             drhov(ji,jj,jk) = 0._wp 
    757          ENDIF 
    758  
    759          IF( cffy > zep ) THEN 
    760             dzv(ji,jj,jk) = 2._wp *   dzy(ji,jj+1,jk) * dzy(ji,jj,jk)   & 
    761                &                  / ( dzy(ji,jj+1,jk) + dzy(ji,jj,jk) ) 
    762          ELSE 
    763             dzv(ji,jj,jk) = 0._wp 
    764          ENDIF 
    765  
    766       END_3D 
     853            zdz_j(ji,jj,jk) = 0._wp 
     854         ENDIF 
     855      END_3D 
     856       
     857!!! Note that zdzx, zdzy, zdzz, zdrhox, zdrhoy and zdrhoz should NOT be used beyond this point       
    767858 
    768859      !---------------------------------------------------------------------------------- 
    769       ! apply boundary conditions at top and bottom using 5.36-5.37 
     860      ! 6B. apply boundary conditions at side boundaries using 5.36-5.37 
    770861      !---------------------------------------------------------------------------------- 
    771       drhow(:,:, 1 ) = 1.5_wp * ( drhoz(:,:, 2 ) - drhoz(:,:,  1  ) ) - 0.5_wp * drhow(:,:,  2  ) 
    772       drhou(:,:, 1 ) = 1.5_wp * ( drhox(:,:, 2 ) - drhox(:,:,  1  ) ) - 0.5_wp * drhou(:,:,  2  ) 
    773       drhov(:,:, 1 ) = 1.5_wp * ( drhoy(:,:, 2 ) - drhoy(:,:,  1  ) ) - 0.5_wp * drhov(:,:,  2  ) 
    774  
    775       drhow(:,:,jpk) = 1.5_wp * ( drhoz(:,:,jpk) - drhoz(:,:,jpkm1) ) - 0.5_wp * drhow(:,:,jpkm1) 
    776       drhou(:,:,jpk) = 1.5_wp * ( drhox(:,:,jpk) - drhox(:,:,jpkm1) ) - 0.5_wp * drhou(:,:,jpkm1) 
    777       drhov(:,:,jpk) = 1.5_wp * ( drhoy(:,:,jpk) - drhoy(:,:,jpkm1) ) - 0.5_wp * drhov(:,:,jpkm1) 
    778  
     862 
     863      DO jk = 1, jpkm1 
     864         zz_drho_i(:,:) = zdrho_i(:,:,jk) 
     865         zz_dz_i  (:,:) = zdz_i  (:,:,jk) 
     866         zz_drho_j(:,:) = zdrho_j(:,:,jk) 
     867         zz_dz_j  (:,:) = zdz_j  (:,:,jk) 
     868         DO_2D( 0, 1, 0, 1) 
     869            ! Walls coming from left: should check from 2 to jpi-1 (and jpj=2-jpj) 
     870            IF (ji < jpi) THEN 
     871               IF ( umask(ji,jj,jk) > 0.5_wp .AND. umask(ji-1,jj,jk) < 0.5_wp .AND. umask(ji+1,jj,jk) > 0.5_wp)  THEN   
     872                  zz_drho_i(ji,jj) = aco_bc_hor * ( rhd_opt    (ji+1,jj,jk) - rhd_opt    (ji,jj,jk) ) - bco_bc_hor * zdrho_i(ji+1,jj,jk)  
     873                  zz_dz_i  (ji,jj) = aco_bc_hor * (-gde3w(ji+1,jj,jk) + gde3w(ji,jj,jk) ) - bco_bc_hor * zdz_i  (ji+1,jj,jk) 
     874               END IF 
     875            END IF 
     876            ! Walls coming from right: should check from 3 to jpi (and jpj=2-jpj) 
     877            IF (ji > 2) THEN 
     878               IF ( umask(ji,jj,jk) < 0.5_wp .AND. umask(ji-1,jj,jk) > 0.5_wp .AND. umask(ji-2,jj,jk) > 0.5_wp) THEN 
     879                  zz_drho_i(ji,jj) = aco_bc_hor * ( rhd_opt    (ji,jj,jk) - rhd_opt    (ji-1,jj,jk) ) - bco_bc_hor * zdrho_i(ji-1,jj,jk)   
     880                  zz_dz_i  (ji,jj) = aco_bc_hor * (-gde3w(ji,jj,jk) + gde3w(ji-1,jj,jk) ) - bco_bc_hor * zdz_i  (ji-1,jj,jk) 
     881               END IF 
     882            END IF 
     883            ! Walls coming from left: should check from 2 to jpj-1 (and jpi=2-jpi) 
     884            IF (jj < jpj) THEN 
     885               IF ( vmask(ji,jj,jk) > 0.5_wp .AND. vmask(ji,jj-1,jk) < 0.5_wp .AND. vmask(ji,jj+1,jk) > 0.5_wp)  THEN 
     886                  zz_drho_j(ji,jj) = aco_bc_hor * ( rhd_opt    (ji,jj+1,jk) - rhd_opt    (ji,jj,jk) ) - bco_bc_hor * zdrho_j(ji,jj+1,jk) 
     887                  zz_dz_j  (ji,jj) = aco_bc_hor * (-gde3w(ji,jj+1,jk) + gde3w(ji,jj,jk) ) - bco_bc_hor * zdz_j  (ji,jj+1,jk) 
     888               END IF 
     889            END IF  
     890            ! Walls coming from right: should check from 3 to jpj (and jpi=2-jpi) 
     891            IF (jj > 2) THEN 
     892               IF ( vmask(ji,jj,jk) < 0.5_wp .AND. vmask(ji,jj-1,jk) > 0.5_wp .AND. vmask(ji,jj-2,jk) > 0.5_wp) THEN  
     893                  zz_drho_j(ji,jj) = aco_bc_hor * ( rhd_opt    (ji,jj,jk) - rhd_opt    (ji,jj-1,jk) ) - bco_bc_hor * zdrho_j(ji,jj-1,jk)  
     894                  zz_dz_j  (ji,jj) = aco_bc_hor * (-gde3w(ji,jj,jk) + gde3w(ji,jj-1,jk) ) - bco_bc_hor * zdz_j  (ji,jj-1,jk) 
     895               END IF 
     896            END IF 
     897         END_2D 
     898         zdrho_i(:,:,jk) = zz_drho_i(:,:) 
     899         zdz_i  (:,:,jk) = zz_dz_i  (:,:) 
     900         zdrho_j(:,:,jk) = zz_drho_j(:,:) 
     901         zdz_j  (:,:,jk) = zz_dz_j  (:,:) 
     902      END DO 
    779903 
    780904      !-------------------------------------------------------------- 
    781       ! Upper half of top-most grid box, compute and store 
     905      ! 7. Calculate integrals on side faces   
    782906      !------------------------------------------------------------- 
    783907 
    784 !!bug gm   :  e3w-gde3w(:,:,:) = 0.5*e3w  ....  and gde3w(:,:,2)-gde3w(:,:,1)=e3w(:,:,2,Kmm) ....   to be verified 
    785 !          true if gde3w(:,:,:) is really defined as the sum of the e3w scale factors as, it seems to me, it should be 
    786  
    787       DO_2D( 0, 0, 0, 0 ) 
    788          rho_k(ji,jj,1) = -grav * ( e3w(ji,jj,1,Kmm) - gde3w(ji,jj,1) )               & 
    789             &                   * (  rhd(ji,jj,1)                                     & 
    790             &                     + 0.5_wp * ( rhd    (ji,jj,2) - rhd    (ji,jj,1) )  & 
    791             &                              * ( e3w  (ji,jj,1,Kmm) - gde3w(ji,jj,1) )  & 
    792             &                              / ( gde3w(ji,jj,2) - gde3w(ji,jj,1) )  ) 
    793       END_2D 
    794  
    795 !!bug gm    : here also, simplification is possible 
    796 !!bug gm    : optimisation: 1/10 and 1/12 the division should be done before the loop 
    797  
    798       DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 
    799  
    800          rho_k(ji,jj,jk) = zcoef0 * ( rhd    (ji,jj,jk) + rhd    (ji,jj,jk-1) )                                   & 
    801             &                     * ( gde3w(ji,jj,jk) - gde3w(ji,jj,jk-1) )                                   & 
    802             &            - grav * z1_10 * (                                                                   & 
    803             &     ( drhow  (ji,jj,jk) - drhow  (ji,jj,jk-1) )                                                     & 
    804             &   * ( gde3w(ji,jj,jk) - gde3w(ji,jj,jk-1) - z1_12 * ( dzw  (ji,jj,jk) + dzw  (ji,jj,jk-1) ) )   & 
    805             &   - ( dzw    (ji,jj,jk) - dzw    (ji,jj,jk-1) )                                                     & 
    806             &   * ( rhd    (ji,jj,jk) - rhd    (ji,jj,jk-1) - z1_12 * ( drhow(ji,jj,jk) + drhow(ji,jj,jk-1) ) )   & 
    807             &                             ) 
    808  
    809          rho_i(ji,jj,jk) = zcoef0 * ( rhd    (ji+1,jj,jk) + rhd    (ji,jj,jk) )                                   & 
    810             &                     * ( gde3w(ji+1,jj,jk) - gde3w(ji,jj,jk) )                                   & 
    811             &            - grav* z1_10 * (                                                                    & 
    812             &     ( drhou  (ji+1,jj,jk) - drhou  (ji,jj,jk) )                                                     & 
    813             &   * ( gde3w(ji+1,jj,jk) - gde3w(ji,jj,jk) - z1_12 * ( dzu  (ji+1,jj,jk) + dzu  (ji,jj,jk) ) )   & 
    814             &   - ( dzu    (ji+1,jj,jk) - dzu    (ji,jj,jk) )                                                     & 
    815             &   * ( rhd    (ji+1,jj,jk) - rhd    (ji,jj,jk) - z1_12 * ( drhou(ji+1,jj,jk) + drhou(ji,jj,jk) ) )   & 
    816             &                            ) 
    817  
    818          rho_j(ji,jj,jk) = zcoef0 * ( rhd    (ji,jj+1,jk) + rhd    (ji,jj,jk) )                                 & 
    819             &                     * ( gde3w(ji,jj+1,jk) - gde3w(ji,jj,jk) )                                   & 
    820             &            - grav* z1_10 * (                                                                    & 
    821             &     ( drhov  (ji,jj+1,jk) - drhov  (ji,jj,jk) )                                                     & 
    822             &   * ( gde3w(ji,jj+1,jk) - gde3w(ji,jj,jk) - z1_12 * ( dzv  (ji,jj+1,jk) + dzv  (ji,jj,jk) ) )   & 
    823             &   - ( dzv    (ji,jj+1,jk) - dzv    (ji,jj,jk) )                                                     & 
    824             &   * ( rhd    (ji,jj+1,jk) - rhd    (ji,jj,jk) - z1_12 * ( drhov(ji,jj+1,jk) + drhov(ji,jj,jk) ) )   & 
    825             &                            ) 
    826  
    827       END_3D 
    828       CALL lbc_lnk_multi( 'dynhpg', rho_k, 'W', 1.0_wp, rho_i, 'U', 1.0_wp, rho_j, 'V', 1.0_wp ) 
     908      DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 
     909! two -ve signs cancel in next two lines (within zcoef0 and because gde3w is a depth not a height) 
     910         z_rho_i(ji,jj,jk) = zcoef0 * ( rhd_opt    (ji+1,jj,jk) + rhd_opt    (ji,jj,jk) )                                       & 
     911             &                    * ( gde3w(ji+1,jj,jk) - gde3w(ji,jj,jk) )                                     
     912         IF ( umask(ji-1, jj, jk) > 0.5 .OR. umask(ji+1, jj, jk) > 0.5 ) THEN 
     913            z_rho_i(ji,jj,jk) = z_rho_i(ji,jj,jk) - z_grav_10 * (                                                               & 
     914             &     (   zdrho_i  (ji+1,jj,jk) - zdrho_i  (ji,jj,jk) )                                                            & 
     915             &   * ( - gde3w(ji+1,jj,jk) + gde3w(ji,jj,jk) - z1_12 * ( zdz_i  (ji+1,jj,jk) + zdz_i  (ji,jj,jk) ) )              & 
     916             &   - (   zdz_i    (ji+1,jj,jk) - zdz_i    (ji,jj,jk) )                                                            & 
     917             &   * (   rhd_opt    (ji+1,jj,jk) - rhd_opt    (ji,jj,jk) - z1_12 * ( zdrho_i(ji+1,jj,jk) + zdrho_i(ji,jj,jk) ) )  & 
     918             &                                               ) 
     919         END IF 
     920   
     921         z_rho_j(ji,jj,jk) = zcoef0 * ( rhd_opt    (ji,jj+1,jk) + rhd_opt    (ji,jj,jk) )                                       & 
     922             &                    * ( gde3w(ji,jj+1,jk) - gde3w(ji,jj,jk) )                                   
     923         IF ( vmask(ji, jj-1, jk) > 0.5 .OR. vmask(ji, jj+1, jk) > 0.5 ) THEN 
     924            z_rho_j(ji,jj,jk) = z_rho_j(ji,jj,jk) - z_grav_10 * (                                                               & 
     925             &     (   zdrho_j  (ji,jj+1,jk) - zdrho_j  (ji,jj,jk) )                                                            & 
     926             &   * ( - gde3w(ji,jj+1,jk) + gde3w(ji,jj,jk) - z1_12 * ( zdz_j  (ji,jj+1,jk) + zdz_j  (ji,jj,jk) ) )              & 
     927             &   - (   zdz_j    (ji,jj+1,jk) - zdz_j    (ji,jj,jk) )                                                            & 
     928             &   * (   rhd_opt    (ji,jj+1,jk) - rhd_opt    (ji,jj,jk) - z1_12 * ( zdrho_j(ji,jj+1,jk) + zdrho_j(ji,jj,jk) ) )  & 
     929             &                                                 ) 
     930         END IF 
     931      END_3D 
     932 
     933      !-------------------------------------------------------------- 
     934      ! 8. Integrate in the vertical    
     935      !------------------------------------------------------------- 
    829936 
    830937      ! --------------- 
     
    832939      ! --------------- 
    833940      DO_2D( 0, 0, 0, 0 ) 
    834          zhpi(ji,jj,1) = ( rho_k(ji+1,jj  ,1) - rho_k(ji,jj,1) - rho_i(ji,jj,1) ) * r1_e1u(ji,jj) 
    835          zhpj(ji,jj,1) = ( rho_k(ji  ,jj+1,1) - rho_k(ji,jj,1) - rho_j(ji,jj,1) ) * r1_e2v(ji,jj) 
     941         zhpi(ji,jj,1) = ( z_rho_k(ji,jj,1) - z_rho_k(ji+1,jj  ,1) - z_rho_i(ji,jj,1) ) * r1_e1u(ji,jj) 
     942         zhpj(ji,jj,1) = ( z_rho_k(ji,jj,1) - z_rho_k(ji  ,jj+1,1) - z_rho_j(ji,jj,1) ) * r1_e2v(ji,jj) 
    836943         IF( ln_wd_il ) THEN 
    837944           zhpi(ji,jj,1) = zhpi(ji,jj,1) * zcpx(ji,jj) 
     
    848955      DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 
    849956         ! hydrostatic pressure gradient along s-surfaces 
    850          zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1)                                & 
    851             &           + (  ( rho_k(ji+1,jj,jk) - rho_k(ji,jj,jk  ) )    & 
    852             &              - ( rho_i(ji  ,jj,jk) - rho_i(ji,jj,jk-1) )  ) * r1_e1u(ji,jj) 
    853          zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1)                                & 
    854             &           + (  ( rho_k(ji,jj+1,jk) - rho_k(ji,jj,jk  ) )    & 
    855             &               -( rho_j(ji,jj  ,jk) - rho_j(ji,jj,jk-1) )  ) * r1_e2v(ji,jj) 
     957         zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1)                                                     & 
     958            &           + (  ( z_rho_k(ji,jj,jk) - z_rho_k(ji+1,jj,jk  ) )                     & 
     959            &              - ( z_rho_i(ji,jj,jk) - z_rho_i(ji  ,jj,jk-1) )  ) * r1_e1u(ji,jj) 
     960         zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1)                                                     & 
     961            &           + (  ( z_rho_k(ji,jj,jk) - z_rho_k(ji,jj+1,jk  ) )                     & 
     962            &               -( z_rho_j(ji,jj,jk) - z_rho_j(ji,jj  ,jk-1) )  ) * r1_e2v(ji,jj) 
    856963         IF( ln_wd_il ) THEN 
    857964           zhpi(ji,jj,jk) = zhpi(ji,jj,jk) * zcpx(ji,jj) 
Note: See TracChangeset for help on using the changeset viewer.