Ignore:
Timestamp:
2018-12-11T21:00:08+01:00 (2 years ago)
Author:
clem
Message:

attempt to correct several bugs in the NESTING tools. With this version, domcfg.nc file should be written correctly and the partial steps should be taken into account.

File:
1 edited

Legend:

Unmodified
Added
Removed
  • utils/tools/NESTING/src/agrif_readwrite.f90

    r10160 r10381  
    209209    TYPE(Coordinates) :: Grid 
    210210    ! 
    211     CALL read_ncdf_var('mbathy',name,Grid%Bathy_level)     
     211    CALL read_ncdf_var(parent_level_name,name,Grid%Bathy_level)     
    212212    ! 
    213213    WRITE(*,*) ' ' 
     
    239239    CALL write_ncdf_dim(dimnames(2),name,nyfin) 
    240240    !       
    241     CALL write_ncdf_var('nav_lon',dimnames,name,Grid%nav_lon    ,'float') 
    242     CALL write_ncdf_var('nav_lat',dimnames,name,Grid%nav_lat    ,'float') 
    243     CALL write_ncdf_var('mbathy' ,dimnames,name,Grid%bathy_level,'float') 
    244     ! 
    245     CALL copy_ncdf_att('nav_lon',TRIM(parent_bathy_file),name,MINVAL(Grid%nav_lon),MAXVAL(Grid%nav_lon)) 
    246     CALL copy_ncdf_att('nav_lat',TRIM(parent_bathy_file),name,MINVAL(Grid%nav_lat),MAXVAL(Grid%nav_lat)) 
    247     CALL copy_ncdf_att('mbathy' ,TRIM(parent_bathy_file),name)        
     241    CALL write_ncdf_var('nav_lon'        ,dimnames,name,Grid%nav_lon    ,'float') 
     242    CALL write_ncdf_var('nav_lat'        ,dimnames,name,Grid%nav_lat    ,'float') 
     243    CALL write_ncdf_var(parent_level_name,dimnames,name,Grid%bathy_level,'float') 
     244    ! 
     245    CALL copy_ncdf_att('nav_lon'        ,TRIM(parent_bathy_level),name,MINVAL(Grid%nav_lon),MAXVAL(Grid%nav_lon)) 
     246    CALL copy_ncdf_att('nav_lat'        ,TRIM(parent_bathy_level),name,MINVAL(Grid%nav_lat),MAXVAL(Grid%nav_lat)) 
     247    CALL copy_ncdf_att(parent_level_name,TRIM(parent_bathy_level),name)        
    248248    ! 
    249249    WRITE(*,*) ' ' 
     
    256256  ! 
    257257  !***************************************************** 
    258   !   function Read_Bathy_meter(name,CoarseGrid,ChildGrid) 
    259   !***************************************************** 
    260   ! 
    261   INTEGER FUNCTION Read_Bathy_meter(name,CoarseGrid,ChildGrid,Pacifique) 
     258  !   function read_bathy_coord(name,CoarseGrid,ChildGrid) 
     259  !***************************************************** 
     260  ! 
     261  INTEGER FUNCTION read_bathy_coord(name,CoarseGrid,ChildGrid,Pacifique) 
    262262    ! 
    263263    USE io_netcdf 
     
    284284       CALL read_ncdf_var('nav_lon',name,CoarseGrid%nav_lon) 
    285285       CALL read_ncdf_var('nav_lat',name,CoarseGrid%nav_lat) 
    286        CALL read_ncdf_var(parent_batmet_name,name,CoarseGrid%Bathy_meter) 
     286       CALL read_ncdf_var(parent_meter_name,name,CoarseGrid%Bathy_meter) 
    287287       !             
    288288       IF ( PRESENT(Pacifique) ) THEN 
     
    294294       ENDIF 
    295295       !       
    296        Read_Bathy_meter = 1 
     296       read_bathy_coord = 1 
    297297       RETURN       
    298298    ELSE 
     
    411411    CoarseGrid%Bathy_meter(:,:) = -1.0 * CoarseGrid%Bathy_meter(:,:) 
    412412    !       
    413     Read_Bathy_meter = 1 
     413    read_bathy_coord = 1 
    414414    RETURN 
    415415    !       
    416   END FUNCTION Read_Bathy_meter 
     416  END FUNCTION read_bathy_coord 
    417417  !  
    418418  ! 
    419419  !***************************************************** 
    420   !   function Read_Bathy_meter(name,CoarseGrid,ChildGrid) 
    421   !***************************************************** 
    422   ! 
    423   INTEGER FUNCTION Read_Bathymeter(name,Grid) 
     420  !   function read_bathy_meter(name,CoarseGrid,ChildGrid) 
     421  !***************************************************** 
     422  ! 
     423  INTEGER FUNCTION read_bathy_meter(name,Grid) 
    424424    ! 
    425425    ! 
     
    429429    TYPE(Coordinates) :: Grid 
    430430    ! 
    431     CALL read_ncdf_var(parent_batmet_name,name,Grid%Bathy_meter)     
     431    CALL read_ncdf_var(parent_meter_name,name,Grid%Bathy_meter)     
    432432    ! 
    433433    WRITE(*,*) ' ' 
     
    435435    WRITE(*,*) ' '       
    436436    ! 
    437     Read_Bathymeter = 1 
    438     !       
    439   END FUNCTION Read_Bathymeter 
     437    read_bathy_meter = 1 
     438    !       
     439  END FUNCTION read_bathy_meter 
    440440  !      
    441441  !***************************************************** 
     
    463463    CALL write_ncdf_dim(dimnames(2),name,ny) 
    464464    !      
    465     CALL write_ncdf_var('nav_lon'         ,dimnames,name,Grid%nav_lon    ,'float') 
    466     CALL write_ncdf_var('nav_lat'         ,dimnames,name,Grid%nav_lat    ,'float') 
    467     CALL write_ncdf_var(parent_batmet_name,dimnames,name,Grid%bathy_meter,'float') 
    468     CALL write_ncdf_var('weight'          ,dimnames,name,Grid%wgt        ,'float') 
    469     ! 
    470     CALL copy_ncdf_att('nav_lon'         ,TRIM(parent_bathy_meter),name,MINVAL(Grid%nav_lon),MAXVAL(Grid%nav_lon)) 
    471     CALL copy_ncdf_att('nav_lat'         ,TRIM(parent_bathy_meter),name,MINVAL(Grid%nav_lat),MAXVAL(Grid%nav_lat))   
    472     CALL copy_ncdf_att(parent_batmet_name,TRIM(parent_bathy_meter),name) 
     465    CALL write_ncdf_var('nav_lon'        ,dimnames,name,Grid%nav_lon    ,'float') 
     466    CALL write_ncdf_var('nav_lat'        ,dimnames,name,Grid%nav_lat    ,'float') 
     467    CALL write_ncdf_var(parent_meter_name,dimnames,name,Grid%bathy_meter,'float') 
     468    CALL write_ncdf_var('weight'         ,dimnames,name,Grid%wgt        ,'float') 
     469    ! 
     470    CALL copy_ncdf_att('nav_lon'        ,TRIM(parent_bathy_meter),name,MINVAL(Grid%nav_lon),MAXVAL(Grid%nav_lon)) 
     471    CALL copy_ncdf_att('nav_lat'        ,TRIM(parent_bathy_meter),name,MINVAL(Grid%nav_lat),MAXVAL(Grid%nav_lat))   
     472    CALL copy_ncdf_att(parent_meter_name,TRIM(parent_bathy_meter),name) 
    473473    ! 
    474474    WRITE(*,*) ' ' 
     
    500500     ! 
    501501     CHARACTER(len=1),     DIMENSION(3)     ::   dimnames 
    502      REAL*8 ,              DIMENSION(N)     ::   e3t_1d, e3w_1d, gdept_1d 
    503      REAL*8 , ALLOCATABLE, DIMENSION(:,:)   ::   ff_t, ff_f 
    504      INTEGER, ALLOCATABLE, DIMENSION(:,:)   ::   bottom_level, top_level, zbathy 
     502     REAL*8 ,              DIMENSION(N)     ::   e3t_1d, e3w_1d, gdept_1d, gdepw_1d 
     503     REAL*8 , ALLOCATABLE, DIMENSION(:,:)   ::   ff_t, ff_f, zbathy 
     504     INTEGER, ALLOCATABLE, DIMENSION(:,:)   ::   bottom_level, top_level, mbathy 
    505505     REAL*8 , ALLOCATABLE, DIMENSION(:,:,:) ::   e3t_0, e3u_0, e3v_0, e3f_0, e3w_0, e3uw_0, e3vw_0 
    506506     ! 
     
    510510 
    511511     ! allocate needed arrays for domain_cfg 
    512      ALLOCATE( ff_t(nx,ny), ff_f(nx,ny) ) 
    513      ALLOCATE( bottom_level(nx,ny), top_level(nx,ny), zbathy(nx,ny) ) 
     512     ALLOCATE( ff_t(nx,ny), ff_f(nx,ny), zbathy(nx,ny) ) 
     513     ALLOCATE( bottom_level(nx,ny), top_level(nx,ny), mbathy(nx,ny) ) 
    514514     ALLOCATE( e3t_0(nx,ny,N), e3u_0 (nx,ny,N), e3v_0 (nx,ny,N), e3f_0(nx,ny,N),  & 
    515515        &      e3w_0(nx,ny,N), e3uw_0(nx,ny,N), e3vw_0(nx,ny,N) ) 
     
    527527     ff_t(:,:) = 2. * omega * SIN( rad * Grid%gphit(:,:) )     !    -        -       -    at t-point 
    528528 
    529      ! top/bottom levels 
    530      bottom_level(:,:) = Grid%bathy_level(:,:) 
    531      top_level(:,:)    = MIN( 1, bottom_level(:,:) ) 
     529     ! bathy in meters 
     530     zbathy(:,:) = Grid%bathy_meter 
    532531 
    533532     ! vertical scale factors 
    534      CALL zgr_z( e3t_1d, e3w_1d, gdept_1d )     
    535      DO jk = 1, N 
    536         e3t_0  (:,:,jk) = e3t_1d  (jk) 
    537         e3u_0  (:,:,jk) = e3t_1d  (jk) 
    538         e3v_0  (:,:,jk) = e3t_1d  (jk) 
    539         e3f_0  (:,:,jk) = e3t_1d  (jk) 
    540         e3w_0  (:,:,jk) = e3w_1d  (jk) 
    541         e3uw_0 (:,:,jk) = e3w_1d  (jk) 
    542         e3vw_0 (:,:,jk) = e3w_1d  (jk) 
    543      END DO 
     533     CALL zgr_z( e3t_1d, e3w_1d, gdept_1d, gdepw_1d )     
     534!     DO jk = 1, N 
     535!        e3t_0  (:,:,jk) = e3t_1d  (jk) 
     536!        e3u_0  (:,:,jk) = e3t_1d  (jk) 
     537!        e3v_0  (:,:,jk) = e3t_1d  (jk) 
     538!        e3f_0  (:,:,jk) = e3t_1d  (jk) 
     539!        e3w_0  (:,:,jk) = e3w_1d  (jk) 
     540!        e3uw_0 (:,:,jk) = e3w_1d  (jk) 
     541!        e3vw_0 (:,:,jk) = e3w_1d  (jk) 
     542!     END DO 
    544543 
    545544     ! logicals and others 
     
    549548        ln_zps = 1 
    550549        ln_zco = 0 
     550        CALL zgr_zps( nx, ny, gdept_1d, gdepw_1d, e3t_1d, e3w_1d, zbathy, &        ! input 
     551        &             mbathy, e3t_0, e3u_0, e3v_0, e3f_0, e3w_0, e3uw_0, e3vw_0 )  ! output 
    551552     ELSE 
    552553        ln_zps = 0 
     
    554555     ENDIF 
    555556 
     557     ! top/bottom levels 
     558     bottom_level(:,:) = mbathy(:,:) 
     559     top_level(:,:)    = MIN( 1, mbathy(:,:) ) 
     560 
    556561     ! closed domain (agrif) 
    557562     jperio = 0 
     
    567572     top_level(  nx,1:ny) = 0 
    568573 
    569      zbathy(:,:) = Grid%bathy_meter 
    570574     zbathy(1:nx,1   ) = 0. 
    571575     zbathy(1:nx,  ny) = 0. 
     
    576580     ! write domain_cfg.nc 
    577581     ! ------------------- 
    578      status = nf90_create(name,NF90_WRITE,ncid) 
     582     status = nf90_create(name,IOR(NF90_64BIT_OFFSET,NF90_WRITE),ncid) 
    579583     status = nf90_close(ncid) 
    580584     ! 
     
    664668     WRITE(*,*) ' ' 
    665669     ! 
    666      DEALLOCATE( ff_t, ff_f ) 
    667      DEALLOCATE( bottom_level, top_level, zbathy ) 
     670     DEALLOCATE( ff_t, ff_f, zbathy ) 
     671     DEALLOCATE( bottom_level, top_level, mbathy ) 
    668672     DEALLOCATE( e3t_0, e3u_0, e3v_0, e3f_0, e3w_0, e3uw_0, e3vw_0 ) 
    669673     !       
     
    672676  END FUNCTION write_domcfg 
    673677  ! 
    674   SUBROUTINE zgr_z(e3t_1d, e3w_1d, gdept_1d) 
     678  SUBROUTINE zgr_z(e3t_1d, e3w_1d, gdept_1d, gdepw_1d ) 
    675679     !!---------------------------------------------------------------------- 
    676680     !!                   ***  ROUTINE zgr_z (from NEMO4) *** 
     
    694698     !!---------------------------------------------------------------------- 
    695699     INTEGER  ::   jk                   ! dummy loop indices 
    696      INTEGER  ::   jpk 
    697700     REAL*8 ::   zt, zw                 ! temporary scalars 
    698701     REAL*8 ::   zsur, za0, za1, zkth   ! Values set from parameters in 
     
    700703     REAL*8 ::   za2, zkth2, zacr2      ! Values for optional double tanh function set from parameters  
    701704 
    702      REAL*8, DIMENSION(:), INTENT(inout) ::  e3t_1d, e3w_1d, gdept_1d 
    703      REAL*8, DIMENSION(N)                ::  gdepw_1d 
     705     REAL*8, DIMENSION(:), INTENT(out) ::  e3t_1d, e3w_1d, gdept_1d, gdepw_1d 
    704706     !!---------------------------------------------------------------------- 
    705707     ! 
     
    788790  END SUBROUTINE zgr_z 
    789791 
     792  SUBROUTINE zgr_zps( nx, ny, gdept_1d, gdepw_1d, e3t_1d, e3w_1d, zbathy, & 
     793     &                mbathy, e3t_0, e3u_0, e3v_0, e3f_0, e3w_0, e3uw_0, e3vw_0 ) 
     794      !!---------------------------------------------------------------------- 
     795      !!                  ***  ROUTINE zgr_zps (from NEMO4)  *** 
     796      !!                      
     797      !! ** Purpose :   the depth and vertical scale factor in partial step 
     798      !!              reference z-coordinate case 
     799      !! 
     800      !! ** Method  :   Partial steps : computes the 3D vertical scale factors 
     801      !!      of T-, U-, V-, W-, UW-, VW and F-points that are associated with 
     802      !!      a partial step representation of bottom topography. 
     803      !! 
     804      !!        The reference depth of model levels is defined from an analytical 
     805      !!      function the derivative of which gives the reference vertical 
     806      !!      scale factors. 
     807      !!        From  depth and scale factors reference, we compute there new value 
     808      !!      with partial steps  on 3d arrays ( i, j, k ). 
     809      !! 
     810      !!              w-level: gdepw_0(i,j,k)  = gdep(k) 
     811      !!                       e3w_0(i,j,k) = dk(gdep)(k)     = e3(i,j,k) 
     812      !!              t-level: gdept_0(i,j,k)  = gdep(k+0.5) 
     813      !!                       e3t_0(i,j,k) = dk(gdep)(k+0.5) = e3(i,j,k+0.5) 
     814      !! 
     815      !!        With the help of the bathymetric file ( bathymetry_depth_ORCA_R2.nc), 
     816      !!      we find the mbathy index of the depth at each grid point. 
     817      !!      This leads us to three cases: 
     818      !! 
     819      !!              - bathy = 0 => mbathy = 0 
     820      !!              - 1 < mbathy < jpkm1     
     821      !!              - bathy > gdepw_0(jpk) => mbathy = jpkm1   
     822      !! 
     823      !!        Then, for each case, we find the new depth at t- and w- levels 
     824      !!      and the new vertical scale factors at t-, u-, v-, w-, uw-, vw-  
     825      !!      and f-points. 
     826      !!  
     827      !!        This routine is given as an example, it must be modified 
     828      !!      following the user s desiderata. nevertheless, the output as 
     829      !!      well as the way to compute the model levels and scale factors 
     830      !!      must be respected in order to insure second order accuracy 
     831      !!      schemes. 
     832      !! 
     833      !!         c a u t i o n : gdept_1d, gdepw_1d and e3._1d are positives 
     834      !!         - - - - - - -   gdept_0, gdepw_0 and e3. are positives 
     835      !!       
     836      !!  Reference :   Pacanowsky & Gnanadesikan 1997, Mon. Wea. Rev., 126, 3248-3270. 
     837      !!---------------------------------------------------------------------- 
     838      INTEGER  ::   ji, jj, jk       ! dummy loop indices 
     839      INTEGER  ::   ik 
     840      INTEGER  ::   ikb, ikt       ! temporary integers 
     841      REAL*8 ::   ze3tp , ze3wp    ! Last ocean level thickness at T- and W-points 
     842      REAL*8 ::   zdepwp, zdepth   ! Ajusted ocean depth to avoid too small e3t 
     843      REAL*8 ::   zmax             ! temporary scalar 
     844      ! 
     845      REAL*8, ALLOCATABLE, DIMENSION(:,:,:) ::  gdept_0, gdepw_0, zprt 
     846      ! 
     847      INTEGER,                   INTENT(in   ) ::  nx, ny 
     848      REAL*8 , DIMENSION(:)    , INTENT(in   ) ::  gdept_1d, gdepw_1d, e3t_1d, e3w_1d 
     849      REAL*8 , DIMENSION(:,:)  , INTENT(inout) ::  zbathy 
     850      INTEGER, DIMENSION(:,:)  , INTENT(  out) ::  mbathy 
     851      REAL*8 , DIMENSION(:,:,:), INTENT(  out) ::  e3t_0, e3u_0, e3v_0, e3f_0, e3w_0, e3uw_0, e3vw_0 
     852      ! 
     853      !!--------------------------------------------------------------------- 
     854       
     855      ALLOCATE( zprt(nx,ny,N), gdept_0(nx,ny,N), gdepw_0(nx,ny,N) ) 
     856      ! 
     857      ! bathymetry in level (from bathy_meter) 
     858      ! =================== 
     859      zmax = gdepw_1d(N) + e3t_1d(N)        ! maximum depth (i.e. the last ocean level thickness <= 2*e3t_1d(N-1) ) 
     860      zbathy(:,:) = MIN( zmax ,  zbathy(:,:) )    ! bounded value of bathy (min already set at the end of zgr_bat) 
     861      WHERE( zbathy(:,:) == 0. )   ;   mbathy(:,:) = 0       ! land  : set mbathy to 0 
     862      ELSEWHERE                    ;   mbathy(:,:) = N-1   ! ocean : initialize mbathy to the max ocean level 
     863      END WHERE 
     864 
     865      ! Compute mbathy for ocean points (i.e. the number of ocean levels) 
     866      ! find the number of ocean levels such that the last level thickness 
     867      ! is larger than the minimum of e3zps_min and e3zps_rat * e3t_1d (where 
     868      ! e3t_1d is the reference level thickness 
     869      DO jk = N-1, 1, -1 
     870         zdepth = gdepw_1d(jk) + MIN( e3zps_min, e3t_1d(jk)*e3zps_rat ) 
     871         WHERE( 0. < zbathy(:,:) .AND. zbathy(:,:) <= zdepth )   mbathy(:,:) = jk-1 
     872      END DO 
     873 
     874      ! Scale factors and depth at T- and W-points 
     875      DO jk = 1, N                        ! intitialization to the reference z-coordinate 
     876         gdept_0(:,:,jk) = gdept_1d(jk) 
     877         gdepw_0(:,:,jk) = gdepw_1d(jk) 
     878         e3t_0  (:,:,jk) = e3t_1d  (jk) 
     879         e3w_0  (:,:,jk) = e3w_1d  (jk) 
     880      END DO 
     881       
     882      ! Bathy, iceshelf draft, scale factor and depth at T- and W- points in case of isf 
     883      !!clem: not implemented yet 
     884      !!      IF ( ln_isfcav == 1 ) CALL zgr_isf 
     885      ! 
     886      ! Scale factors and depth at T- and W-points 
     887!      IF ( ln_isfcav == 0 ) THEN 
     888         DO jj = 1, ny 
     889            DO ji = 1, nx 
     890               ik = mbathy(ji,jj) 
     891               IF( ik > 0 ) THEN               ! ocean point only 
     892                  ! max ocean level case 
     893                  IF( ik == N-1 ) THEN 
     894                     zdepwp = zbathy(ji,jj) 
     895                     ze3tp  = zbathy(ji,jj) - gdepw_1d(ik) 
     896                     ze3wp = 0.5 * e3w_1d(ik) * ( 1. + ( ze3tp/e3t_1d(ik) ) ) 
     897                     e3t_0(ji,jj,ik  ) = ze3tp 
     898                     e3t_0(ji,jj,ik+1) = ze3tp 
     899                     e3w_0(ji,jj,ik  ) = ze3wp 
     900                     e3w_0(ji,jj,ik+1) = ze3tp 
     901                     gdepw_0(ji,jj,ik+1) = zdepwp 
     902                     gdept_0(ji,jj,ik  ) = gdept_1d(ik-1) + ze3wp 
     903                     gdept_0(ji,jj,ik+1) = gdept_0(ji,jj,ik) + ze3tp 
     904                     ! 
     905                  ELSE                         ! standard case 
     906                     IF( zbathy(ji,jj) <= gdepw_1d(ik+1) ) THEN  ;   gdepw_0(ji,jj,ik+1) = zbathy(ji,jj) 
     907                     ELSE                                        ;   gdepw_0(ji,jj,ik+1) = gdepw_1d(ik+1) 
     908                     ENDIF 
     909   !gm Bug?  check the gdepw_1d 
     910                     !       ... on ik 
     911                     gdept_0(ji,jj,ik) = gdepw_1d(ik) + ( gdepw_0(ji,jj,ik+1) - gdepw_1d(ik) )   & 
     912                        &                             * ((gdept_1d(     ik  ) - gdepw_1d(ik) )   & 
     913                        &                             / ( gdepw_1d(     ik+1) - gdepw_1d(ik) )) 
     914                     e3t_0  (ji,jj,ik) = e3t_1d  (ik) * ( gdepw_0 (ji,jj,ik+1) - gdepw_1d(ik) )   &  
     915                        &                             / ( gdepw_1d(      ik+1) - gdepw_1d(ik) )  
     916                     e3w_0(ji,jj,ik) = 0.5 * ( gdepw_0(ji,jj,ik+1) + gdepw_1d(ik+1) - 2. * gdepw_1d(ik) )   & 
     917                        &                     * ( e3w_1d(ik) / ( gdepw_1d(ik+1) - gdepw_1d(ik) ) ) 
     918                     !       ... on ik+1 
     919                     e3w_0  (ji,jj,ik+1) = e3t_0  (ji,jj,ik) 
     920                     e3t_0  (ji,jj,ik+1) = e3t_0  (ji,jj,ik) 
     921                     gdept_0(ji,jj,ik+1) = gdept_0(ji,jj,ik) + e3t_0(ji,jj,ik) 
     922                  ENDIF 
     923               ENDIF 
     924            END DO 
     925         END DO 
     926         ! 
     927!         DO jj = 1, ny 
     928!            DO ji = 1, nx 
     929!               ik = mbathy(ji,jj) 
     930!               IF( ik > 0 ) THEN               ! ocean point only 
     931!                  e3tp (ji,jj) = e3t_0(ji,jj,ik) 
     932!                  e3wp (ji,jj) = e3w_0(ji,jj,ik) 
     933!               ENDIF 
     934!            END DO 
     935!         END DO 
     936!      END IF 
     937      ! 
     938      ! Scale factors and depth at U-, V-, UW and VW-points 
     939      DO jk = 1, N                        ! initialisation to z-scale factors 
     940         e3u_0 (:,:,jk) = e3t_1d(jk) 
     941         e3v_0 (:,:,jk) = e3t_1d(jk) 
     942         e3uw_0(:,:,jk) = e3w_1d(jk) 
     943         e3vw_0(:,:,jk) = e3w_1d(jk) 
     944      END DO 
     945 
     946      DO jk = 1,N                         ! Computed as the minimum of neighbooring scale factors 
     947         DO jj = 1, ny-1 
     948            DO ji = 1, nx-1   ! vector opt. 
     949               e3u_0 (ji,jj,jk) = MIN( e3t_0(ji,jj,jk), e3t_0(ji+1,jj,jk) ) 
     950               e3v_0 (ji,jj,jk) = MIN( e3t_0(ji,jj,jk), e3t_0(ji,jj+1,jk) ) 
     951               e3uw_0(ji,jj,jk) = MIN( e3w_0(ji,jj,jk), e3w_0(ji+1,jj,jk) ) 
     952               e3vw_0(ji,jj,jk) = MIN( e3w_0(ji,jj,jk), e3w_0(ji,jj+1,jk) ) 
     953            END DO 
     954         END DO 
     955      END DO 
     956       
     957!      IF ( ln_isfcav == 1 ) THEN 
     958!      ! (ISF) define e3uw (adapted for 2 cells in the water column) 
     959!         DO jj = 2, ny-1  
     960!            DO ji = 2, nx-1   ! vector opt.  
     961!               ikb = MAX(mbathy (ji,jj),mbathy (ji+1,jj)) 
     962!               ikt = MAX(misfdep(ji,jj),misfdep(ji+1,jj)) 
     963!               IF (ikb == ikt+1) e3uw_0(ji,jj,ikb) =  MIN( gdept_0(ji,jj,ikb  ), gdept_0(ji+1,jj  ,ikb  ) ) & 
     964!                                       &            - MAX( gdept_0(ji,jj,ikb-1), gdept_0(ji+1,jj  ,ikb-1) ) 
     965!               ikb = MAX(mbathy (ji,jj),mbathy (ji,jj+1)) 
     966!               ikt = MAX(misfdep(ji,jj),misfdep(ji,jj+1)) 
     967!               IF (ikb == ikt+1) e3vw_0(ji,jj,ikb) =  MIN( gdept_0(ji,jj,ikb  ), gdept_0(ji  ,jj+1,ikb  ) ) & 
     968!                                       &            - MAX( gdept_0(ji,jj,ikb-1), gdept_0(ji  ,jj+1,ikb-1) ) 
     969!            END DO 
     970!         END DO 
     971!      END IF 
     972      ! 
     973 
     974      DO jk = 1, N                        ! set to z-scale factor if zero (i.e. along closed boundaries) 
     975         WHERE( e3u_0 (:,:,jk) == 0. )   e3u_0 (:,:,jk) = e3t_1d(jk) 
     976         WHERE( e3v_0 (:,:,jk) == 0. )   e3v_0 (:,:,jk) = e3t_1d(jk) 
     977         WHERE( e3uw_0(:,:,jk) == 0. )   e3uw_0(:,:,jk) = e3w_1d(jk) 
     978         WHERE( e3vw_0(:,:,jk) == 0. )   e3vw_0(:,:,jk) = e3w_1d(jk) 
     979      END DO 
     980       
     981      ! Scale factor at F-point 
     982      DO jk = 1, N                        ! initialisation to z-scale factors 
     983         e3f_0(:,:,jk) = e3t_1d(jk) 
     984      END DO 
     985      DO jk = 1, N                        ! Computed as the minimum of neighbooring V-scale factors 
     986         DO jj = 1, ny-1 
     987            DO ji = 1, nx-1   ! vector opt. 
     988               e3f_0(ji,jj,jk) = MIN( e3v_0(ji,jj,jk), e3v_0(ji+1,jj,jk) ) 
     989            END DO 
     990         END DO 
     991      END DO 
     992      ! 
     993      DO jk = 1, N                        ! set to z-scale factor if zero (i.e. along closed boundaries) 
     994         WHERE( e3f_0(:,:,jk) == 0. )   e3f_0(:,:,jk) = e3t_1d(jk) 
     995      END DO 
     996 
     997      ! Control of the sign 
     998      IF( MINVAL( e3t_0  (:,:,:) ) <= 0. )   STOP '    zgr_zps :   e r r o r   e3t_0 <= 0'  
     999      IF( MINVAL( e3w_0  (:,:,:) ) <= 0. )   STOP '    zgr_zps :   e r r o r   e3w_0 <= 0'  
     1000      IF( MINVAL( gdept_0(:,:,:) ) <  0. )   STOP '    zgr_zps :   e r r o r   gdept_0 <  0'  
     1001      IF( MINVAL( gdepw_0(:,:,:) ) <  0. )   STOP '    zgr_zps :   e r r o r   gdepw_0 <  0'  
     1002      
     1003      ! Compute gde3w_0 (vertical sum of e3w) 
     1004!      IF ( ln_isfcav ==1 ) THEN ! if cavity 
     1005!         WHERE( misfdep == 0 )   misfdep = 1 
     1006!         DO jj = 1,ny 
     1007!            DO ji = 1,nx 
     1008!               gde3w_0(ji,jj,1) = 0.5 * e3w_0(ji,jj,1) 
     1009!               DO jk = 2, misfdep(ji,jj) 
     1010!                  gde3w_0(ji,jj,jk) = gde3w_0(ji,jj,jk-1) + e3w_0(ji,jj,jk)  
     1011!               END DO 
     1012!               IF( misfdep(ji,jj) >= 2 )   gde3w_0(ji,jj,misfdep(ji,jj)) = risfdep(ji,jj) + 0.5 * e3w_0(ji,jj,misfdep(ji,jj)) 
     1013!               DO jk = misfdep(ji,jj) + 1, N 
     1014!                  gde3w_0(ji,jj,jk) = gde3w_0(ji,jj,jk-1) + e3w_0(ji,jj,jk)  
     1015!               END DO 
     1016!            END DO 
     1017!         END DO 
     1018!      ELSE ! no cavity 
     1019!         gde3w_0(:,:,1) = 0.5 * e3w_0(:,:,1) 
     1020!         DO jk = 2, N 
     1021!            gde3w_0(:,:,jk) = gde3w_0(:,:,jk-1) + e3w_0(:,:,jk) 
     1022!         END DO 
     1023!      END IF 
     1024      ! 
     1025      DEALLOCATE( zprt, gdept_0, gdepw_0 ) 
     1026      ! 
     1027   END SUBROUTINE zgr_zps 
    7901028 
    7911029  !***************************************************** 
     
    8281066    LOGICAL :: find         
    8291067    i=1 
    830     DO WHILE ( TRIM(VAR_INTERP(i)) .NE. 'NULL' )      
     1068    DO WHILE ( TRIM(VAR_INTERP(i)) /= '' )      
    8311069       pos = INDEX( TRIM(VAR_INTERP(i)) , TRIM(varname) ) 
    8321070       IF ( pos .NE. 0 ) THEN       
Note: See TracChangeset for help on using the changeset viewer.