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 6904 for branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/CONFIG/OVERFLOW/MY_SRC/usrdef_zgr.F90 – NEMO

Ignore:
Timestamp:
2016-09-01T12:17:31+02:00 (8 years ago)
Author:
gm
Message:

#1692 - branch SIMPLIF_2_usrdef: OVERFLOW configuration (zco & sco) + small bug corrections

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/CONFIG/OVERFLOW/MY_SRC/usrdef_zgr.F90

    r6900 r6904  
    1111 
    1212   !!---------------------------------------------------------------------- 
    13    !!   usr_def_zgr      : user defined vertical coordinate system 
    14    !!       zgr_z        : reference 1D z-coordinate  
    15    !!       zgr_top_bot  : ocean top and bottom level indices 
    16    !!       zgr_zco      : 3D verticl coordinate in pure z-coordinate case 
     13   !!   usr_def_zgr      : user defined vertical coordinate system (required) 
     14   !!       zgr_z1d      : reference 1D z-coordinate  
     15   !!       zgr_zps      : 3D vertical coordinate in z-partial cell coordinate 
    1716   !!--------------------------------------------------------------------- 
    1817   USE oce               ! ocean variables 
    19    USE dom_oce           ! ocean domain 
     18   USE dom_oce  ,  ONLY: ln_zco, ln_zps, ln_sco   ! ocean space and time domain 
     19   USE dom_oce  ,  ONLY: nimpp, njmpp             ! ocean space and time domain 
     20   USE dom_oce  ,  ONLY: glamt                    ! ocean space and time domain 
     21   USE usrdef_nam        ! User defined : namelist variables 
    2022   ! 
    2123   USE in_out_manager    ! I/O manager 
     
    4244      &                    pdept_1d, pdepw_1d, pe3t_1d , pe3w_1d  ,    &   ! 1D reference vertical coordinate 
    4345      &                    pdept , pdepw ,                             &   ! 3D t & w-points depth 
    44       &                    pe3t  , pe3u  , pe3v   , pe3f ,             &   ! vertical scale factors 
    45       &                    pe3w  , pe3uw , pe3vw         ,             &   !     -      -      - 
     46      &                    pe3t  , pe3u  , pe3v , pe3f ,               &   ! vertical scale factors 
     47      &                    pe3w  , pe3uw , pe3vw,                      &   !     -      -      - 
    4648      &                    k_top  , k_bot    )                             ! top & bottom ocean level 
    4749      !!--------------------------------------------------------------------- 
     
    6062      INTEGER , DIMENSION(:,:)  , INTENT(out) ::   k_top, k_bot                ! first & last ocean level 
    6163      ! 
    62       INTEGER  ::   inum   ! local logical unit 
    63       REAL(WP) ::   z_zco, z_zps, z_sco, z_cav 
    64       REAL(wp), DIMENSION(jpi,jpj) ::   z2d   ! 2D workspace 
     64      INTEGER  ::   ji, jj, jk        ! dummy indices 
     65      REAL(wp) ::   zfact, z1_jpkm1   ! local scalar 
     66      REAL(wp) ::   ze3min            ! local scalar 
     67      REAL(wp), DIMENSION(jpi,jpj) ::   zht, zhu, z2d   ! 2D workspace 
     68 
     69 
     70      INTEGER  ::   ik, it, ikb, ikt ! temporary integers 
     71      REAL(wp) ::   ze3tp , ze3wp    ! Last ocean level thickness at T- and W-points 
     72      REAL(wp) ::   zdepwp           ! Ajusted ocean depth to avoid too small e3t 
     73      REAL(wp) ::   zdiff            ! temporary scalar 
     74      REAL(wp) ::   zmax             ! temporary scalar 
     75 
     76 
    6577      !!---------------------------------------------------------------------- 
    6678      ! 
     
    7284      ! type of vertical coordinate 
    7385      ! --------------------------- 
    74       !SF ld_zco    = .TRUE.         ! OVERFLOW case:  z-coordinate & no ocean cavities 
    75       !SF ld_zps    = .FALSE. 
    76       !SF ld_sco    = .FALSE. 
     86      ! already set in usrdef_nam.F90 by reading the namusr_def namelist except for ISF 
    7787      ld_isfcav = .FALSE. 
    7888      ! 
     
    8090      ! Build the vertical coordinate system 
    8191      ! ------------------------------------ 
    82       CALL zgr_z  ( pdept_1d, pdepw_1d, pe3t_1d , pe3w_1d )  ! Reference z-coordinate system 
    83       ! 
    84       CALL zgr_top_bot( k_top , k_bot )                      ! top and bottom ocean level indices 
    85       ! 
    86       !                                                      ! z-coordinate (3D arrays) from the 1D z-coord. 
    87       CALL zgr_zco( pdept_1d, pdepw_1d, pe3t_1d, pe3w_1d,   &   ! in  : 1D reference vertical coordinate 
    88          &          pdept   , pdepw   ,                     &   ! out : 3D t & w-points depth 
    89          &          pe3t    , pe3u    , pe3v   , pe3f   ,   &   !       vertical scale factors 
    90          &          pe3w    , pe3uw   , pe3vw             )     !           -      -      - 
     92      ! 
     93      !                       !==  UNmasked meter bathymetry  ==! 
     94      ! 
     95      ! western continental shelf (500m deep) and eastern deep ocean (2000m deep) 
     96      ! with a hyperbolic tangent transition centered at 40km 
     97      ! NB: here glamt is in kilometers 
     98      ! 
     99      zht(:,:) = + (  500. + 0.5 * 1500. * ( 1.0 + tanh( (glamt(:,:) - 40.) / 7. ) )  ) 
     100      ! 
     101      ! at u-point: recompute from glamu (see usrdef_hgr.F90) to avoid the masking when using lbc_lnk 
     102      zfact = rn_dx * 1.e-3         ! conversion in km 
     103      DO ji = 1, jpi 
     104         zhu(ji,:) = + (  500. + 0.5 * 1500. * ( 1.0 + tanh( ( zfact * REAL( ji-1 + nimpp-1 , wp ) - 40.) / 7. ) )  ) 
     105      END DO 
     106      ! 
     107      CALL zgr_z1d( pdept_1d, pdepw_1d, pe3t_1d , pe3w_1d )   ! Reference z-coordinate system 
     108      ! 
     109      ! 
     110      !                       !==  top masked level bathymetry  ==!  (all coordinates) 
     111      ! 
     112      ! no ocean cavities : top ocean level is ONE, except over land 
     113      ! the ocean basin surrounded by land (1 grid-point) set through lbc_lnk call as jperio=0  
     114      z2d(:,:) = 1._wp                    ! surface ocean is the 1st level 
     115      CALL lbc_lnk( z2d, 'T', 1. )        ! closed basin  
     116      k_top(:,:) = z2d(:,:) 
     117      ! 
     118      !                               
     119      ! 
     120      IF ( ln_sco ) THEN      !==  s-coordinate  ==!   (terrain-following coordinate) 
     121         ! 
     122         k_bot(:,:) = jpkm1 * k_top(:,:)  !* bottom ocean = jpk-1 (here use k_top as a land mask) 
     123         ! 
     124         !                                !* terrain-following coordinate with e3.(k)=cst) 
     125         z1_jpkm1 = 1._wp / REAL( jpkm1 , wp) 
     126         DO jk = 1, jpk 
     127            pdept(:,:,jk) = zht(:,:) * z1_jpkm1 * ( REAL( jk   , wp ) - 0.5_wp ) 
     128            pdepw(:,:,jk) = zht(:,:) * z1_jpkm1 * ( REAL( jk-1 , wp )          ) 
     129            pe3t (:,:,jk) = zht(:,:) * z1_jpkm1 
     130            pe3u (:,:,jk) = zhu(:,:) * z1_jpkm1 
     131            pe3v (:,:,jk) = zht(:,:) * z1_jpkm1 
     132            pe3f (:,:,jk) = zhu(:,:) * z1_jpkm1 
     133            pe3w (:,:,jk) = zht(:,:) * z1_jpkm1 
     134            pe3uw(:,:,jk) = zhu(:,:) * z1_jpkm1 
     135            pe3vw(:,:,jk) = zht(:,:) * z1_jpkm1 
     136         END DO       
     137      ENDIF 
     138      ! 
     139      ! 
     140      IF ( ln_zco ) THEN      !==  z-coordinate  ==!   (step-like topography) 
     141         ! 
     142         !                                !* bottom ocean compute from the depth of grid-points 
     143         k_bot(:,:) = jpkm1 * k_top(:,:)     ! here use k_top as a land mask 
     144         DO jk = 1, jpkm1 
     145            WHERE( pdept_1d(jk) < zht(:,:) .AND. zht(:,:) <= pdept_1d(jk+1) )   k_bot(:,:) = jk * k_top(:,:) 
     146         END DO 
     147         !                                !* horizontally uniform coordinate (reference z-co everywhere) 
     148         DO jk = 1, jpk 
     149            pdept(:,:,jk) = pdept_1d(jk) 
     150            pdepw(:,:,jk) = pdepw_1d(jk) 
     151            pe3t (:,:,jk) = pe3t_1d (jk) 
     152            pe3u (:,:,jk) = pe3t_1d (jk) 
     153            pe3v (:,:,jk) = pe3t_1d (jk) 
     154            pe3f (:,:,jk) = pe3t_1d (jk) 
     155            pe3w (:,:,jk) = pe3w_1d (jk) 
     156            pe3uw(:,:,jk) = pe3w_1d (jk) 
     157            pe3vw(:,:,jk) = pe3w_1d (jk) 
     158         END DO 
     159      ENDIF 
     160      ! 
     161      ! 
     162      IF ( ln_zps ) THEN      !==  zps-coordinate  ==!   (partial bottom-steps) 
     163       
     164       
     165       
     166         CALL ctl_stop( 'STOP', ' zps coordinate not yet coded' ) 
     167 
     168 
     169 
     170       
     171      !!---------------------------------------------------------------------- 
     172      !!                  ***  ROUTINE zgr_zps  *** 
     173      !!                      
     174      !! ** Purpose :   the depth and vertical scale factor in partial step 
     175      !!              reference z-coordinate case 
     176      !! 
     177      !! ** Method  :   Partial steps : computes the 3D vertical scale factors 
     178      !!      of T-, U-, V-, W-, UW-, VW and F-points that are associated with 
     179      !!      a partial step representation of bottom topography. 
     180      !! 
     181      !!       
     182      !!  Reference :   Pacanowsky & Gnanadesikan 1997, Mon. Wea. Rev., 126, 3248-3270. 
     183      !!---------------------------------------------------------------------- 
     184      !!--------------------------------------------------------------------- 
     185         ! 
     186         ze3min = 0.1_wp * rn_dz 
     187         IF(lwp) WRITE(numout,*) '   minimum thickness of the partial cells = 10 % of e3 = ', ze3min 
     188         ! 
     189         ! 
     190         !                                !* bottom ocean compute from the depth of grid-points 
     191         k_bot(:,:) = jpkm1 
     192         DO jk = jpkm1, 1, -1 
     193            WHERE( zht(:,:) <= pdepw_1d(jk) + ze3min )   k_bot(:,:) = jk-1 
     194         END DO 
     195 
     196         !                                !* vertical coordinate system 
     197         ! 
     198         DO jk = 1, jpk                         ! initialization to the reference z-coordinate 
     199            pdept(:,:,jk) = pdept_1d(jk) 
     200            pdepw(:,:,jk) = pdepw_1d(jk) 
     201            pe3t (:,:,jk) = pe3t_1d (jk) 
     202            pe3u (:,:,jk) = pe3t_1d (jk) 
     203            pe3v (:,:,jk) = pe3t_1d (jk) 
     204            pe3f (:,:,jk) = pe3t_1d (jk) 
     205            pe3w (:,:,jk) = pe3w_1d (jk) 
     206            pe3uw(:,:,jk) = pe3w_1d (jk) 
     207            pe3vw(:,:,jk) = pe3w_1d (jk) 
     208         END DO 
     209         ! 
     210         DO jj = 1, jpj                         ! bottom scale factors and depth at T- and W-points 
     211            DO ji = 1, jpi 
     212               ik = k_bot(ji,jj) 
     213               IF( ik /= jpkm1 ) THEN                 ! last level ==> ref 1d z-coordinate 
     214                  pdepw(ji,jj,ik+1) = MIN( zht(ji,jj) , pdepw_1d(ik+1) ) 
     215                  pe3t (ji,jj,ik  ) = pdepw(ji,jj,ik) - pdepw(ji,jj,ik+1) 
     216                  pe3t (ji,jj,ik+1) = pe3t (ji,jj,ik)  
     217                
     218                  pdept(ji,jj,ik  ) = pdept(ji,jj,ik-1) + pe3t(ji,jj,ik) * 0.5_wp 
     219                  pe3w (ji,jj,ik+1) = pdepw(ji,jj,ik) - pdepw(ji,jj,ik+1) 
     220               ENDIF 
     221            END DO 
     222         END DO          
     223          
     224          
     225         ! 
     226         DO jj = 1, jpj                         ! bottom scale factors and depth at T- and W-points 
     227            DO ji = 1, jpi 
     228               ik = k_bot(ji,jj) 
     229                  ! 
     230                  IF( zht(ji,jj) <= pdepw_1d(ik+1) ) THEN  ;   pdepw(ji,jj,ik+1) = zht(ji,jj) 
     231                  ELSE                                     ;   pdepw(ji,jj,ik+1) = pdepw_1d(ik+1) 
     232                  ENDIF 
     233   !gm Bug?  check the gdepw_1d 
     234                  !       ... on ik 
     235                  pdept(ji,jj,ik) = pdepw_1d(ik) + ( pdepw   (ji,jj,ik+1) - pdepw_1d(ik) )   & 
     236                     &                           * ( pdept_1d(      ik  ) - pdepw_1d(ik) )   & 
     237                     &                           / ( pdepw_1d(      ik+1) - pdepw_1d(ik) ) 
     238                  pe3t (ji,jj,ik) = pe3t_1d (ik) * ( pdepw   (ji,jj,ik+1) - pdepw_1d(ik) )   &  
     239                     &                           / ( pdepw_1d(      ik+1) - pdepw_1d(ik) )  
     240                  pe3w (ji,jj,ik) = 0.5_wp * ( pdepw(ji,jj,ik+1) + pdepw_1d(ik+1) - 2._wp * pdepw_1d(ik) )   & 
     241                     &                     * ( pe3w_1d(ik) / ( pdepw_1d(ik+1) - pdepw_1d(ik) ) ) 
     242                  !       ... on ik+1 
     243                  pe3w (ji,jj,ik+1) = pe3t (ji,jj,ik) 
     244                  pe3t (ji,jj,ik+1) = pe3t (ji,jj,ik) 
     245                  pdept(ji,jj,ik+1) = pdept(ji,jj,ik) + pe3t(ji,jj,ik) 
     246            END DO 
     247         END DO 
     248         ! 
     249      ! 
     250      ! 
     251      !                                      ! bottom scale factors and depth at  U-, V-, UW and VW-points 
     252      ! 
     253      DO jk = 1,jpk                          ! Computed as the minimum of neighbooring scale factors 
     254         DO jj = 1, jpjm1 
     255            DO ji = 1, jpim1 
     256               pe3u (ji,jj,jk) = MIN( pe3t(ji,jj,jk), pe3t(ji+1,jj,jk) ) 
     257               pe3v (ji,jj,jk) = MIN( pe3t(ji,jj,jk), pe3t(ji,jj+1,jk) ) 
     258               pe3uw(ji,jj,jk) = MIN( pe3w(ji,jj,jk), pe3w(ji+1,jj,jk) ) 
     259               pe3vw(ji,jj,jk) = MIN( pe3w(ji,jj,jk), pe3w(ji,jj+1,jk) ) 
     260            END DO 
     261         END DO 
     262      END DO 
     263      !                                      ! lateral boundary conditions 
     264      CALL lbc_lnk( pe3u , 'U', 1._wp )   ;   CALL lbc_lnk( pe3uw, 'U', 1._wp )    
     265      CALL lbc_lnk( pe3v , 'V', 1._wp )   ;   CALL lbc_lnk( pe3vw, 'V', 1._wp ) 
     266      ! 
     267 
     268      DO jk = 1, jpk                         ! set to z-scale factor if zero (i.e. along closed boundaries) 
     269         WHERE( pe3u (:,:,jk) == 0._wp )   pe3u (:,:,jk) = pe3t_1d(jk) 
     270         WHERE( pe3v (:,:,jk) == 0._wp )   pe3v (:,:,jk) = pe3t_1d(jk) 
     271         WHERE( pe3uw(:,:,jk) == 0._wp )   pe3uw(:,:,jk) = pe3w_1d(jk) 
     272         WHERE( pe3vw(:,:,jk) == 0._wp )   pe3vw(:,:,jk) = pe3w_1d(jk) 
     273      END DO 
     274       
     275      ! Scale factor at F-point 
     276      DO jk = 1, jpk                        ! initialisation to z-scale factors 
     277         pe3f(:,:,jk) = pe3t_1d(jk) 
     278      END DO 
     279      DO jk = 1, jpk                        ! Computed as the minimum of neighbooring V-scale factors 
     280         DO jj = 1, jpjm1 
     281            DO ji = 1, fs_jpim1   ! vector opt. 
     282               pe3f(ji,jj,jk) = MIN( pe3v(ji,jj,jk), pe3v(ji+1,jj,jk) ) 
     283            END DO 
     284         END DO 
     285      END DO 
     286      CALL lbc_lnk( pe3f, 'F', 1._wp )       ! Lateral boundary conditions 
     287      ! 
     288      DO jk = 1, jpk                        ! set to z-scale factor if zero (i.e. along closed boundaries) 
     289         WHERE( pe3f(:,:,jk) == 0._wp )   pe3f(:,:,jk) = pe3t_1d(jk) 
     290      END DO 
     291!!gm  bug ? :  must be a do loop with mj0,mj1 
     292      !  
     293       
     294!!gm  DO it differently ! 
     295!      pe3t(:,mj0(1),:) = e3t(:,mj0(2),:)     ! we duplicate factor scales for jj = 1 and jj = 2 
     296!      pe3w(:,mj0(1),:) = e3w(:,mj0(2),:)  
     297!      pe3u(:,mj0(1),:) = e3u(:,mj0(2),:)  
     298!      pe3v(:,mj0(1),:) = e3v(:,mj0(2),:)  
     299!      pe3f(:,mj0(1),:) = e3f(:,mj0(2),:)  
     300 
     301      ! Control of the sign 
     302      IF( MINVAL( pe3t (:,:,:) ) <= 0._wp )   CALL ctl_stop( '    zgr_zps :   e r r o r   e3t_0 <= 0' ) 
     303      IF( MINVAL( pe3w (:,:,:) ) <= 0._wp )   CALL ctl_stop( '    zgr_zps :   e r r o r   e3w_0 <= 0' ) 
     304      IF( MINVAL( pdept(:,:,:) ) <  0._wp )   CALL ctl_stop( '    zgr_zps :   e r r o r   gdept_0 <  0' ) 
     305      IF( MINVAL( pdepw(:,:,:) ) <  0._wp )   CALL ctl_stop( '    zgr_zps :   e r r o r   gdepw_0 <  0' ) 
     306      
     307      ! Compute gde3w_0 (vertical sum of e3w) 
     308!      IF ( ln_isfcav ) THEN ! if cavity 
     309!         WHERE( misfdep == 0 )   misfdep = 1 
     310!         DO jj = 1,jpj 
     311!            DO ji = 1,jpi 
     312!               gde3w_0(ji,jj,1) = 0.5_wp * e3w_0(ji,jj,1) 
     313!               DO jk = 2, misfdep(ji,jj) 
     314!                  gde3w_0(ji,jj,jk) = gde3w_0(ji,jj,jk-1) + e3w_0(ji,jj,jk)  
     315!               END DO 
     316!               IF( misfdep(ji,jj) >= 2 )   gde3w_0(ji,jj,misfdep(ji,jj)) = risfdep(ji,jj) + 0.5_wp * e3w_0(ji,jj,misfdep(ji,jj)) 
     317!               DO jk = misfdep(ji,jj) + 1, jpk 
     318!                  gde3w_0(ji,jj,jk) = gde3w_0(ji,jj,jk-1) + e3w_0(ji,jj,jk)  
     319!               END DO 
     320!            END DO 
     321!         END DO 
     322!      ELSE ! no cavity 
     323!         gde3w_0(:,:,1) = 0.5_wp * e3w_0(:,:,1) 
     324!         DO jk = 2, jpk 
     325!            gde3w_0(:,:,jk) = gde3w_0(:,:,jk-1) + e3w_0(:,:,jk) 
     326!         END DO 
     327!      END IF 
     328      ! 
     329      ! 
     330       
     331       
     332       
     333       
     334      ENDIF 
    91335      ! 
    92336   END SUBROUTINE usr_def_zgr 
    93337 
    94338 
    95    SUBROUTINE zgr_z( pdept_1d, pdepw_1d, pe3t_1d , pe3w_1d )   ! 1D reference vertical coordinate 
    96       !!---------------------------------------------------------------------- 
    97       !!                   ***  ROUTINE zgr_z  *** 
     339   SUBROUTINE zgr_z1d( pdept_1d, pdepw_1d, pe3t_1d , pe3w_1d )   ! 1D reference vertical coordinate 
     340      !!---------------------------------------------------------------------- 
     341      !!                   ***  ROUTINE zgr_z1d  *** 
    98342      !! 
    99343      !! ** Purpose :   set the depth of model levels and the resulting  
     
    109353      !!                       pe3t_1d(k) = dk(pdep)(k+0.5) = e3(k+0.5) 
    110354      !! 
    111       !!      Here the Madec & Imbard (1996) function is used 
     355      !!            ===    Here constant vertical resolution   === 
    112356      !! 
    113357      !! ** Action  : - pdept_1d, pdepw_1d : depth of T- and W-point (m) 
    114358      !!              - pe3t_1d , pe3w_1d  : scale factors at T- and W-levels (m) 
    115       !! 
    116       !! Reference : Marti, Madec & Delecluse, 1992, JGR, 97, No8, 12,763-12,766. 
    117       !!             Madec and Imbard, 1996, Clim. Dyn. 
    118       !!---------------------------------------------------------------------- 
    119       REAL(wp), DIMENSION(:)    , INTENT(out) ::   pdept_1d, pdepw_1d   ! 1D grid-point depth        [m] 
    120       REAL(wp), DIMENSION(:)    , INTENT(out) ::   pe3t_1d , pe3w_1d    ! 1D vertical scale factors  [m] 
     359      !!---------------------------------------------------------------------- 
     360      REAL(wp), DIMENSION(:), INTENT(out) ::   pdept_1d, pdepw_1d   ! 1D grid-point depth        [m] 
     361      REAL(wp), DIMENSION(:), INTENT(out) ::   pe3t_1d , pe3w_1d    ! 1D vertical scale factors  [m] 
    121362      ! 
    122363      INTEGER  ::   jk       ! dummy loop indices 
    123       REAL(wp) ::   zt, zw   ! local scalars 
    124       REAL(wp) ::   zsur, za0, za1, zkth, zacr   ! Values for the Madec & Imbard (1996) function   
    125       !!---------------------------------------------------------------------- 
    126       ! 
    127       IF( nn_timing == 1 )  CALL timing_start('zgr_z') 
    128       ! 
    129       ! Set variables from parameters 
    130       ! ------------------------------ 
    131       zsur = -2033.194295283385_wp        
    132       za0  =   155.8325369664153_wp  
    133       za1  =   146.3615918601890_wp 
    134       zkth =    17.28520372419791_wp    
    135       zacr =     5.0_wp        
     364      REAL(wp) ::   zt, zw   ! local scalar 
     365      !!---------------------------------------------------------------------- 
    136366      ! 
    137367      IF(lwp) THEN                         ! Parameter print 
    138368         WRITE(numout,*) 
    139          WRITE(numout,*) '    zgr_z   : Reference vertical z-coordinates ' 
     369         WRITE(numout,*) '    zgr_z1d : Reference vertical z-coordinates: uniform dz = ', rn_dz 
    140370         WRITE(numout,*) '    ~~~~~~~' 
    141          WRITE(numout,*) '       OVERFLOW case : MI96 function with the following coefficients :' 
    142          WRITE(numout,*) '                 zsur = ', zsur 
    143          WRITE(numout,*) '                 za0  = ', za0 
    144          WRITE(numout,*) '                 za1  = ', za1 
    145          WRITE(numout,*) '                 zkth = ', zkth 
    146          WRITE(numout,*) '                 zacr = ', zacr 
    147       ENDIF 
    148  
    149  
     371      ENDIF 
     372      ! 
    150373      ! Reference z-coordinate (depth - scale factor at T- and W-points)   ! Madec & Imbard 1996 function 
    151374      ! ---------------------- 
     
    153376         zw = REAL( jk , wp ) 
    154377         zt = REAL( jk , wp ) + 0.5_wp 
    155          pdepw_1d(jk) = ( zsur + za0 * zw + za1 * zacr *  LOG( COSH( (zw-zkth) / zacr ) )  ) 
    156          pdept_1d(jk) = ( zsur + za0 * zt + za1 * zacr *  LOG( COSH( (zt-zkth) / zacr ) )  ) 
    157          pe3w_1d (jk) =          za0      + za1        * TANH(       (zw-zkth) / zacr   ) 
    158          pe3t_1d (jk) =          za0      + za1        * TANH(       (zt-zkth) / zacr   ) 
    159       END DO 
    160       pdepw_1d(1) = 0._wp                    ! force first w-level to be exactly at zero 
    161  
    162  
     378         pdepw_1d(jk) =    rn_dz *   REAL( jk-1 , wp ) 
     379         pdept_1d(jk) =    rn_dz * ( REAL( jk-1 , wp ) + 0.5_wp ) 
     380         pe3w_1d (jk) =    rn_dz 
     381         pe3t_1d (jk) =    rn_dz 
     382      END DO 
     383      ! 
    163384      IF(lwp) THEN                        ! control print 
    164385         WRITE(numout,*) 
     
    167388         WRITE(numout, "(10x, i4, 4f9.2)" ) ( jk, pdept_1d(jk), pdepw_1d(jk), pe3t_1d(jk), pe3w_1d(jk), jk = 1, jpk ) 
    168389      ENDIF 
    169       DO jk = 1, jpk                      ! control positivity 
    170          IF( pe3w_1d (jk) <= 0._wp .OR. pe3t_1d (jk) <= 0._wp )   CALL ctl_stop( 'dom:zgr_z: e3w_1d or e3t_1d =< 0 '    ) 
    171          IF( pdepw_1d(jk) <  0._wp .OR. pdept_1d(jk) <  0._wp )   CALL ctl_stop( 'dom:zgr_z: gdepw_1d or gdept_1d < 0 ' ) 
    172       END DO 
    173       ! 
    174       IF( nn_timing == 1 )  CALL timing_stop('zgr_z') 
    175       ! 
    176    END SUBROUTINE zgr_z 
    177  
    178  
    179    SUBROUTINE zgr_top_bot( k_top , k_bot ) 
    180       !!---------------------------------------------------------------------- 
    181       !!                    ***  ROUTINE zgr_top_bot  *** 
    182       !! 
    183       !! ** Purpose :   set the top and bottom ocean levels 
    184       !! 
    185       !! ** Method  :   OVERFLOW case = closed box ocean without ocean cavities 
    186       !!                   k_top = 1     except along north, south, east and west boundaries 
    187       !!                   k_bot = jpk-1 except along north, south, east and west boundaries 
    188       !! 
    189       !! ** Action  : - k_top : first ocean level index 
    190       !!              - k_bot : last  ocean level index 
    191       !!---------------------------------------------------------------------- 
    192       INTEGER , DIMENSION(:,:), INTENT(out) ::   k_top , k_bot   ! first & last ocean level 
    193       ! 
    194       INTEGER  ::   ji, jj               ! dummy loop indices 
    195       INTEGER  ::   ii0, ii1, ij0, ij1   ! local indices 
    196       REAL(wp), DIMENSION(jpi,jpj) ::   z2d   ! 2D local workspace 
    197       !!---------------------------------------------------------------------- 
    198       ! 
    199       IF(lwp) WRITE(numout,*) 
    200       IF(lwp) WRITE(numout,*) '    zgr_top_bot : defines the top and bottom ocean levels.' 
    201       IF(lwp) WRITE(numout,*) '    ~~~~~~~~~~~' 
    202       IF(lwp) WRITE(numout,*) '       OVERFLOW case : closed box ocean without ocean cavities' 
    203       ! 
    204       z2d(:,:) = REAL( jpkm1 , wp )          ! bottom 
    205       ! 
    206       CALL lbc_lnk( z2d, 'T', 1. )           ! applied the lateral boundary condition (here jperio=0 ==>> closed) 
    207       ! 
    208       k_bot(:,:) = INT( z2d(:,:) )           ! =jpkm1 over the ocean point, =0 elsewhere 
    209       ! 
    210       k_top(:,:) = MIN( 1 , k_bot(:,:) )     ! =1     over the ocean point, =0 elsewhere 
    211       ! 
    212    END SUBROUTINE zgr_top_bot 
     390      ! 
     391   END SUBROUTINE zgr_z1d 
    213392    
    214  
    215    SUBROUTINE zgr_zco( pdept_1d, pdepw_1d, pe3t_1d, pe3w_1d,   &   ! in : 1D reference vertical coordinate 
    216       &                pdept   , pdepw   ,                     &   ! out: 3D t & w-points depth 
    217       &                pe3t    , pe3u    , pe3v   , pe3f   ,   &   !      vertical scale factors 
    218       &                pe3w    , pe3uw   , pe3vw             )     !          -      -      - 
    219       !!---------------------------------------------------------------------- 
    220       !!                  ***  ROUTINE zgr_zco  *** 
    221       !! 
    222       !! ** Purpose :   define the reference z-coordinate system 
    223       !! 
    224       !! ** Method  :   set 3D coord. arrays to reference 1D array  
    225       !!---------------------------------------------------------------------- 
    226       REAL(wp), DIMENSION(:)    , INTENT(in   ) ::   pdept_1d, pdepw_1d          ! 1D grid-point depth       [m] 
    227       REAL(wp), DIMENSION(:)    , INTENT(in   ) ::   pe3t_1d , pe3w_1d           ! 1D vertical scale factors [m] 
    228       REAL(wp), DIMENSION(:,:,:), INTENT(  out) ::   pdept, pdepw                ! grid-point depth          [m] 
    229       REAL(wp), DIMENSION(:,:,:), INTENT(  out) ::   pe3t , pe3u , pe3v , pe3f   ! vertical scale factors    [m] 
    230       REAL(wp), DIMENSION(:,:,:), INTENT(  out) ::   pe3w , pe3uw, pe3vw         !    -       -      - 
    231       ! 
    232       INTEGER  ::   jk 
    233       !!---------------------------------------------------------------------- 
    234       ! 
    235       IF( nn_timing == 1 )  CALL timing_start('zgr_zco') 
    236       ! 
    237       DO jk = 1, jpk 
    238          pdept(:,:,jk) = pdept_1d(jk) 
    239          pdepw(:,:,jk) = pdepw_1d(jk) 
    240          pe3t (:,:,jk) = pe3t_1d (jk) 
    241          pe3u (:,:,jk) = pe3t_1d (jk) 
    242          pe3v (:,:,jk) = pe3t_1d (jk) 
    243          pe3f (:,:,jk) = pe3t_1d (jk) 
    244          pe3w (:,:,jk) = pe3w_1d (jk) 
    245          pe3uw(:,:,jk) = pe3w_1d (jk) 
    246          pe3vw(:,:,jk) = pe3w_1d (jk) 
    247       END DO 
    248       ! 
    249       IF( nn_timing == 1 )  CALL timing_stop('zgr_zco') 
    250       ! 
    251    END SUBROUTINE zgr_zco 
    252  
    253393   !!====================================================================== 
    254394END MODULE usrdef_zgr 
Note: See TracChangeset for help on using the changeset viewer.