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 13729 – NEMO

Changeset 13729


Ignore:
Timestamp:
2020-11-05T17:12:35+01:00 (4 years ago)
Author:
ayoung
Message:

Revised hpg_djc scheme. See ticket #2480.

File:
1 edited

Legend:

Unmodified
Added
Removed

  • NEMO/branches/2020/dev_r13723_KERNEL-01_Amy_Mike_newHPGschemes/src/OCE/DYN/dynhpg.F90

    r13295 r13729  
    7373   ! 
    7474   INTEGER, PUBLIC ::   nhpg         !: type of pressure gradient scheme used ! (deduced from ln_hpg_... flags) (PUBLIC for TAM) 
     75   LOGICAL,  PUBLIC ::   ln_hpg_djc_vN_hor, ln_hpg_djc_vN_vrt 
     76   REAL(wp), PUBLIC ::   aco_bc_hor, bco_bc_hor, aco_bc_vrt, bco_bc_vrt 
    7577 
    7678   !! * Substitutions 
    … …  
    156158      !! 
    157159      NAMELIST/namdyn_hpg/ ln_hpg_zco, ln_hpg_zps, ln_hpg_sco,     & 
    158          &                 ln_hpg_djc, ln_hpg_prj, ln_hpg_isf 
     160         &                 ln_hpg_djc, ln_hpg_prj, ln_hpg_isf,     & 
     161         &                 ln_hpg_djc_vN_hor, ln_hpg_djc_vN_vrt 
    159162      !!---------------------------------------------------------------------- 
    160163      ! 
    … …  
    179182      ENDIF 
    180183      ! 
    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) )          & 
     184      IF( .NOT.ln_linssh .AND. .NOT.(ln_hpg_sco.OR.ln_hpg_prj.OR.ln_hpg_isf.OR.ln_hpg_djc) )          & 
    187185         &   CALL ctl_stop('dyn_hpg_init : non-linear free surface requires either ',   & 
    188186         &                 '   the standard jacobian formulation hpg_sco    or '    ,   & 
    … …  
    220218      ENDIF 
    221219      !                           
     220      IF (ln_hpg_djc_vN_hor) THEN ! Von Neumann boundary condition 
     221        aco_bc_hor = 6.0_wp/5.0_wp 
     222        bco_bc_hor = 7.0_wp/15.0_wp 
     223      ELSE ! Linear extrapolation 
     224        aco_bc_hor = 3.0_wp/2.0_wp 
     225        bco_bc_hor = 1.0_wp/2.0_wp 
     226      END IF 
     227      IF (ln_hpg_djc_vN_vrt) THEN ! Von Neumann boundary condition 
     228        aco_bc_vrt = 6.0_wp/5.0_wp 
     229        bco_bc_vrt = 7.0_wp/15.0_wp 
     230      ELSE ! Linear extrapolation 
     231        aco_bc_vrt = 3.0_wp/2.0_wp 
     232        bco_bc_vrt = 1.0_wp/2.0_wp 
     233      END IF 
    222234   END SUBROUTINE dyn_hpg_init 
    223235 
    … …  
    630642      !! 
    631643      INTEGER  ::   ji, jj, jk          ! dummy loop indices 
     644      INTEGER  ::   iktb, iktt          ! jk indices at tracer points for top and bottom points  
    632645      REAL(wp) ::   zcoef0, zep, cffw   ! temporary scalars 
    633       REAL(wp) ::   z1_10, cffu, cffx   !    "         " 
    634       REAL(wp) ::   z1_12, cffv, cffy   !    "         " 
     646      REAL(wp) ::   z_grav_10, z1_12 
     647      REAL(wp) ::   cffu, cffx          !    "         " 
     648      REAL(wp) ::   cffv, cffy          !    "         " 
    635649      LOGICAL  ::   ll_tmp1, ll_tmp2    ! local logical variables 
    636650      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zhpi, zhpj 
    637       REAL(wp), DIMENSION(jpi,jpj,jpk) ::   dzx, dzy, dzz, dzu, dzv, dzw 
    638       REAL(wp), DIMENSION(jpi,jpj,jpk) ::   drhox, drhoy, drhoz, drhou, drhov, drhow 
    639       REAL(wp), DIMENSION(jpi,jpj,jpk) ::   rho_i, rho_j, rho_k 
     651      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   rhd_opt 
     652  
     653      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zdzx, zdzy, zdzz                          ! Primitive grid differences ('delta_xyz') 
     654      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zdz_i, zdz_j, zdz_k                       ! Harmonic average of primitive grid differences ('d_xyz') 
     655      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zdrhox, zdrhoy, zdrhoz                    ! Primitive rho differences ('delta_rho') 
     656      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zdrho_i, zdrho_j, zdrho_k                 ! Harmonic average of primitive rho differences ('d_rho') 
     657      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   z_rho_i, z_rho_j, z_rho_k                 ! Face intergrals 
     658      REAL(wp), DIMENSION(jpi,jpj)     ::   zz_dz_i, zz_dz_j, zz_drho_i, zz_drho_j    ! temporary arrays  
    640659      REAL(wp), DIMENSION(:,:), ALLOCATABLE ::   zcpx, zcpy   !W/D pressure filter 
    641660      !!---------------------------------------------------------------------- 
    … …  
    690709      ! Local constant initialization 
    691710      zcoef0 = - grav * 0.5_wp 
    692       z1_10  = 1._wp / 10._wp 
    693       z1_12  = 1._wp / 12._wp 
     711      z_grav_10  = grav / 10._wp 
     712      z1_12  = 1.0_wp / 12._wp 
     713 
     714 
     715!! 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(:,:,:) 
     716      IF( .NOT. ln_linssh ) THEN 
     717         rhd_opt(:,:,:) = rhd(:,:,:) + 1.0_wp ! for vvl option 
     718      ELSE 
     719         rhd_opt(:,:,:) = rhd(:,:,:) 
     720      END IF 
    694721 
    695722      !---------------------------------------------------------------------------------------- 
    696       !  compute and store in provisional arrays elementary vertical and horizontal differences 
     723      !  1. compute and store elementary vertical differences in provisional arrays  
    697724      !---------------------------------------------------------------------------------------- 
    698725 
    699 !!bug gm   Not a true bug, but... dzz=e3w  for dzx, dzy verify what it is really 
    700  
    701       DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 
    702          drhoz(ji,jj,jk) = rhd    (ji  ,jj  ,jk) - rhd    (ji,jj,jk-1) 
    703          dzz  (ji,jj,jk) = gde3w(ji  ,jj  ,jk) - gde3w(ji,jj,jk-1) 
    704          drhox(ji,jj,jk) = rhd    (ji+1,jj  ,jk) - rhd    (ji,jj,jk  ) 
    705          dzx  (ji,jj,jk) = gde3w(ji+1,jj  ,jk) - gde3w(ji,jj,jk  ) 
    706          drhoy(ji,jj,jk) = rhd    (ji  ,jj+1,jk) - rhd    (ji,jj,jk  ) 
    707          dzy  (ji,jj,jk) = gde3w(ji  ,jj+1,jk) - gde3w(ji,jj,jk  ) 
     726!!bug gm   Not a true bug, but... zdzz=e3w  for zdzx, zdzy verify what it is really 
     727!! zdzz, zdzx and zdzy changed to heights rather than depths; lower bounds of jj and ji changed from 2 to 1  
     728 
     729      DO_3D( 1, 1, 1, 1, 2, jpkm1 )   
     730         zdrhoz(ji,jj,jk) =   rhd_opt    (ji  ,jj  ,jk) - rhd_opt    (ji,jj,jk-1) 
     731         zdzz  (ji,jj,jk) = - gde3w(ji  ,jj  ,jk) + gde3w(ji,jj,jk-1) 
    708732      END_3D 
    709733 
    710734      !------------------------------------------------------------------------- 
    711       ! compute harmonic averages using eq. 5.18 
     735      ! 2. compute harmonic averages for vertical differences using eq. 5.18 
    712736      !------------------------------------------------------------------------- 
    713737      zep = 1.e-15 
    714738 
    715 !!bug  gm  drhoz not defined at level 1 and used (jk-1 with jk=2) 
    716 !!bug  gm  idem for drhox, drhoy et ji=jpi and jj=jpj 
    717  
    718       DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 
    719          cffw = 2._wp * drhoz(ji  ,jj  ,jk) * drhoz(ji,jj,jk-1) 
    720  
    721          cffu = 2._wp * drhox(ji+1,jj  ,jk) * drhox(ji,jj,jk  ) 
    722          cffx = 2._wp * dzx  (ji+1,jj  ,jk) * dzx  (ji,jj,jk  ) 
    723  
    724          cffv = 2._wp * drhoy(ji  ,jj+1,jk) * drhoy(ji,jj,jk  ) 
    725          cffy = 2._wp * dzy  (ji  ,jj+1,jk) * dzy  (ji,jj,jk  ) 
    726  
     739!!bug  gm  zdrhoz not defined at level 1 and used (jk-1 with jk=2) ; this issue has now been addressed 
     740!!bug  gm  idem for zdrhox, zdrhoy et ji=jpi and jj=jpj; idem 
     741 
     742!! mb zdrho_k, zdz_k, zdrho_i, zdz_i, zdrho_j, zdz_j re-centred about the point (ji,jj,jk)  
     743!! mb DO loops broken up so that zdrho_k and zdz_k are calculated only for jk = 2 to jpk - 2  
     744      zdrho_k(:,:,:) = 0._wp 
     745      zdz_k  (:,:,:) = 0._wp 
     746 
     747      DO_3D( 1, 1, 1, 1, 2, jpk-2 )  
     748         cffw = 2._wp * zdrhoz(ji  ,jj  ,jk) * zdrhoz(ji,jj,jk+1) 
    727749         IF( cffw > zep) THEN 
    728             drhow(ji,jj,jk) = 2._wp *   drhoz(ji,jj,jk) * drhoz(ji,jj,jk-1)   & 
    729                &                    / ( drhoz(ji,jj,jk) + drhoz(ji,jj,jk-1) ) 
     750            zdrho_k(ji,jj,jk) = cffw / ( zdrhoz(ji,jj,jk) + zdrhoz(ji,jj,jk+1) ) 
     751         ENDIF 
     752         zdz_k(ji,jj,jk) = 2._wp *   zdzz(ji,jj,jk) * zdzz(ji,jj,jk+1)   & 
     753            &                  / ( zdzz(ji,jj,jk) + zdzz(ji,jj,jk+1) ) 
     754      END_3D 
     755 
     756      !---------------------------------------------------------------------------------- 
     757      ! 3. apply boundary conditions at top and bottom using 5.36-5.37 
     758      !---------------------------------------------------------------------------------- 
     759 
     760! mb for sea-ice shelves we will need to re-write this upper boundary condition in the same form as the lower boundary condition 
     761      zdrho_k(:,:,1) = aco_bc_vrt * ( rhd_opt    (:,:,2) - rhd_opt    (:,:,1) ) - bco_bc_vrt * zdrho_k(:,:,2) 
     762      zdz_k  (:,:,1) = aco_bc_vrt * (-gde3w(:,:,2) + gde3w(:,:,1) ) - bco_bc_vrt * zdz_k  (:,:,2) 
     763 
     764      DO_2D( 1, 1, 1, 1 ) 
     765         IF ( mbkt(ji,jj)>1 ) THEN 
     766            iktb = mbkt(ji,jj) 
     767            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) 
     768            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)  
     769         END IF 
     770      END_2D 
     771 
     772      !-------------------------------------------------------------- 
     773      ! 4. Compute side face integrals 
     774      !------------------------------------------------------------- 
     775 
     776!! ssh replaces e3w_n ; gde3w is a depth; the formulae involve heights   
     777!! rho_k stores grav * FX / rho_0   
     778 
     779      !-------------------------------------------------------------- 
     780      ! 4. a) Upper half of top-most grid box, compute and store 
     781      !------------------------------------------------------------- 
     782! *** AY note: ssh(ji,jj,Kmm) + gde3w(ji,jj,1) = e3w(ji,jj,1) 
     783      DO_2D( 0, 1, 0, 1) 
     784         z_rho_k(ji,jj,1) =  grav * ( ssh(ji,jj,Kmm) + gde3w(ji,jj,1) )                        &  
     785            &                     * (  rhd_opt(ji,jj,1)                                        & 
     786            &                     + 0.5_wp * (   rhd_opt    (ji,jj,2) - rhd_opt    (ji,jj,1) ) & 
     787            &                              * (   ssh   (ji,jj,Kmm) + gde3w(ji,jj,1) )          & 
     788            &                              / ( - gde3w(ji,jj,2) + gde3w(ji,jj,1) )  ) 
     789      END_2D 
     790 
     791      !-------------------------------------------------------------- 
     792      ! 4. b) Interior faces, compute and store 
     793      !------------------------------------------------------------- 
     794 
     795      DO_3D( 0, 1, 0, 1, 2, jpkm1 ) 
     796         z_rho_k(ji,jj,jk) = zcoef0 * (   rhd_opt    (ji,jj,jk) + rhd_opt    (ji,jj,jk-1) )                                   & 
     797            &                       * ( - gde3w(ji,jj,jk) + gde3w(ji,jj,jk-1) )                                               & 
     798            &                       + z_grav_10 * (                                                                           & 
     799            &     (   zdrho_k  (ji,jj,jk) - zdrho_k  (ji,jj,jk-1) )                                                           & 
     800            &   * ( - gde3w(ji,jj,jk) + gde3w(ji,jj,jk-1) - z1_12 * ( zdz_k  (ji,jj,jk) + zdz_k  (ji,jj,jk-1) ) )             & 
     801            &   - ( zdz_k    (ji,jj,jk) - zdz_k    (ji,jj,jk-1) )                                                             & 
     802            &   * ( rhd_opt    (ji,jj,jk) - rhd_opt    (ji,jj,jk-1) - z1_12 * ( zdrho_k(ji,jj,jk) + zdrho_k(ji,jj,jk-1) ) )   & 
     803            &                             ) 
     804      END_3D 
     805 
     806      !---------------------------------------------------------------------------------------- 
     807      !  5. compute and store elementary horizontal differences in provisional arrays  
     808      !---------------------------------------------------------------------------------------- 
     809 
     810      DO_3D( 1, 0, 1, 0, 1, jpkm1 ) 
     811         zdrhox(ji,jj,jk) =   rhd_opt    (ji+1,jj  ,jk) - rhd_opt    (ji,jj,jk  ) 
     812         zdzx  (ji,jj,jk) = - gde3w(ji+1,jj  ,jk) + gde3w(ji,jj,jk  ) 
     813         zdrhoy(ji,jj,jk) =   rhd_opt    (ji  ,jj+1,jk) - rhd_opt    (ji,jj,jk  ) 
     814         zdzy  (ji,jj,jk) = - gde3w(ji  ,jj+1,jk) + gde3w(ji,jj,jk  ) 
     815      END_3D 
     816 
     817      CALL lbc_lnk_multi( 'dynhpg', zdrhox, 'U', 1., zdzx, 'U', 1., zdrhoy, 'V', 1., zdzy, 'V', 1. )  
     818 
     819      !------------------------------------------------------------------------- 
     820      ! 6. compute harmonic averages using eq. 5.18 
     821      !------------------------------------------------------------------------- 
     822 
     823      DO_3D( 0, 1, 0, 1, 1, jpkm1 ) 
     824         cffu = 2._wp * zdrhox(ji-1,jj  ,jk) * zdrhox(ji,jj,jk  ) 
     825         IF( cffu > zep ) THEN 
     826            zdrho_i(ji,jj,jk) = cffu / ( zdrhox(ji-1,jj,jk) + zdrhox(ji,jj,jk) ) 
    730827         ELSE 
    731             drhow(ji,jj,jk) = 0._wp 
    732          ENDIF 
    733  
    734          dzw(ji,jj,jk) = 2._wp *   dzz(ji,jj,jk) * dzz(ji,jj,jk-1)   & 
    735             &                  / ( dzz(ji,jj,jk) + dzz(ji,jj,jk-1) ) 
    736  
    737          IF( cffu > zep ) THEN 
    738             drhou(ji,jj,jk) = 2._wp *   drhox(ji+1,jj,jk) * drhox(ji,jj,jk)   & 
    739                &                    / ( drhox(ji+1,jj,jk) + drhox(ji,jj,jk) ) 
     828            zdrho_i(ji,jj,jk ) = 0._wp 
     829         ENDIF 
     830 
     831         cffx = 2._wp * zdzx  (ji-1,jj  ,jk) * zdzx  (ji,jj,jk  ) 
     832         IF( cffx > zep ) THEN 
     833            zdz_i(ji,jj,jk) = cffx / ( zdzx(ji-1,jj,jk) + zdzx(ji,jj,jk) ) 
    740834         ELSE 
    741             drhou(ji,jj,jk ) = 0._wp 
    742          ENDIF 
    743  
    744          IF( cffx > zep ) THEN 
    745             dzu(ji,jj,jk) = 2._wp *   dzx(ji+1,jj,jk) * dzx(ji,jj,jk)   & 
    746                &                  / ( dzx(ji+1,jj,jk) + dzx(ji,jj,jk) ) 
     835            zdz_i(ji,jj,jk) = 0._wp 
     836         ENDIF 
     837 
     838         cffv = 2._wp * zdrhoy(ji  ,jj-1,jk) * zdrhoy(ji,jj,jk  ) 
     839         IF( cffv > zep ) THEN 
     840            zdrho_j(ji,jj,jk) = cffv / ( zdrhoy(ji,jj-1,jk) + zdrhoy(ji,jj,jk) ) 
    747841         ELSE 
    748             dzu(ji,jj,jk) = 0._wp 
    749          ENDIF 
    750  
    751          IF( cffv > zep ) THEN 
    752             drhov(ji,jj,jk) = 2._wp *   drhoy(ji,jj+1,jk) * drhoy(ji,jj,jk)   & 
    753                &                    / ( drhoy(ji,jj+1,jk) + drhoy(ji,jj,jk) ) 
     842            zdrho_j(ji,jj,jk) = 0._wp 
     843         ENDIF 
     844 
     845         cffy = 2._wp * zdzy  (ji  ,jj-1,jk) * zdzy  (ji,jj,jk  ) 
     846         IF( cffy > zep ) THEN 
     847            zdz_j(ji,jj,jk) = cffy / ( zdzy(ji,jj-1,jk) + zdzy(ji,jj,jk) ) 
    754848         ELSE 
    755             drhov(ji,jj,jk) = 0._wp 
    756          ENDIF 
    757  
    758          IF( cffy > zep ) THEN 
    759             dzv(ji,jj,jk) = 2._wp *   dzy(ji,jj+1,jk) * dzy(ji,jj,jk)   & 
    760                &                  / ( dzy(ji,jj+1,jk) + dzy(ji,jj,jk) ) 
    761          ELSE 
    762             dzv(ji,jj,jk) = 0._wp 
    763          ENDIF 
    764  
    765       END_3D 
     849            zdz_j(ji,jj,jk) = 0._wp 
     850         ENDIF 
     851      END_3D 
     852       
     853!!! Note that zdzx, zdzy, zdzz, zdrhox, zdrhoy and zdrhoz should NOT be used beyond this point       
    766854 
    767855      !---------------------------------------------------------------------------------- 
    768       ! apply boundary conditions at top and bottom using 5.36-5.37 
     856      ! 6B. apply boundary conditions at side boundaries using 5.36-5.37 
    769857      !---------------------------------------------------------------------------------- 
    770       drhow(:,:, 1 ) = 1.5_wp * ( drhoz(:,:, 2 ) - drhoz(:,:,  1  ) ) - 0.5_wp * drhow(:,:,  2  ) 
    771       drhou(:,:, 1 ) = 1.5_wp * ( drhox(:,:, 2 ) - drhox(:,:,  1  ) ) - 0.5_wp * drhou(:,:,  2  ) 
    772       drhov(:,:, 1 ) = 1.5_wp * ( drhoy(:,:, 2 ) - drhoy(:,:,  1  ) ) - 0.5_wp * drhov(:,:,  2  ) 
    773  
    774       drhow(:,:,jpk) = 1.5_wp * ( drhoz(:,:,jpk) - drhoz(:,:,jpkm1) ) - 0.5_wp * drhow(:,:,jpkm1) 
    775       drhou(:,:,jpk) = 1.5_wp * ( drhox(:,:,jpk) - drhox(:,:,jpkm1) ) - 0.5_wp * drhou(:,:,jpkm1) 
    776       drhov(:,:,jpk) = 1.5_wp * ( drhoy(:,:,jpk) - drhoy(:,:,jpkm1) ) - 0.5_wp * drhov(:,:,jpkm1) 
    777  
     858 
     859      DO jk = 1, jpkm1 
     860         zz_drho_i(:,:) = zdrho_i(:,:,jk) 
     861         zz_dz_i  (:,:) = zdz_i  (:,:,jk) 
     862         zz_drho_j(:,:) = zdrho_j(:,:,jk) 
     863         zz_dz_j  (:,:) = zdz_j  (:,:,jk) 
     864         DO_2D( 0, 1, 0, 1) 
     865            ! Walls coming from left: should check from 2 to jpi-1 (and jpj=2-jpj) 
     866            IF (ji < jpi) THEN 
     867               IF ( umask(ji,jj,jk) > 0.5_wp .AND. tmask(ji-1,jj,jk) < 0.5_wp .AND. umask(ji+1,jj,jk) > 0.5_wp)  THEN   
     868                  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)  
     869                  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) 
     870               END IF 
     871            END IF 
     872            ! Walls coming from right: should check from 3 to jpi (and jpj=2-jpj) 
     873            IF (ji > 2) THEN 
     874               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 
     875                  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)   
     876                  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) 
     877               END IF 
     878            END IF 
     879            ! Walls coming from left: should check from 2 to jpj-1 (and jpi=2-jpi) 
     880            IF (jj < jpj) THEN 
     881               IF ( vmask(ji,jj,jk) > 0.5_wp .AND. tmask(ji,jj-1,jk) < 0.5_wp .AND. vmask(ji,jj+1,jk) > 0.5_wp)  THEN 
     882                  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) 
     883                  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) 
     884               END IF 
     885            END IF  
     886            ! Walls coming from right: should check from 3 to jpj (and jpi=2-jpi) 
     887            IF (jj > 2) THEN 
     888               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  
     889                  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)  
     890                  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) 
     891               END IF 
     892            END IF 
     893         END_2D 
     894         zdrho_i(:,:,jk) = zz_drho_i(:,:) 
     895         zdz_i  (:,:,jk) = zz_dz_i  (:,:) 
     896         zdrho_j(:,:,jk) = zz_drho_j(:,:) 
     897         zdz_j  (:,:,jk) = zz_dz_j  (:,:) 
     898      END DO 
    778899 
    779900      !-------------------------------------------------------------- 
    780       ! Upper half of top-most grid box, compute and store 
     901      ! 7. Calculate integrals on side faces   
    781902      !------------------------------------------------------------- 
    782903 
    783 !!bug gm   :  e3w-gde3w(:,:,:) = 0.5*e3w  ....  and gde3w(:,:,2)-gde3w(:,:,1)=e3w(:,:,2,Kmm) ....   to be verified 
    784 !          true if gde3w(:,:,:) is really defined as the sum of the e3w scale factors as, it seems to me, it should be 
    785  
    786       DO_2D( 0, 0, 0, 0 ) 
    787          rho_k(ji,jj,1) = -grav * ( e3w(ji,jj,1,Kmm) - gde3w(ji,jj,1) )               & 
    788             &                   * (  rhd(ji,jj,1)                                     & 
    789             &                     + 0.5_wp * ( rhd    (ji,jj,2) - rhd    (ji,jj,1) )  & 
    790             &                              * ( e3w  (ji,jj,1,Kmm) - gde3w(ji,jj,1) )  & 
    791             &                              / ( gde3w(ji,jj,2) - gde3w(ji,jj,1) )  ) 
    792       END_2D 
    793  
    794 !!bug gm    : here also, simplification is possible 
    795 !!bug gm    : optimisation: 1/10 and 1/12 the division should be done before the loop 
    796  
    797       DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 
    798  
    799          rho_k(ji,jj,jk) = zcoef0 * ( rhd    (ji,jj,jk) + rhd    (ji,jj,jk-1) )                                   & 
    800             &                     * ( gde3w(ji,jj,jk) - gde3w(ji,jj,jk-1) )                                   & 
    801             &            - grav * z1_10 * (                                                                   & 
    802             &     ( drhow  (ji,jj,jk) - drhow  (ji,jj,jk-1) )                                                     & 
    803             &   * ( gde3w(ji,jj,jk) - gde3w(ji,jj,jk-1) - z1_12 * ( dzw  (ji,jj,jk) + dzw  (ji,jj,jk-1) ) )   & 
    804             &   - ( dzw    (ji,jj,jk) - dzw    (ji,jj,jk-1) )                                                     & 
    805             &   * ( rhd    (ji,jj,jk) - rhd    (ji,jj,jk-1) - z1_12 * ( drhow(ji,jj,jk) + drhow(ji,jj,jk-1) ) )   & 
    806             &                             ) 
    807  
    808          rho_i(ji,jj,jk) = zcoef0 * ( rhd    (ji+1,jj,jk) + rhd    (ji,jj,jk) )                                   & 
    809             &                     * ( gde3w(ji+1,jj,jk) - gde3w(ji,jj,jk) )                                   & 
    810             &            - grav* z1_10 * (                                                                    & 
    811             &     ( drhou  (ji+1,jj,jk) - drhou  (ji,jj,jk) )                                                     & 
    812             &   * ( gde3w(ji+1,jj,jk) - gde3w(ji,jj,jk) - z1_12 * ( dzu  (ji+1,jj,jk) + dzu  (ji,jj,jk) ) )   & 
    813             &   - ( dzu    (ji+1,jj,jk) - dzu    (ji,jj,jk) )                                                     & 
    814             &   * ( rhd    (ji+1,jj,jk) - rhd    (ji,jj,jk) - z1_12 * ( drhou(ji+1,jj,jk) + drhou(ji,jj,jk) ) )   & 
    815             &                            ) 
    816  
    817          rho_j(ji,jj,jk) = zcoef0 * ( rhd    (ji,jj+1,jk) + rhd    (ji,jj,jk) )                                 & 
    818             &                     * ( gde3w(ji,jj+1,jk) - gde3w(ji,jj,jk) )                                   & 
    819             &            - grav* z1_10 * (                                                                    & 
    820             &     ( drhov  (ji,jj+1,jk) - drhov  (ji,jj,jk) )                                                     & 
    821             &   * ( gde3w(ji,jj+1,jk) - gde3w(ji,jj,jk) - z1_12 * ( dzv  (ji,jj+1,jk) + dzv  (ji,jj,jk) ) )   & 
    822             &   - ( dzv    (ji,jj+1,jk) - dzv    (ji,jj,jk) )                                                     & 
    823             &   * ( rhd    (ji,jj+1,jk) - rhd    (ji,jj,jk) - z1_12 * ( drhov(ji,jj+1,jk) + drhov(ji,jj,jk) ) )   & 
    824             &                            ) 
    825  
    826       END_3D 
    827       CALL lbc_lnk_multi( 'dynhpg', rho_k, 'W', 1.0_wp, rho_i, 'U', 1.0_wp, rho_j, 'V', 1.0_wp ) 
     904      DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 
     905! two -ve signs cancel in next two lines (within zcoef0 and because gde3w is a depth not a height) 
     906         z_rho_i(ji,jj,jk) = zcoef0 * ( rhd_opt    (ji+1,jj,jk) + rhd_opt    (ji,jj,jk) )                                       & 
     907             &                    * ( gde3w(ji+1,jj,jk) - gde3w(ji,jj,jk) )                                     
     908         IF ( umask(ji-1, jj, jk) > 0.5 .OR. umask(ji+1, jj, jk) > 0.5 ) THEN 
     909            z_rho_i(ji,jj,jk) = z_rho_i(ji,jj,jk) - z_grav_10 * (                                                               & 
     910             &     (   zdrho_i  (ji+1,jj,jk) - zdrho_i  (ji,jj,jk) )                                                            & 
     911             &   * ( - gde3w(ji+1,jj,jk) + gde3w(ji,jj,jk) - z1_12 * ( zdz_i  (ji+1,jj,jk) + zdz_i  (ji,jj,jk) ) )              & 
     912             &   - (   zdz_i    (ji+1,jj,jk) - zdz_i    (ji,jj,jk) )                                                            & 
     913             &   * (   rhd_opt    (ji+1,jj,jk) - rhd_opt    (ji,jj,jk) - z1_12 * ( zdrho_i(ji+1,jj,jk) + zdrho_i(ji,jj,jk) ) )  & 
     914             &                                               ) 
     915         END IF 
     916   
     917         z_rho_j(ji,jj,jk) = zcoef0 * ( rhd_opt    (ji,jj+1,jk) + rhd_opt    (ji,jj,jk) )                                       & 
     918             &                    * ( gde3w(ji,jj+1,jk) - gde3w(ji,jj,jk) )                                   
     919         IF ( vmask(ji, jj-1, jk) > 0.5 .OR. vmask(ji, jj+1, jk) > 0.5 ) THEN 
     920            z_rho_j(ji,jj,jk) = z_rho_j(ji,jj,jk) - z_grav_10 * (                                                               & 
     921             &     (   zdrho_j  (ji,jj+1,jk) - zdrho_j  (ji,jj,jk) )                                                            & 
     922             &   * ( - gde3w(ji,jj+1,jk) + gde3w(ji,jj,jk) - z1_12 * ( zdz_j  (ji,jj+1,jk) + zdz_j  (ji,jj,jk) ) )              & 
     923             &   - (   zdz_j    (ji,jj+1,jk) - zdz_j    (ji,jj,jk) )                                                            & 
     924             &   * (   rhd_opt    (ji,jj+1,jk) - rhd_opt    (ji,jj,jk) - z1_12 * ( zdrho_j(ji,jj+1,jk) + zdrho_j(ji,jj,jk) ) )  & 
     925             &                                                 ) 
     926         END IF 
     927      END_3D 
     928 
     929      !-------------------------------------------------------------- 
     930      ! 8. Integrate in the vertical    
     931      !------------------------------------------------------------- 
    828932 
    829933      ! --------------- 
    … …  
    831935      ! --------------- 
    832936      DO_2D( 0, 0, 0, 0 ) 
    833          zhpi(ji,jj,1) = ( rho_k(ji+1,jj  ,1) - rho_k(ji,jj,1) - rho_i(ji,jj,1) ) * r1_e1u(ji,jj) 
    834          zhpj(ji,jj,1) = ( rho_k(ji  ,jj+1,1) - rho_k(ji,jj,1) - rho_j(ji,jj,1) ) * r1_e2v(ji,jj) 
     937         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) 
     938         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) 
    835939         IF( ln_wd_il ) THEN 
    836940           zhpi(ji,jj,1) = zhpi(ji,jj,1) * zcpx(ji,jj) 
    … …  
    847951      DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 
    848952         ! hydrostatic pressure gradient along s-surfaces 
    849          zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1)                                & 
    850             &           + (  ( rho_k(ji+1,jj,jk) - rho_k(ji,jj,jk  ) )    & 
    851             &              - ( rho_i(ji  ,jj,jk) - rho_i(ji,jj,jk-1) )  ) * r1_e1u(ji,jj) 
    852          zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1)                                & 
    853             &           + (  ( rho_k(ji,jj+1,jk) - rho_k(ji,jj,jk  ) )    & 
    854             &               -( rho_j(ji,jj  ,jk) - rho_j(ji,jj,jk-1) )  ) * r1_e2v(ji,jj) 
     953         zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1)                                                     & 
     954            &           + (  ( z_rho_k(ji,jj,jk) - z_rho_k(ji+1,jj,jk  ) )                     & 
     955            &              - ( z_rho_i(ji,jj,jk) - z_rho_i(ji  ,jj,jk-1) )  ) * r1_e1u(ji,jj) 
     956         zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1)                                                     & 
     957            &           + (  ( z_rho_k(ji,jj,jk) - z_rho_k(ji,jj+1,jk  ) )                     & 
     958            &               -( z_rho_j(ji,jj,jk) - z_rho_j(ji,jj  ,jk-1) )  ) * r1_e2v(ji,jj) 
    855959         IF( ln_wd_il ) THEN 
    856960           zhpi(ji,jj,jk) = zhpi(ji,jj,jk) * zcpx(ji,jj) 
Note: See TracChangeset for help on using the changeset viewer.