Changeset 15035


Ignore:
Timestamp:
2021-06-21T13:07:02+02:00 (4 months ago)
Author:
ayoung
Message:

Backporting djc into r4.0-HEAD. #2480.

File:
1 edited

Legend:

Unmodified
Added
Removed
  • NEMO/branches/UKMO/NEMO_4.0-HEAD_DJC_backwards_port/src/OCE/DYN/dynhpg.F90

    r13095 r15035  
    7171   ! 
    7272   INTEGER, PUBLIC ::   nhpg         !: type of pressure gradient scheme used ! (deduced from ln_hpg_... flags) (PUBLIC for TAM) 
     73   LOGICAL,  PUBLIC ::   ln_hpg_djc_vN_hor, ln_hpg_djc_vN_vrt 
     74   REAL(wp), PUBLIC ::   aco_bc_hor, bco_bc_hor, aco_bc_vrt, bco_bc_vrt 
    7375 
    7476   !! * Substitutions 
     
    147149      !! 
    148150      NAMELIST/namdyn_hpg/ ln_hpg_zco, ln_hpg_zps, ln_hpg_sco,     & 
    149          &                 ln_hpg_djc, ln_hpg_prj, ln_hpg_isf 
     151         &                 ln_hpg_djc, ln_hpg_prj, ln_hpg_isf,     & 
     152         &                 ln_hpg_djc_vN_hor, ln_hpg_djc_vN_vrt 
    150153      !!---------------------------------------------------------------------- 
    151154      ! 
     
    172175      ENDIF 
    173176      ! 
    174       IF( ln_hpg_djc )   & 
    175          &   CALL ctl_stop('dyn_hpg_init : Density Jacobian: Cubic polynominal method',   & 
    176          &                 '   currently disabled (bugs under investigation).'        ,   & 
    177          &                 '   Please select either  ln_hpg_sco or  ln_hpg_prj instead' ) 
    178          ! 
    179       IF( .NOT.ln_linssh .AND. .NOT.(ln_hpg_sco.OR.ln_hpg_prj.OR.ln_hpg_isf) )          & 
     177      IF( .NOT.ln_linssh .AND. .NOT.(ln_hpg_sco.OR.ln_hpg_prj.OR.ln_hpg_isf.OR.ln_hpg_djc) )          & 
    180178         &   CALL ctl_stop('dyn_hpg_init : non-linear free surface requires either ',   & 
    181179         &                 '   the standard jacobian formulation hpg_sco    or '    ,   & 
     
    213211      ENDIF 
    214212      !                           
     213      IF (ln_hpg_djc_vN_hor) THEN ! Von Neumann boundary condition 
     214        aco_bc_hor = 6.0_wp/5.0_wp 
     215        bco_bc_hor = 7.0_wp/15.0_wp 
     216      ELSE ! Linear extrapolation 
     217        aco_bc_hor = 3.0_wp/2.0_wp 
     218        bco_bc_hor = 1.0_wp/2.0_wp 
     219      END IF 
     220      IF (ln_hpg_djc_vN_vrt) THEN ! Von Neumann boundary condition 
     221        aco_bc_vrt = 6.0_wp/5.0_wp 
     222        bco_bc_vrt = 7.0_wp/15.0_wp 
     223      ELSE ! Linear extrapolation 
     224        aco_bc_vrt = 3.0_wp/2.0_wp 
     225        bco_bc_vrt = 1.0_wp/2.0_wp 
     226      END IF 
    215227      IF ( .NOT. ln_isfcav ) THEN     !--- no ice shelf load 
    216228         riceload(:,:) = 0._wp 
     
    673685      !! 
    674686      INTEGER  ::   ji, jj, jk          ! dummy loop indices 
     687      INTEGER  ::   iktb, iktt          ! jk indices at tracer points for top and bottom points  
    675688      REAL(wp) ::   zcoef0, zep, cffw   ! temporary scalars 
    676       REAL(wp) ::   z1_10, cffu, cffx   !    "         " 
    677       REAL(wp) ::   z1_12, cffv, cffy   !    "         " 
     689      REAL(wp) ::   z_grav_10, z1_12 
     690      REAL(wp) ::   cffu, cffx          !    "         " 
     691      REAL(wp) ::   cffv, cffy          !    "         " 
    678692      LOGICAL  ::   ll_tmp1, ll_tmp2    ! local logical variables 
    679693      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zhpi, zhpj 
    680       REAL(wp), DIMENSION(jpi,jpj,jpk) ::   dzx, dzy, dzz, dzu, dzv, dzw 
    681       REAL(wp), DIMENSION(jpi,jpj,jpk) ::   drhox, drhoy, drhoz, drhou, drhov, drhow 
    682       REAL(wp), DIMENSION(jpi,jpj,jpk) ::   rho_i, rho_j, rho_k 
     694      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   rhd_opt 
     695  
     696      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zdzx, zdzy, zdzz                          ! Primitive grid differences ('delta_xyz') 
     697      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zdz_i, zdz_j, zdz_k                       ! Harmonic average of primitive grid differences ('d_xyz') 
     698      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zdrhox, zdrhoy, zdrhoz                    ! Primitive rho differences ('delta_rho') 
     699      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zdrho_i, zdrho_j, zdrho_k                 ! Harmonic average of primitive rho differences ('d_rho') 
     700      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   z_rho_i, z_rho_j, z_rho_k                 ! Face intergrals 
     701      REAL(wp), DIMENSION(jpi,jpj)     ::   zz_dz_i, zz_dz_j, zz_drho_i, zz_drho_j    ! temporary arrays  
    683702      REAL(wp), DIMENSION(:,:), ALLOCATABLE ::   zcpx, zcpy   !W/D pressure filter 
    684703      !!---------------------------------------------------------------------- 
     
    735754      ! Local constant initialization 
    736755      zcoef0 = - grav * 0.5_wp 
    737       z1_10  = 1._wp / 10._wp 
    738       z1_12  = 1._wp / 12._wp 
    739  
     756      z_grav_10  = grav / 10._wp 
     757      z1_12  = 1.0_wp / 12._wp 
     758 
     759      IF( .NOT. ln_linssh ) THEN 
     760         rhd_opt(:,:,:) = rhd(:,:,:) + 1.0_wp ! for vvl option 
     761      ELSE 
     762         rhd_opt(:,:,:) = rhd(:,:,:) 
     763      END IF 
    740764      !---------------------------------------------------------------------------------------- 
    741       !  compute and store in provisional arrays elementary vertical and horizontal differences 
     765      !  1. compute and store elementary vertical differences in provisional arrays  
    742766      !---------------------------------------------------------------------------------------- 
    743767 
    744 !!bug gm   Not a true bug, but... dzz=e3w  for dzx, dzy verify what it is really 
     768!!bug gm   Not a true bug, but... zdzz=e3w  for zdzx, zdzy verify what it is really 
     769!! zdzz, zdzx and zdzy changed to heights rather than depths; lower bounds of jj and ji changed from 2 to 1  
    745770 
    746771      DO jk = 2, jpkm1 
     772         DO jj = 1, jpj 
     773            DO ji = 1, jpi    
     774               zdrhoz(ji,jj,jk) =   rhd_opt    (ji  ,jj  ,jk) - rhd_opt    (ji,jj,jk-1) 
     775               zdzz  (ji,jj,jk) = - gde3w_n(ji  ,jj  ,jk) + gde3w_n(ji,jj,jk-1) 
     776            END DO 
     777         END DO 
     778      END DO 
     779 
     780      !------------------------------------------------------------------------- 
     781      ! 2. compute harmonic averages for vertical differences using eq. 5.18 
     782      !------------------------------------------------------------------------- 
     783      zep = 1.e-15 
     784 
     785!!bug  gm  zdrhoz not defined at level 1 and used (jk-1 with jk=2) ; this issue has now been addressed 
     786!!bug  gm  idem for zdrhox, zdrhoy et ji=jpi and jj=jpj; idem 
     787 
     788!! zdrho_k, zdz_k, zdrho_i, zdz_i, zdrho_j, zdz_j re-centred about the point (ji,jj,jk)  
     789!! DO loops broken up so that zdrho_k and zdz_k are calculated only for jk = 2 to jpk - 2  
     790      zdrho_k(:,:,:) = 0._wp 
     791      zdz_k  (:,:,:) = 0._wp 
     792 
     793      DO jk = 2, jpk - 2 
     794         DO jj = 1, jpj 
     795            DO ji = 1, jpi    
     796               cffw = 2._wp * zdrhoz(ji  ,jj  ,jk) * zdrhoz(ji,jj,jk+1) 
     797               IF( cffw > zep) THEN 
     798                  zdrho_k(ji,jj,jk) = cffw / ( zdrhoz(ji,jj,jk) + zdrhoz(ji,jj,jk+1) ) 
     799               ENDIF 
     800 
     801               zdz_k(ji,jj,jk) = 2._wp *   zdzz(ji,jj,jk) * zdzz(ji,jj,jk+1)   & 
     802                  &                  / ( zdzz(ji,jj,jk) + zdzz(ji,jj,jk+1) ) 
     803 
     804            END DO 
     805         END DO 
     806      END DO 
     807 
     808      !---------------------------------------------------------------------------------- 
     809      ! 3. apply boundary conditions at top and bottom using 5.36-5.37 
     810      !---------------------------------------------------------------------------------- 
     811 
     812! for sea-ice shelves we will need to re-write this upper boundary condition in the same form as the lower boundary condition 
     813      zdrho_k(:,:,1) = aco_bc_vrt * ( rhd_opt    (:,:,2) - rhd_opt    (:,:,1) ) - bco_bc_vrt * zdrho_k(:,:,2) 
     814      zdz_k  (:,:,1) = aco_bc_vrt * (-gde3w_n(:,:,2) + gde3w_n(:,:,1) ) - bco_bc_vrt * zdz_k  (:,:,2) 
     815 
     816      DO jj = 1, jpj 
     817         DO ji = 1, jpi 
     818            IF ( mbkt(ji,jj)>1 ) THEN 
     819               iktb = mbkt(ji,jj) 
     820               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) 
     821               zdz_k  (ji,jj,iktb) = aco_bc_vrt * (-gde3w_n(ji,jj,iktb) + gde3w_n(ji,jj,iktb-1) ) - bco_bc_vrt * zdz_k  (ji,jj,iktb-1)  
     822            END IF 
     823         END DO 
     824      END DO 
     825 
     826      !-------------------------------------------------------------- 
     827      ! 4. Compute side face integrals 
     828      !------------------------------------------------------------- 
     829 
     830!! sshn replaces e3w_n ; gde3w_n is a depth; the formulae involve heights   
     831!! rho_k stores grav * FX / rho_0   
     832 
     833      !-------------------------------------------------------------- 
     834      ! 4. a) Upper half of top-most grid box, compute and store 
     835      !------------------------------------------------------------- 
     836      ! Concerns that zdrho_k might be oddly defined (just -1.5rho) for single celled columns are resolved by the fact that z_rho_k is defined explicity for the surface layer 
     837      DO jj = 2, jpj 
     838         DO ji = 2, jpi  
     839            z_rho_k(ji,jj,1) =  grav * ( sshn(ji,jj) + gde3w_n(ji,jj,1) )                  & ! *** AY sshn included in djc but not in sco 
     840               &                     * (  rhd_opt(ji,jj,1)                                     & 
     841               &                     + 0.5_wp * (   rhd_opt    (ji,jj,2) - rhd_opt    (ji,jj,1) )  & 
     842               &                              * (   sshn   (ji,jj  ) + gde3w_n(ji,jj,1) )  & 
     843               &                              / ( - gde3w_n(ji,jj,2) + gde3w_n(ji,jj,1) )  ) 
     844         END DO 
     845      END DO 
     846 
     847      !-------------------------------------------------------------- 
     848      ! 4. b) Interior faces, compute and store 
     849      !------------------------------------------------------------- 
     850 
     851      DO jk = 2, jpkm1 
     852         DO jj = 2, jpj 
     853            DO ji = 2, jpi  
     854               z_rho_k(ji,jj,jk) = zcoef0 * (   rhd_opt    (ji,jj,jk) + rhd_opt    (ji,jj,jk-1) )                                   & 
     855                  &                       * ( - gde3w_n(ji,jj,jk) + gde3w_n(ji,jj,jk-1) )                                   & 
     856                  &                       + z_grav_10 * (                                                                   & 
     857                  &     (   zdrho_k  (ji,jj,jk) - zdrho_k  (ji,jj,jk-1) )                                                   & 
     858                  &   * ( - gde3w_n(ji,jj,jk) + gde3w_n(ji,jj,jk-1) - z1_12 * ( zdz_k  (ji,jj,jk) + zdz_k  (ji,jj,jk-1) ) ) & 
     859                  &   - ( zdz_k    (ji,jj,jk) - zdz_k    (ji,jj,jk-1) )                                                     & 
     860                  &   * ( rhd_opt    (ji,jj,jk) - rhd_opt    (ji,jj,jk-1) - z1_12 * ( zdrho_k(ji,jj,jk) + zdrho_k(ji,jj,jk-1) ) )   & 
     861                  &                             ) 
     862            END DO 
     863         END DO 
     864      END DO 
     865 
     866      !---------------------------------------------------------------------------------------- 
     867      !  5. compute and store elementary horizontal differences in provisional arrays  
     868      !---------------------------------------------------------------------------------------- 
     869 
     870      DO jk = 1, jpkm1 
     871         DO jj = 1, jpjm1 
     872            DO ji = 1, fs_jpim1   ! vector opt. 
     873               zdrhox(ji,jj,jk) =   rhd_opt    (ji+1,jj  ,jk) - rhd_opt    (ji,jj,jk  ) 
     874               zdzx  (ji,jj,jk) = - gde3w_n(ji+1,jj  ,jk) + gde3w_n(ji,jj,jk  ) 
     875               zdrhoy(ji,jj,jk) =   rhd_opt    (ji  ,jj+1,jk) - rhd_opt    (ji,jj,jk  ) 
     876               zdzy  (ji,jj,jk) = - gde3w_n(ji  ,jj+1,jk) + gde3w_n(ji,jj,jk  ) 
     877            END DO 
     878         END DO 
     879      END DO 
     880 
     881      CALL lbc_lnk_multi( 'dynhpg', zdrhox, 'U', 1., zdzx, 'U', 1., zdrhoy, 'V', 1., zdzy, 'V', 1. )  
     882 
     883      !------------------------------------------------------------------------- 
     884      ! 6. compute harmonic averages using eq. 5.18 
     885      !------------------------------------------------------------------------- 
     886 
     887      DO jk = 1, jpkm1 
     888         DO jj = 2, jpj 
     889            DO ji = fs_2, jpi   ! vector opt. 
     890 
     891               cffu = 2._wp * zdrhox(ji-1,jj  ,jk) * zdrhox(ji,jj,jk  ) 
     892               IF( cffu > zep ) THEN 
     893                  zdrho_i(ji,jj,jk) = cffu / ( zdrhox(ji-1,jj,jk) + zdrhox(ji,jj,jk) ) 
     894               ELSE 
     895                  zdrho_i(ji,jj,jk ) = 0._wp 
     896               ENDIF 
     897 
     898               cffx = 2._wp * zdzx  (ji-1,jj  ,jk) * zdzx  (ji,jj,jk  ) 
     899               IF( cffx > zep ) THEN 
     900                  zdz_i(ji,jj,jk) = cffx / ( zdzx(ji-1,jj,jk) + zdzx(ji,jj,jk) ) 
     901               ELSE 
     902                  zdz_i(ji,jj,jk) = 0._wp 
     903               ENDIF 
     904 
     905               cffv = 2._wp * zdrhoy(ji  ,jj-1,jk) * zdrhoy(ji,jj,jk  ) 
     906               IF( cffv > zep ) THEN 
     907                  zdrho_j(ji,jj,jk) = cffv / ( zdrhoy(ji,jj-1,jk) + zdrhoy(ji,jj,jk) ) 
     908               ELSE 
     909                  zdrho_j(ji,jj,jk) = 0._wp 
     910               ENDIF 
     911 
     912               cffy = 2._wp * zdzy  (ji  ,jj-1,jk) * zdzy  (ji,jj,jk  ) 
     913               IF( cffy > zep ) THEN 
     914                  zdz_j(ji,jj,jk) = cffy / ( zdzy(ji,jj-1,jk) + zdzy(ji,jj,jk) ) 
     915               ELSE 
     916                  zdz_j(ji,jj,jk) = 0._wp 
     917               ENDIF 
     918 
     919            END DO 
     920         END DO 
     921      END DO 
     922       
     923!!! Note that zdzx, zdzy, zdzz, zdrhox, zdrhoy and zdrhoz should NOT be used beyond this point       
     924 
     925      !---------------------------------------------------------------------------------- 
     926      ! 6B. apply boundary conditions at side boundaries using 5.36-5.37 
     927      !---------------------------------------------------------------------------------- 
     928 
     929      DO jk = 1, jpkm1 
     930         DO jj = 2, jpj 
     931            DO ji = 2, jpi   ! vector opt. 
     932               zz_drho_i(ji,jj) = zdrho_i(ji,jj,jk) 
     933               zz_dz_i  (ji,jj) = zdz_i  (ji,jj,jk) 
     934               zz_drho_j(ji,jj) = zdrho_j(ji,jj,jk) 
     935               zz_dz_j  (ji,jj) = zdz_j  (ji,jj,jk) 
     936               ! Walls coming from left should check from 2 to jpi-1 (and jpj=2-jpj) 
     937               IF (ji < jpi) THEN 
     938                  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   
     939                     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)  
     940                     zz_dz_i  (ji,jj) = aco_bc_hor * (-gde3w_n(ji+1,jj,jk) + gde3w_n(ji,jj,jk) ) - bco_bc_hor * zdz_i  (ji+1,jj,jk) 
     941                  END IF 
     942               END IF 
     943               ! Walls coming from right should check from 3 to jpi (and jpj=2-jpj) 
     944               IF (ji > 2) THEN 
     945                  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 
     946                     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)   
     947                     zz_dz_i  (ji,jj) = aco_bc_hor * (-gde3w_n(ji,jj,jk) + gde3w_n(ji-1,jj,jk) ) - bco_bc_hor * zdz_i  (ji-1,jj,jk) 
     948                  END IF 
     949               END IF 
     950               ! Walls coming from left should check from 2 to jpj-1 (and jpi=2-jpi) 
     951               IF (jj < jpj) THEN 
     952                  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 
     953                     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) 
     954                     zz_dz_j  (ji,jj) = aco_bc_hor * (-gde3w_n(ji,jj+1,jk) + gde3w_n(ji,jj,jk) ) - bco_bc_hor * zdz_j  (ji,jj+1,jk) 
     955                  END IF 
     956               END IF  
     957               ! Walls coming from right should check from 3 to jpj (and jpi=2-jpi) 
     958               IF (jj > 2) THEN 
     959                  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  
     960                     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)  
     961                     zz_dz_j  (ji,jj) = aco_bc_hor * (-gde3w_n(ji,jj,jk) + gde3w_n(ji,jj-1,jk) ) - bco_bc_hor * zdz_j  (ji,jj-1,jk) 
     962                  END IF 
     963               END IF 
     964 
     965            END DO 
     966         END DO 
     967          
     968         DO jj = 2, jpj 
     969            DO ji = 2, jpi   ! vector opt. 
     970              zdrho_i(ji,jj,jk) = zz_drho_i(ji,jj) 
     971              zdz_i  (ji,jj,jk) = zz_dz_i  (ji,jj) 
     972              zdrho_j(ji,jj,jk) = zz_drho_j(ji,jj) 
     973              zdz_j  (ji,jj,jk) = zz_dz_j  (ji,jj) 
     974            END DO 
     975         END DO 
     976 
     977      END DO 
     978 
     979      !-------------------------------------------------------------- 
     980      ! 7. Calculate integrals on side faces   
     981      !------------------------------------------------------------- 
     982 
     983      DO jk = 1, jpkm1 
    747984         DO jj = 2, jpjm1 
    748985            DO ji = fs_2, fs_jpim1   ! vector opt. 
    749                drhoz(ji,jj,jk) = rhd    (ji  ,jj  ,jk) - rhd    (ji,jj,jk-1) 
    750                dzz  (ji,jj,jk) = gde3w_n(ji  ,jj  ,jk) - gde3w_n(ji,jj,jk-1) 
    751                drhox(ji,jj,jk) = rhd    (ji+1,jj  ,jk) - rhd    (ji,jj,jk  ) 
    752                dzx  (ji,jj,jk) = gde3w_n(ji+1,jj  ,jk) - gde3w_n(ji,jj,jk  ) 
    753                drhoy(ji,jj,jk) = rhd    (ji  ,jj+1,jk) - rhd    (ji,jj,jk  ) 
    754                dzy  (ji,jj,jk) = gde3w_n(ji  ,jj+1,jk) - gde3w_n(ji,jj,jk  ) 
    755             END DO 
    756          END DO 
    757       END DO 
    758  
    759       !------------------------------------------------------------------------- 
    760       ! compute harmonic averages using eq. 5.18 
    761       !------------------------------------------------------------------------- 
    762       zep = 1.e-15 
    763  
    764 !!bug  gm  drhoz not defined at level 1 and used (jk-1 with jk=2) 
    765 !!bug  gm  idem for drhox, drhoy et ji=jpi and jj=jpj 
    766  
    767       DO jk = 2, jpkm1 
    768          DO jj = 2, jpjm1 
    769             DO ji = fs_2, fs_jpim1   ! vector opt. 
    770                cffw = 2._wp * drhoz(ji  ,jj  ,jk) * drhoz(ji,jj,jk-1) 
    771  
    772                cffu = 2._wp * drhox(ji+1,jj  ,jk) * drhox(ji,jj,jk  ) 
    773                cffx = 2._wp * dzx  (ji+1,jj  ,jk) * dzx  (ji,jj,jk  ) 
    774  
    775                cffv = 2._wp * drhoy(ji  ,jj+1,jk) * drhoy(ji,jj,jk  ) 
    776                cffy = 2._wp * dzy  (ji  ,jj+1,jk) * dzy  (ji,jj,jk  ) 
    777  
    778                IF( cffw > zep) THEN 
    779                   drhow(ji,jj,jk) = 2._wp *   drhoz(ji,jj,jk) * drhoz(ji,jj,jk-1)   & 
    780                      &                    / ( drhoz(ji,jj,jk) + drhoz(ji,jj,jk-1) ) 
    781                ELSE 
    782                   drhow(ji,jj,jk) = 0._wp 
    783                ENDIF 
    784  
    785                dzw(ji,jj,jk) = 2._wp *   dzz(ji,jj,jk) * dzz(ji,jj,jk-1)   & 
    786                   &                  / ( dzz(ji,jj,jk) + dzz(ji,jj,jk-1) ) 
    787  
    788                IF( cffu > zep ) THEN 
    789                   drhou(ji,jj,jk) = 2._wp *   drhox(ji+1,jj,jk) * drhox(ji,jj,jk)   & 
    790                      &                    / ( drhox(ji+1,jj,jk) + drhox(ji,jj,jk) ) 
    791                ELSE 
    792                   drhou(ji,jj,jk ) = 0._wp 
    793                ENDIF 
    794  
    795                IF( cffx > zep ) THEN 
    796                   dzu(ji,jj,jk) = 2._wp *   dzx(ji+1,jj,jk) * dzx(ji,jj,jk)   & 
    797                      &                  / ( dzx(ji+1,jj,jk) + dzx(ji,jj,jk) ) 
    798                ELSE 
    799                   dzu(ji,jj,jk) = 0._wp 
    800                ENDIF 
    801  
    802                IF( cffv > zep ) THEN 
    803                   drhov(ji,jj,jk) = 2._wp *   drhoy(ji,jj+1,jk) * drhoy(ji,jj,jk)   & 
    804                      &                    / ( drhoy(ji,jj+1,jk) + drhoy(ji,jj,jk) ) 
    805                ELSE 
    806                   drhov(ji,jj,jk) = 0._wp 
    807                ENDIF 
    808  
    809                IF( cffy > zep ) THEN 
    810                   dzv(ji,jj,jk) = 2._wp *   dzy(ji,jj+1,jk) * dzy(ji,jj,jk)   & 
    811                      &                  / ( dzy(ji,jj+1,jk) + dzy(ji,jj,jk) ) 
    812                ELSE 
    813                   dzv(ji,jj,jk) = 0._wp 
    814                ENDIF 
    815  
    816             END DO 
    817          END DO 
    818       END DO 
    819  
    820       !---------------------------------------------------------------------------------- 
    821       ! apply boundary conditions at top and bottom using 5.36-5.37 
    822       !---------------------------------------------------------------------------------- 
    823       drhow(:,:, 1 ) = 1.5_wp * ( drhoz(:,:, 2 ) - drhoz(:,:,  1  ) ) - 0.5_wp * drhow(:,:,  2  ) 
    824       drhou(:,:, 1 ) = 1.5_wp * ( drhox(:,:, 2 ) - drhox(:,:,  1  ) ) - 0.5_wp * drhou(:,:,  2  ) 
    825       drhov(:,:, 1 ) = 1.5_wp * ( drhoy(:,:, 2 ) - drhoy(:,:,  1  ) ) - 0.5_wp * drhov(:,:,  2  ) 
    826  
    827       drhow(:,:,jpk) = 1.5_wp * ( drhoz(:,:,jpk) - drhoz(:,:,jpkm1) ) - 0.5_wp * drhow(:,:,jpkm1) 
    828       drhou(:,:,jpk) = 1.5_wp * ( drhox(:,:,jpk) - drhox(:,:,jpkm1) ) - 0.5_wp * drhou(:,:,jpkm1) 
    829       drhov(:,:,jpk) = 1.5_wp * ( drhoy(:,:,jpk) - drhoy(:,:,jpkm1) ) - 0.5_wp * drhov(:,:,jpkm1) 
    830  
     986 
     987! two -ve signs cancel in next two lines (within zcoef0 and because gde3w is a depth not a height) 
     988               z_rho_i(ji,jj,jk) = zcoef0 * ( rhd_opt    (ji+1,jj,jk) + rhd_opt    (ji,jj,jk) )                                     & 
     989                   &                    * ( gde3w_n(ji+1,jj,jk) - gde3w_n(ji,jj,jk) )                                     
     990               IF ( umask(ji-1, jj, jk) > 0.5 .OR. umask(ji+1, jj, jk) > 0.5 ) THEN 
     991                  z_rho_i(ji,jj,jk) = z_rho_i(ji,jj,jk) - z_grav_10 * (                                                       & 
     992                   &     (   zdrho_i  (ji+1,jj,jk) - zdrho_i  (ji,jj,jk) )                                                    & 
     993                   &   * ( - gde3w_n(ji+1,jj,jk) + gde3w_n(ji,jj,jk) - z1_12 * ( zdz_i  (ji+1,jj,jk) + zdz_i  (ji,jj,jk) ) )  & 
     994                   &   - (   zdz_i    (ji+1,jj,jk) - zdz_i    (ji,jj,jk) )                                                    & 
     995                   &   * (   rhd_opt    (ji+1,jj,jk) - rhd_opt    (ji,jj,jk) - z1_12 * ( zdrho_i(ji+1,jj,jk) + zdrho_i(ji,jj,jk) ) )  & 
     996                   &                                               ) 
     997               END IF 
     998   
     999               z_rho_j(ji,jj,jk) = zcoef0 * ( rhd_opt    (ji,jj+1,jk) + rhd_opt    (ji,jj,jk) )                                     & 
     1000                   &                    * ( gde3w_n(ji,jj+1,jk) - gde3w_n(ji,jj,jk) )                                   
     1001               IF ( vmask(ji, jj-1, jk) > 0.5 .OR. vmask(ji, jj+1, jk) > 0.5 ) THEN 
     1002                  z_rho_j(ji,jj,jk) = z_rho_j(ji,jj,jk) - z_grav_10 * (                                                       & 
     1003                   &     (   zdrho_j  (ji,jj+1,jk) - zdrho_j  (ji,jj,jk) )                                                    & 
     1004                   &   * ( - gde3w_n(ji,jj+1,jk) + gde3w_n(ji,jj,jk) - z1_12 * ( zdz_j  (ji,jj+1,jk) + zdz_j  (ji,jj,jk) ) )  & 
     1005                   &   - (   zdz_j    (ji,jj+1,jk) - zdz_j    (ji,jj,jk) )                                                    & 
     1006                   &   * (   rhd_opt    (ji,jj+1,jk) - rhd_opt    (ji,jj,jk) - z1_12 * ( zdrho_j(ji,jj+1,jk) + zdrho_j(ji,jj,jk) ) )  & 
     1007                   &                                                 ) 
     1008               END IF 
     1009 
     1010            END DO 
     1011         END DO 
     1012      END DO 
    8311013 
    8321014      !-------------------------------------------------------------- 
    833       ! Upper half of top-most grid box, compute and store 
     1015      ! 8. Integrate in the vertical    
    8341016      !------------------------------------------------------------- 
    835  
    836 !!bug gm   :  e3w-gde3w = 0.5*e3w  ....  and gde3w(2)-gde3w(1)=e3w(2) ....   to be verified 
    837 !          true if gde3w is really defined as the sum of the e3w scale factors as, it seems to me, it should be 
    838  
    839       DO jj = 2, jpjm1 
    840          DO ji = fs_2, fs_jpim1   ! vector opt. 
    841             rho_k(ji,jj,1) = -grav * ( e3w_n(ji,jj,1) - gde3w_n(ji,jj,1) )               & 
    842                &                   * (  rhd(ji,jj,1)                                     & 
    843                &                     + 0.5_wp * ( rhd    (ji,jj,2) - rhd    (ji,jj,1) )  & 
    844                &                              * ( e3w_n  (ji,jj,1) - gde3w_n(ji,jj,1) )  & 
    845                &                              / ( gde3w_n(ji,jj,2) - gde3w_n(ji,jj,1) )  ) 
    846          END DO 
    847       END DO 
    848  
    849 !!bug gm    : here also, simplification is possible 
    850 !!bug gm    : optimisation: 1/10 and 1/12 the division should be done before the loop 
    851  
    852       DO jk = 2, jpkm1 
    853          DO jj = 2, jpjm1 
    854             DO ji = fs_2, fs_jpim1   ! vector opt. 
    855  
    856                rho_k(ji,jj,jk) = zcoef0 * ( rhd    (ji,jj,jk) + rhd    (ji,jj,jk-1) )                                   & 
    857                   &                     * ( gde3w_n(ji,jj,jk) - gde3w_n(ji,jj,jk-1) )                                   & 
    858                   &            - grav * z1_10 * (                                                                   & 
    859                   &     ( drhow  (ji,jj,jk) - drhow  (ji,jj,jk-1) )                                                     & 
    860                   &   * ( gde3w_n(ji,jj,jk) - gde3w_n(ji,jj,jk-1) - z1_12 * ( dzw  (ji,jj,jk) + dzw  (ji,jj,jk-1) ) )   & 
    861                   &   - ( dzw    (ji,jj,jk) - dzw    (ji,jj,jk-1) )                                                     & 
    862                   &   * ( rhd    (ji,jj,jk) - rhd    (ji,jj,jk-1) - z1_12 * ( drhow(ji,jj,jk) + drhow(ji,jj,jk-1) ) )   & 
    863                   &                             ) 
    864  
    865                rho_i(ji,jj,jk) = zcoef0 * ( rhd    (ji+1,jj,jk) + rhd    (ji,jj,jk) )                                   & 
    866                   &                     * ( gde3w_n(ji+1,jj,jk) - gde3w_n(ji,jj,jk) )                                   & 
    867                   &            - grav* z1_10 * (                                                                    & 
    868                   &     ( drhou  (ji+1,jj,jk) - drhou  (ji,jj,jk) )                                                     & 
    869                   &   * ( gde3w_n(ji+1,jj,jk) - gde3w_n(ji,jj,jk) - z1_12 * ( dzu  (ji+1,jj,jk) + dzu  (ji,jj,jk) ) )   & 
    870                   &   - ( dzu    (ji+1,jj,jk) - dzu    (ji,jj,jk) )                                                     & 
    871                   &   * ( rhd    (ji+1,jj,jk) - rhd    (ji,jj,jk) - z1_12 * ( drhou(ji+1,jj,jk) + drhou(ji,jj,jk) ) )   & 
    872                   &                            ) 
    873  
    874                rho_j(ji,jj,jk) = zcoef0 * ( rhd    (ji,jj+1,jk) + rhd    (ji,jj,jk) )                                 & 
    875                   &                     * ( gde3w_n(ji,jj+1,jk) - gde3w_n(ji,jj,jk) )                                   & 
    876                   &            - grav* z1_10 * (                                                                    & 
    877                   &     ( drhov  (ji,jj+1,jk) - drhov  (ji,jj,jk) )                                                     & 
    878                   &   * ( gde3w_n(ji,jj+1,jk) - gde3w_n(ji,jj,jk) - z1_12 * ( dzv  (ji,jj+1,jk) + dzv  (ji,jj,jk) ) )   & 
    879                   &   - ( dzv    (ji,jj+1,jk) - dzv    (ji,jj,jk) )                                                     & 
    880                   &   * ( rhd    (ji,jj+1,jk) - rhd    (ji,jj,jk) - z1_12 * ( drhov(ji,jj+1,jk) + drhov(ji,jj,jk) ) )   & 
    881                   &                            ) 
    882  
    883             END DO 
    884          END DO 
    885       END DO 
    886       CALL lbc_lnk_multi( 'dynhpg', rho_k, 'W', 1., rho_i, 'U', 1., rho_j, 'V', 1. ) 
    8871017 
    8881018      ! --------------- 
     
    8911021      DO jj = 2, jpjm1 
    8921022         DO ji = fs_2, fs_jpim1   ! vector opt. 
    893             zhpi(ji,jj,1) = ( rho_k(ji+1,jj  ,1) - rho_k(ji,jj,1) - rho_i(ji,jj,1) ) * r1_e1u(ji,jj) 
    894             zhpj(ji,jj,1) = ( rho_k(ji  ,jj+1,1) - rho_k(ji,jj,1) - rho_j(ji,jj,1) ) * r1_e2v(ji,jj) 
     1023            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) 
     1024            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) 
    8951025            IF( ln_wd_il ) THEN 
    8961026              zhpi(ji,jj,1) = zhpi(ji,jj,1) * zcpx(ji,jj) 
     
    9111041               ! hydrostatic pressure gradient along s-surfaces 
    9121042               zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1)                                & 
    913                   &           + (  ( rho_k(ji+1,jj,jk) - rho_k(ji,jj,jk  ) )    & 
    914                   &              - ( rho_i(ji  ,jj,jk) - rho_i(ji,jj,jk-1) )  ) * r1_e1u(ji,jj) 
     1043                  &           + (  ( z_rho_k(ji,jj,jk) - z_rho_k(ji+1,jj,jk  ) )    & 
     1044                  &              - ( z_rho_i(ji,jj,jk) - z_rho_i(ji  ,jj,jk-1) )  ) * r1_e1u(ji,jj) 
    9151045               zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1)                                & 
    916                   &           + (  ( rho_k(ji,jj+1,jk) - rho_k(ji,jj,jk  ) )    & 
    917                   &               -( rho_j(ji,jj  ,jk) - rho_j(ji,jj,jk-1) )  ) * r1_e2v(ji,jj) 
     1046                  &           + (  ( z_rho_k(ji,jj,jk) - z_rho_k(ji,jj+1,jk  ) )    & 
     1047                  &               -( z_rho_j(ji,jj,jk) - z_rho_j(ji,jj  ,jk-1) )  ) * r1_e2v(ji,jj) 
    9181048               IF( ln_wd_il ) THEN 
    9191049                 zhpi(ji,jj,jk) = zhpi(ji,jj,jk) * zcpx(ji,jj) 
Note: See TracChangeset for help on using the changeset viewer.