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 6440 for branches/UKMO/dev_r5518_GC3p0_package/NEMOGCM/TOOLS/SIREN/src/vgrid.f90 – NEMO

Ignore:
Timestamp:
2016-04-07T16:32:24+02:00 (8 years ago)
Author:
dancopsey
Message:

Merged in nemo_v3_6_STABLE_copy up to revision 6436.

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/UKMO/dev_r5518_GC3p0_package/NEMOGCM/TOOLS/SIREN/src/vgrid.f90

    r5037 r6440  
    7070!> @date Spetember, 2014 
    7171!> - add header 
     72!> @date June, 2015 - update subroutine with NEMO 3.6 
    7273!> 
    7374!> @note Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
     
    118119   !> 
    119120   !> @author G. Madec 
    120    !> - 03,08-  G. Madec: F90: Free form and module 
     121   !> @date Marsh,2008 - F90: Free form and module 
    121122   ! 
    122123   !> @note Reference : Marti, Madec & Delecluse, 1992, JGR, 97, No8, 12,763-12,766. 
     
    139140   !------------------------------------------------------------------- 
    140141   SUBROUTINE vgrid_zgr_z( dd_gdepw, dd_gdept, dd_e3w, dd_e3t,          & 
     142   &                       dd_e3w_1d, dd_e3t_1d, & 
    141143   &                       dd_ppkth, dd_ppkth2, dd_ppacr, dd_ppacr2,    & 
    142144   &                       dd_ppdzmin, dd_pphmax, dd_pp_to_be_computed, & 
     
    148150      REAL(dp), DIMENSION(:), INTENT(INOUT) :: dd_e3w 
    149151      REAL(dp), DIMENSION(:), INTENT(INOUT) :: dd_e3t 
     152      REAL(dp), DIMENSION(:), INTENT(INOUT) :: dd_e3w_1d 
     153      REAL(dp), DIMENSION(:), INTENT(INOUT) :: dd_e3t_1d 
    150154 
    151155      REAL(dp)              , INTENT(IN   ) :: dd_ppkth 
     
    226230         DO jk = 1, il_jpk 
    227231            dl_zw = REAL(jk,dp) 
    228             dl_zt = REAL(jk,dp) + 0.5 
     232            dl_zt = REAL(jk,dp) + 0.5_dp 
    229233            dd_gdepw(jk) = ( dl_zw - 1.0 ) * dl_za1 
    230234            dd_gdept(jk) = ( dl_zt - 1.0 ) * dl_za1 
     
    237241         DO jk = 1, il_jpk 
    238242            dl_zw = REAL( jk,dp) 
    239             dl_zt = REAL( jk,dp) + 0.5 
     243            dl_zt = REAL( jk,dp) + 0.5_dp 
    240244            dd_gdepw(jk) = ( dl_zsur + dl_za0 * dl_zw + & 
    241245            &                dl_za1 * dl_zacr * LOG( COSH( (dl_zw-dl_zkth)/dl_zacr ) ) + & 
     
    255259      ENDIF 
    256260 
     261   ! need to be like this to compute the pressure gradient with ISF. 
     262   ! If not, level beneath the ISF are not aligned (sum(e3t) /= depth) 
     263   ! define e3t_0 and e3w_0 as the differences between gdept and gdepw respectively 
     264      DO jk = 1, il_jpk-1 
     265         dd_e3t_1d(jk) = dd_gdepw(jk+1)-dd_gdepw(jk)  
     266      END DO 
     267      dd_e3t_1d(il_jpk) = dd_e3t_1d(il_jpk-1) ! we don't care because this level is masked in NEMO 
     268 
     269      DO jk = 2, il_jpk 
     270         dd_e3w_1d(jk) = dd_gdept(jk) - dd_gdept(jk-1)  
     271      END DO 
     272      dd_e3w_1d(1  ) = 2._dp * (dd_gdept(1) - dd_gdepw(1)) 
     273 
    257274      ! Control and  print 
    258275      ! ================== 
     
    260277      DO jk = 1, il_jpk 
    261278         IF( dd_e3w(jk)  <= 0. .OR. dd_e3t(jk)  <= 0. )then 
    262             CALL logger_debug("VGRID ZGR Z: e3w or e3t =< 0 ") 
     279            CALL logger_debug("VGRID ZGR Z: e3w or e3t <= 0 ") 
    263280         ENDIF    
     281 
     282         IF( dd_e3w_1d(jk)  <= 0. .OR. dd_e3t_1d(jk)  <= 0. )then 
     283            CALL logger_debug("VGRID ZGR Z: e3w_1d or e3t_1d <= 0 ") 
     284         ENDIF 
    264285 
    265286         IF( dd_gdepw(jk) < 0. .OR. dd_gdept(jk) < 0. )then 
     
    269290 
    270291   END SUBROUTINE vgrid_zgr_z 
     292   !------------------------------------------------------------------- 
     293   !> @brief This subroutine 
     294   !> 
     295   !> @todo add subroutine description 
     296   !> 
     297   !> @param[inout] dd_bathy 
     298   !> @param[in] dd_gdepw 
     299   !> @param[in] dd_hmin 
     300   !> @param[in] dd_fill 
     301   !------------------------------------------------------------------- 
     302   SUBROUTINE vgrid_zgr_bat( dd_bathy, dd_gdepw, dd_hmin, dd_fill ) 
     303      IMPLICIT NONE 
     304      ! Argument 
     305      REAL(dp), DIMENSION(:,:), INTENT(INOUT) :: dd_bathy  
     306      REAL(dp), DIMENSION(:)  , INTENT(IN   ) :: dd_gdepw  
     307      REAL(dp)                , INTENT(IN   ) :: dd_hmin 
     308      REAL(dp)                , INTENT(IN   ), OPTIONAL :: dd_fill 
     309 
     310      ! local 
     311      INTEGER(i4) :: il_jpk 
     312       
     313      REAL(dp)    :: dl_hmin 
     314      REAL(dp)    :: dl_fill 
     315 
     316      ! loop indices 
     317      INTEGER(i4) :: jk 
     318      !---------------------------------------------------------------- 
     319      il_jpk = SIZE(dd_gdepw(:)) 
     320 
     321      dl_fill=0._dp 
     322      IF( PRESENT(dd_fill) ) dl_fill=dd_fill 
     323 
     324      IF( dd_hmin < 0._dp ) THEN 
     325         jk = - INT( dd_hmin )     ! from a nb of level 
     326      ELSE 
     327         jk = MINLOC( dd_gdepw, mask = dd_gdepw > dd_hmin, dim = 1 )  ! from a depth 
     328      ENDIF 
     329       
     330      dl_hmin = dd_gdepw(jk+1) ! minimum depth = ik+1 w-levels  
     331      WHERE( dd_bathy(:,:) <= 0._wp .OR. dd_bathy(:,:) == dl_fill ) 
     332         dd_bathy(:,:) = dl_fill                         ! min=0     over the lands 
     333      ELSE WHERE 
     334         dd_bathy(:,:) = MAX(  dl_hmin , dd_bathy(:,:)  )   ! min=dl_hmin over the oceans 
     335      END WHERE 
     336      WRITE(*,*) 'Minimum ocean depth: ', dl_hmin, ' minimum number of ocean levels : ', jk       
     337 
     338   END SUBROUTINE vgrid_zgr_bat 
    271339   !------------------------------------------------------------------- 
    272340   !> @brief This subroutine set the depth and vertical scale factor in partial step 
     
    311379   !>         - gdept, gdepw and e3 are positives 
    312380   !>         - gdept_ps, gdepw_ps and e3_ps are positives 
    313    ! 
     381   !> 
    314382   !> @author A. Bozec, G. Madec 
    315    !> - 02-09 (A. Bozec, G. Madec) F90: Free form and module 
    316    !> - 02-09 (A. de Miranda)  rigid-lid + islands 
     383   !> @date February, 2009 - F90: Free form and module 
     384   !> @date February, 2009  
     385   !> - A. de Miranda : rigid-lid + islands 
    317386   !> 
    318387   !> @note Reference : Pacanowsky & Gnanadesikan 1997, Mon. Wea. Rev., 126, 3248-3270. 
     
    325394   !> @param[in] dd_e3zps_min 
    326395   !> @param[in] dd_e3zps_rat 
     396   !> @param[in] dd_fill 
    327397   !------------------------------------------------------------------- 
    328398   SUBROUTINE vgrid_zgr_zps( id_mbathy, dd_bathy, id_jpkmax, & 
    329    &                         dd_gdepw, dd_e3t,               & 
    330    &                         dd_e3zps_min, dd_e3zps_rat ) 
     399   &                          dd_gdepw, dd_e3t,               & 
     400   &                          dd_e3zps_min, dd_e3zps_rat,     & 
     401   &                          dd_fill ) 
    331402      IMPLICIT NONE 
    332403      ! Argument       
     
    336407      REAL(dp)   , DIMENSION(:)  , INTENT(IN   ) :: dd_gdepw 
    337408      REAL(dp)   , DIMENSION(:)  , INTENT(IN   ) :: dd_e3t 
    338       REAL(dp)                                   :: dd_e3zps_min 
    339       REAL(dp)                                   :: dd_e3zps_rat 
     409      REAL(dp)                   , INTENT(IN   ) :: dd_e3zps_min 
     410      REAL(dp)                   , INTENT(IN   ) :: dd_e3zps_rat 
     411      REAL(dp)                   , INTENT(IN   ), OPTIONAL :: dd_fill 
    340412 
    341413      ! local variable 
    342414      REAL(dp) :: dl_zmax     ! Maximum depth 
    343       REAL(dp) :: dl_zmin     ! Minimum depth 
     415      !REAL(dp) :: dl_zmin     ! Minimum depth 
    344416      REAL(dp) :: dl_zdepth   ! Ajusted ocean depth to avoid too small e3t  
     417      REAL(dp) :: dl_fill      
    345418 
    346419      INTEGER(i4) :: il_jpk 
     
    359432      il_jpjglo=SIZE(id_mbathy(:,:),DIM=2) 
    360433 
     434      dl_fill=0._dp 
     435      IF( PRESENT(dd_fill) ) dl_fill=dd_fill 
     436 
    361437      ! Initialization of constant 
    362       dl_zmax = dd_gdepw(il_jpk) + dd_e3t(il_jpk) 
    363       dl_zmin = dd_gdepw(4) 
     438      dl_zmax = dd_gdepw(il_jpk) + dd_e3t(il_jpk) ! maximum depth (i.e. the last ocean level thickness <= 2*e3t_1d(jpkm1) ) 
     439 
     440      ! bounded value of bathy (min already set at the end of zgr_bat) 
     441      WHERE( dd_bathy(:,:) /= dl_fill ) 
     442         dd_bathy(:,:) = MIN( dl_zmax ,  dd_bathy(:,:) ) 
     443      END WHERE 
    364444 
    365445      ! bathymetry in level (from bathy_meter) 
     
    372452      DO jj = 1, il_jpjglo 
    373453         DO ji= 1, il_jpiglo 
    374             IF( dd_bathy(ji,jj) <= 0. )   id_mbathy(ji,jj) = INT(dd_bathy(ji,jj),i4) 
    375          END DO 
    376       END DO 
    377  
    378       ! bounded value of bathy 
    379       ! minimum depth == 3 levels 
    380       ! maximum depth == gdepw(jpk)+e3t(jpk)  
    381       ! i.e. the last ocean level thickness cannot exceed e3t(jpkm1)+e3t(jpk) 
    382       DO jj = 1, il_jpjglo 
    383          DO ji= 1, il_jpiglo 
    384             IF( dd_bathy(ji,jj) <= 0. ) THEN 
    385                dd_bathy(ji,jj) = 0.e0 
    386             ELSE 
    387                dd_bathy(ji,jj) = MAX( dd_bathy(ji,jj), dl_zmin ) 
    388                dd_bathy(ji,jj) = MIN( dd_bathy(ji,jj), dl_zmax ) 
     454            IF( dd_bathy(ji,jj) <= 0._dp )THEN 
     455               id_mbathy(ji,jj) = INT(dd_bathy(ji,jj),i4) 
     456            ELSEIF( dd_bathy(ji,jj) == dl_fill )THEN 
     457               id_mbathy(ji,jj) = 0_i4 
    389458            ENDIF 
    390459         END DO 
     
    401470         DO jj = 1, il_jpjglo 
    402471            DO ji = 1, il_jpiglo 
    403                IF( 0. < dd_bathy(ji,jj) .AND. dd_bathy(ji,jj) <= dl_zdepth ) id_mbathy(ji,jj) = jk-1 
     472               IF( dd_bathy(ji,jj) /= dl_fill )THEN 
     473                  IF( 0. < dd_bathy(ji,jj) .AND. & 
     474                  &       dd_bathy(ji,jj) <= dl_zdepth ) id_mbathy(ji,jj) = jk-1 
     475               ENDIF 
    404476            END DO 
    405477         END DO 
     
    432504   !> ** Action  : - update mbathy: level bathymetry (in level index) 
    433505   !>              - update bathy : meter bathymetry (in meters) 
    434  
     506   !> 
    435507   !> @author G.Madec 
    436    !> - 03-08 Original code 
     508   !> @date Marsh, 2008 - Original code 
    437509   ! 
    438510   !> @param[in] id_mbathy  
     
    543615   !> 
    544616   !> @author J.Paul 
    545    !> - November, 2013- Initial Version 
     617   !> @date November, 2013 - Initial Version 
    546618   ! 
    547619   !> @param[in] td_bathy     Bathymetry file structure  
     
    567639      REAL(dp)   , DIMENSION(:)      , ALLOCATABLE :: dl_e3w  
    568640      REAL(dp)   , DIMENSION(:)      , ALLOCATABLE :: dl_e3t 
     641      REAL(dp)   , DIMENSION(:)      , ALLOCATABLE :: dl_e3w_1d  
     642      REAL(dp)   , DIMENSION(:)      , ALLOCATABLE :: dl_e3t_1d 
    569643 
    570644      INTEGER(i4)                                  :: il_status 
     
    710784      ALLOCATE( dl_gdepw(in_nlevel), dl_gdept(in_nlevel) )  
    711785      ALLOCATE(   dl_e3w(in_nlevel),   dl_e3t(in_nlevel) )  
     786      ALLOCATE(   dl_e3w_1d(in_nlevel),   dl_e3t_1d(in_nlevel) )  
    712787      CALL vgrid_zgr_z( dl_gdepw(:), dl_gdept(:), dl_e3w(:), dl_e3t(:), & 
     788      &                 dl_e3w_1d, dl_e3t_1d, & 
    713789      &                 dn_ppkth, dn_ppkth2, dn_ppacr, dn_ppacr2,       & 
    714790      &                 dn_ppdzmin, dn_pphmax, dn_pp_to_be_computed,    & 
Note: See TracChangeset for help on using the changeset viewer.