Ignore:
Timestamp:
2016-11-04T06:54:44+01:00 (4 years ago)
Author:
gm
Message:

#1692 - branch SIMPLIF_2_usrdef: e3.=dk[dep.] (discret derivative)

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/NEMO/OPA_SRC/USR/usrdef_zgr.F90

    r7164 r7188  
    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 
     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 
    1717   !!--------------------------------------------------------------------- 
    18    USE oce               ! ocean variables 
    19    USE dom_oce           ! ocean domain 
     18   USE oce            ! ocean variables 
     19   USE dom_oce        ! ocean domain 
     20   USE depth_e3       ! depth <=> e3 
    2021   ! 
    21    USE in_out_manager    ! I/O manager 
    22    USE lbclnk            ! ocean lateral boundary conditions (or mpp link) 
    23    USE lib_mpp           ! distributed memory computing library 
    24    USE wrk_nemo          ! Memory allocation 
    25    USE timing            ! Timing 
     22   USE in_out_manager ! I/O manager 
     23   USE lbclnk         ! ocean lateral boundary conditions (or mpp link) 
     24   USE lib_mpp        ! distributed memory computing library 
     25   USE wrk_nemo       ! Memory allocation 
     26   USE timing         ! Timing 
    2627 
    2728   IMPLICIT NONE 
     
    7273      ! type of vertical coordinate 
    7374      ! --------------------------- 
    74       ld_zco    = .TRUE.         ! GYRE case:  z-coordinate & no ocean cavities 
     75      ld_zco    = .TRUE.         ! GYRE case:  z-coordinate without ocean cavities 
    7576      ld_zps    = .FALSE. 
    7677      ld_sco    = .FALSE. 
     
    8081      ! Build the vertical coordinate system 
    8182      ! ------------------------------------ 
    82       CALL zgr_z  ( pdept_1d, pdepw_1d, pe3t_1d , pe3w_1d )  ! Reference z-coordinate system 
    83       ! 
    84       CALL zgr_msk_top_bot( k_top , k_bot )                  ! masked top and bottom ocean t-level indices 
    85       ! 
    86       !                                                      ! z-coordinate (3D arrays) from the 1D z-coord. 
     83      CALL zgr_z( pdept_1d, pdepw_1d, pe3t_1d , pe3w_1d )   ! Reference z-coordinate system 
     84      ! 
     85      CALL zgr_msk_top_bot( k_top , k_bot )                 ! masked top and bottom ocean t-level indices 
     86      ! 
     87      !                                                     ! z-coordinate (3D arrays) from the 1D z-coord. 
    8788      CALL zgr_zco( pdept_1d, pdepw_1d, pe3t_1d, pe3w_1d,   &   ! in  : 1D reference vertical coordinate 
    8889         &          pdept   , pdepw   ,                     &   ! out : 3D t & w-points depth 
     
    100101      !!              vertical scale factors. 
    101102      !! 
    102       !! ** Method  :   z-coordinate system (use in all type of coordinate) 
    103       !!              The depth of model levels is defined from an analytical 
    104       !!              function the derivative of which gives the scale factors. 
    105       !!              both depth and scale factors only depend on k (1d arrays). 
    106       !!                 w-level: pdepw_1d  = pdep(k) 
    107       !!                          pe3w_1d(k) = dk(pdep)(k)     = e3(k) 
    108       !!                 t-level: pdept_1d  = pdep(k+0.5) 
    109       !!                          pe3t_1d(k) = dk(pdep)(k+0.5) = e3(k+0.5) 
    110       !! 
    111       !!      Here the Madec & Imbard (1996) function is used 
     103      !! ** Method  :   1D z-coordinate system (use in all type of coordinate) 
     104      !!       The depth of model levels is set from dep(k), an analytical function: 
     105      !!                   w-level: depw_1d  = dep(k) 
     106      !!                   t-level: dept_1d  = dep(k+0.5) 
     107      !!       The scale factors are the discrete derivative of the depth: 
     108      !!                   e3w_1d(jk) = dk[ dept_1d ]  
     109      !!                   e3t_1d(jk) = dk[ depw_1d ] 
     110      !!           with at top and bottom : 
     111      !!                   e3w_1d( 1 ) = 2 * ( dept_1d( 1 ) - depw_1d( 1 ) ) 
     112      !!                   e3t_1d(jpk) = 2 * ( dept_1d(jpk) - depw_1d(jpk) ) 
     113      !!       The depth are then re-computed from the sum of e3. This ensures  
     114      !!    that depths are identical when reading domain_cfg.nc file. Indeed, 
     115      !!    Only e3. are saved in this file, depth are compute by a call to  
     116      !!    the e3_to_depth subroutine. 
     117      !! 
     118      !!       Here the Madec & Imbard (1996) function is used. 
    112119      !! 
    113120      !! ** Action  : - pdept_1d, pdepw_1d : depth of T- and W-point (m) 
     
    156163         pdepw_1d(jk) = ( zsur + za0 * zw + za1 * zacr *  LOG( COSH( (zw-zkth) / zacr ) )  ) 
    157164         pdept_1d(jk) = ( zsur + za0 * zt + za1 * zacr *  LOG( COSH( (zt-zkth) / zacr ) )  ) 
    158          pe3w_1d (jk) =          za0      + za1        * TANH(       (zw-zkth) / zacr   ) 
    159          pe3t_1d (jk) =          za0      + za1        * TANH(       (zt-zkth) / zacr   ) 
    160165      END DO 
    161       pdepw_1d(1) = 0._wp                    ! force first w-level to be exactly at zero 
    162  
    163  
    164 !!gm   This should become the reference ! 
    165 !      IF ( ln_isfcav ) THEN 
    166 ! need to be like this to compute the pressure gradient with ISF. If not, level beneath the ISF are not aligned (sum(e3t) /= depth) 
    167 ! define pe3t_0 and pe3w_0 as the differences between pdept and pdepw respectively 
    168 !         DO jk = 1, jpkm1 
    169 !            pe3t_1d(jk) = pdepw_1d(jk+1)-pdepw_1d(jk)  
    170 !         END DO 
    171 !         pe3t_1d(jpk) = pe3t_1d(jpk-1)   ! we don't care because this level is masked in NEMO 
    172 ! 
    173 !         DO jk = 2, jpk 
    174 !            pe3w_1d(jk) = pdept_1d(jk) - pdept_1d(jk-1)  
    175 !         END DO 
    176 !         pe3w_1d(1  ) = 2._wp * (pdept_1d(1) - pdepw_1d(1))  
    177 !      END IF 
    178 !!gm end 
    179  
     166      ! 
     167      !                       ! e3t and e3w from depth 
     168      CALL depth_to_e3( pdept_1d, pdepw_1d, pe3t_1d, pe3w_1d )  
     169      ! 
     170      !                       ! recompute depths from SUM(e3)  <== needed 
     171      CALL e3_to_depth( pe3t_1d, pe3w_1d, pdept_1d, pdepw_1d )  
     172      ! 
    180173      IF(lwp) THEN                        ! control print 
    181174         WRITE(numout,*) 
     
    184177         WRITE(numout, "(10x, i4, 4f9.2)" ) ( jk, pdept_1d(jk), pdepw_1d(jk), pe3t_1d(jk), pe3w_1d(jk), jk = 1, jpk ) 
    185178      ENDIF 
    186       DO jk = 1, jpk                      ! control positivity 
    187          IF( pe3w_1d (jk) <= 0._wp .OR. pe3t_1d (jk) <= 0._wp )   CALL ctl_stop( 'dom:zgr_z: e3w_1d or e3t_1d =< 0 '    ) 
    188          IF( pdepw_1d(jk) <  0._wp .OR. pdept_1d(jk) <  0._wp )   CALL ctl_stop( 'dom:zgr_z: gdepw_1d or gdept_1d < 0 ' ) 
    189       END DO 
    190179      ! 
    191180      IF( nn_timing == 1 )  CALL timing_stop('zgr_z') 
Note: See TracChangeset for help on using the changeset viewer.