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 7646 for trunk/NEMOGCM/NEMO/OPA_SRC/DOM/domwri.F90 – NEMO

Ignore:
Timestamp:
2017-02-06T10:25:03+01:00 (7 years ago)
Author:
timgraham
Message:

Merge of dev_merge_2016 into trunk. UPDATE TO ARCHFILES NEEDED for XIOS2.
LIM_SRC_s/limrhg.F90 to follow in next commit due to change of kind (I'm unable to do it in this commit).
Merged using the following steps:

1) svn merge --reintegrate ​svn+ssh://forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/trunk .
2) Resolve minor conflicts in sette.sh and namelist_cfg for ORCA2LIM3 (due to a change in trunk after branch was created)
3) svn commit
4) svn switch ​svn+ssh://forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/trunk
5) svn merge ​svn+ssh://forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/branches/2016/dev_merge_2016 .
6) At this stage I checked out a clean copy of the branch to compare against what is about to be committed to the trunk.
6) svn commit #Commit code to the trunk

In this commit I have also reverted a change to Fcheck_archfile.sh which was causing problems on the Paris machine.

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/NEMOGCM/NEMO/OPA_SRC/DOM/domwri.F90

    r6689 r7646  
    88   !!   NEMO     1.0  ! 2002-08  (G. Madec)  F90 and several file 
    99   !!            3.0  ! 2008-01  (S. Masson)  add dom_uniq  
     10   !!            4.0  ! 2016-01  (G. Madec)  simplified mesh_mask.nc file 
    1011   !!---------------------------------------------------------------------- 
    1112 
    … …  
    1314   !!   dom_wri        : create and write mesh and mask file(s) 
    1415   !!   dom_uniq       : identify unique point of a grid (TUVF) 
     16   !!   dom_stiff      : diagnose maximum grid stiffness/hydrostatic consistency (s-coordinate) 
    1517   !!---------------------------------------------------------------------- 
    1618   USE dom_oce         ! ocean space and time domain 
     19   USE phycst ,   ONLY :   rsmall 
     20   USE wet_dry,   ONLY :   ln_wd, ht_wd 
     21   ! 
    1722   USE in_out_manager  ! I/O manager 
    1823   USE iom             ! I/O library 
    … …  
    2631 
    2732   PUBLIC   dom_wri              ! routine called by inidom.F90 
    28    PUBLIC   dom_wri_coordinate   ! routine called by domhgr.F90 
     33   PUBLIC   dom_stiff            ! routine called by inidom.F90 
     34 
    2935   !! * Substitutions 
    3036#  include "vectopt_loop_substitute.h90" 
    3137   !!---------------------------------------------------------------------- 
    32    !! NEMO/OPA 3.3 , NEMO Consortium (2010) 
     38   !! NEMO/OPA 4.0 , NEMO Consortium (2016) 
    3339   !! $Id$  
    3440   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
    3541   !!---------------------------------------------------------------------- 
    3642CONTAINS 
    37  
    38    SUBROUTINE dom_wri_coordinate 
    39       !!---------------------------------------------------------------------- 
    40       !!                  ***  ROUTINE dom_wri_coordinate  *** 
    41       !!                    
    42       !! ** Purpose :   Create the NetCDF file which contains all the 
    43       !!              standard coordinate information plus the surface, 
    44       !!              e1e2u and e1e2v. By doing so, those surface will 
    45       !!              not be changed by the reduction of e1u or e2v scale  
    46       !!              factors in some straits.  
    47       !!                 NB: call just after the read of standard coordinate 
    48       !!              and the reduction of scale factors in some straits 
    49       !! 
    50       !! ** output file :   coordinate_e1e2u_v.nc 
    51       !!---------------------------------------------------------------------- 
    52       INTEGER           ::   inum0    ! temprary units for 'coordinate_e1e2u_v.nc' file 
    53       CHARACTER(len=21) ::   clnam0   ! filename (mesh and mask informations) 
    54       !                                   !  workspaces 
    55       REAL(wp), POINTER, DIMENSION(:,:  ) :: zprt, zprw  
    56       REAL(wp), POINTER, DIMENSION(:,:,:) :: zdepu, zdepv 
    57       !!---------------------------------------------------------------------- 
    58       ! 
    59       IF( nn_timing == 1 )  CALL timing_start('dom_wri_coordinate') 
    60       ! 
    61       IF(lwp) WRITE(numout,*) 
    62       IF(lwp) WRITE(numout,*) 'dom_wri_coordinate : create NetCDF coordinate file' 
    63       IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~~~~~' 
    64        
    65       clnam0 = 'coordinate_e1e2u_v'  ! filename (mesh and mask informations) 
    66        
    67       !  create 'coordinate_e1e2u_v.nc' file 
    68       ! ============================ 
    69       ! 
    70       CALL iom_open( TRIM(clnam0), inum0, ldwrt = .TRUE., kiolib = jprstlib ) 
    71       ! 
    72       !                                                         ! horizontal mesh (inum3) 
    73       CALL iom_rstput( 0, 0, inum0, 'glamt', glamt, ktype = jp_r8 )     !    ! latitude 
    74       CALL iom_rstput( 0, 0, inum0, 'glamu', glamu, ktype = jp_r8 ) 
    75       CALL iom_rstput( 0, 0, inum0, 'glamv', glamv, ktype = jp_r8 ) 
    76       CALL iom_rstput( 0, 0, inum0, 'glamf', glamf, ktype = jp_r8 ) 
    77        
    78       CALL iom_rstput( 0, 0, inum0, 'gphit', gphit, ktype = jp_r8 )     !    ! longitude 
    79       CALL iom_rstput( 0, 0, inum0, 'gphiu', gphiu, ktype = jp_r8 ) 
    80       CALL iom_rstput( 0, 0, inum0, 'gphiv', gphiv, ktype = jp_r8 ) 
    81       CALL iom_rstput( 0, 0, inum0, 'gphif', gphif, ktype = jp_r8 ) 
    82        
    83       CALL iom_rstput( 0, 0, inum0, 'e1t', e1t, ktype = jp_r8 )         !    ! e1 scale factors 
    84       CALL iom_rstput( 0, 0, inum0, 'e1u', e1u, ktype = jp_r8 ) 
    85       CALL iom_rstput( 0, 0, inum0, 'e1v', e1v, ktype = jp_r8 ) 
    86       CALL iom_rstput( 0, 0, inum0, 'e1f', e1f, ktype = jp_r8 ) 
    87        
    88       CALL iom_rstput( 0, 0, inum0, 'e2t', e2t, ktype = jp_r8 )         !    ! e2 scale factors 
    89       CALL iom_rstput( 0, 0, inum0, 'e2u', e2u, ktype = jp_r8 ) 
    90       CALL iom_rstput( 0, 0, inum0, 'e2v', e2v, ktype = jp_r8 ) 
    91       CALL iom_rstput( 0, 0, inum0, 'e2f', e2f, ktype = jp_r8 ) 
    92        
    93       CALL iom_rstput( 0, 0, inum0, 'e1e2u', e1e2u, ktype = jp_r8 ) 
    94       CALL iom_rstput( 0, 0, inum0, 'e1e2v', e1e2v, ktype = jp_r8 ) 
    95  
    96       CALL iom_close( inum0 ) 
    97       ! 
    98       IF( nn_timing == 1 )  CALL timing_stop('dom_wri_coordinate') 
    99       ! 
    100    END SUBROUTINE dom_wri_coordinate 
    101  
    10243 
    10344   SUBROUTINE dom_wri 
    … …  
    11354      !!      domhgr, domzgr, and dommsk. Note: the file contain depends on 
    11455      !!      the vertical coord. used (z-coord, partial steps, s-coord) 
    115       !!            MOD(nmsh, 3) = 1  :   'mesh_mask.nc' file 
     56      !!            MOD(nn_msh, 3) = 1  :   'mesh_mask.nc' file 
    11657      !!                         = 2  :   'mesh.nc' and mask.nc' files 
    11758      !!                         = 0  :   'mesh_hgr.nc', 'mesh_zgr.nc' and 
    … …  
    12061      !!      vertical coordinate. 
    12162      !! 
    122       !!      if     nmsh <= 3: write full 3D arrays for e3[tuvw] and gdep[tuvw] 
    123       !!      if 3 < nmsh <= 6: write full 3D arrays for e3[tuvw] and 2D arrays  
     63      !!      if     nn_msh <= 3: write full 3D arrays for e3[tuvw] and gdep[tuvw] 
     64      !!      if 3 < nn_msh <= 6: write full 3D arrays for e3[tuvw] and 2D arrays  
    12465      !!                        corresponding to the depth of the bottom t- and w-points 
    125       !!      if 6 < nmsh <= 9: write 2D arrays corresponding to the depth and the 
     66      !!      if 6 < nn_msh <= 9: write 2D arrays corresponding to the depth and the 
    12667      !!                        thickness (e3[tw]_ps) of the bottom points  
    12768      !! 
    … …  
    12970      !!                                   masks, depth and vertical scale factors 
    13071      !!---------------------------------------------------------------------- 
    131       !! 
    132       INTEGER           ::   inum0    ! temprary units for 'mesh_mask.nc' file 
    133       INTEGER           ::   inum1    ! temprary units for 'mesh.nc'      file 
    134       INTEGER           ::   inum2    ! temprary units for 'mask.nc'      file 
    135       INTEGER           ::   inum3    ! temprary units for 'mesh_hgr.nc'  file 
    136       INTEGER           ::   inum4    ! temprary units for 'mesh_zgr.nc'  file 
    137       CHARACTER(len=21) ::   clnam0   ! filename (mesh and mask informations) 
    138       CHARACTER(len=21) ::   clnam1   ! filename (mesh informations) 
    139       CHARACTER(len=21) ::   clnam2   ! filename (mask informations) 
    140       CHARACTER(len=21) ::   clnam3   ! filename (horizontal mesh informations) 
    141       CHARACTER(len=21) ::   clnam4   ! filename (vertical   mesh informations) 
     72      INTEGER           ::   inum    ! temprary units for 'mesh_mask.nc' file 
     73      CHARACTER(len=21) ::   clnam   ! filename (mesh and mask informations) 
    14274      INTEGER           ::   ji, jj, jk   ! dummy loop indices 
    143       !                                   !  workspaces 
    144       REAL(wp), POINTER, DIMENSION(:,:  ) :: zprt, zprw  
    145       REAL(wp), POINTER, DIMENSION(:,:,:) :: zdepu, zdepv 
     75      INTEGER           ::   izco, izps, isco, icav 
     76      !                                
     77      REAL(wp), POINTER, DIMENSION(:,:)   ::   zprt, zprw     ! 2D workspace 
     78      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zdepu, zdepv   ! 3D workspace 
    14679      !!---------------------------------------------------------------------- 
    14780      ! 
    14881      IF( nn_timing == 1 )  CALL timing_start('dom_wri') 
    14982      ! 
    150       CALL wrk_alloc( jpi, jpj, zprt, zprw ) 
    151       CALL wrk_alloc( jpi, jpj, jpk, zdepu, zdepv ) 
     83      CALL wrk_alloc( jpi,jpj,       zprt , zprw ) 
     84      CALL wrk_alloc( jpi,jpj,jpk,  zdepu, zdepv ) 
    15285      ! 
    15386      IF(lwp) WRITE(numout,*) 
    … …  
    15588      IF(lwp) WRITE(numout,*) '~~~~~~~' 
    15689       
    157       clnam0 = 'mesh_mask'  ! filename (mesh and mask informations) 
    158       clnam1 = 'mesh'       ! filename (mesh informations) 
    159       clnam2 = 'mask'       ! filename (mask informations) 
    160       clnam3 = 'mesh_hgr'   ! filename (horizontal mesh informations) 
    161       clnam4 = 'mesh_zgr'   ! filename (vertical   mesh informations) 
    162        
    163       SELECT CASE ( MOD(nmsh, 3) ) 
    164          !                                  ! ============================ 
    165       CASE ( 1 )                            !  create 'mesh_mask.nc' file 
    166          !                                  ! ============================ 
    167          CALL iom_open( TRIM(clnam0), inum0, ldwrt = .TRUE., kiolib = jprstlib ) 
    168          inum2 = inum0                                            ! put all the informations 
    169          inum3 = inum0                                            ! in unit inum0 
    170          inum4 = inum0 
    171           
    172          !                                  ! ============================ 
    173       CASE ( 2 )                            !  create 'mesh.nc' and  
    174          !                                  !         'mask.nc' files 
    175          !                                  ! ============================ 
    176          CALL iom_open( TRIM(clnam1), inum1, ldwrt = .TRUE., kiolib = jprstlib ) 
    177          CALL iom_open( TRIM(clnam2), inum2, ldwrt = .TRUE., kiolib = jprstlib ) 
    178          inum3 = inum1                                            ! put mesh informations  
    179          inum4 = inum1                                            ! in unit inum1  
    180          !                                  ! ============================ 
    181       CASE ( 0 )                            !  create 'mesh_hgr.nc' 
    182          !                                  !         'mesh_zgr.nc' and 
    183          !                                  !         'mask.nc'     files 
    184          !                                  ! ============================ 
    185          CALL iom_open( TRIM(clnam2), inum2, ldwrt = .TRUE., kiolib = jprstlib ) 
    186          CALL iom_open( TRIM(clnam3), inum3, ldwrt = .TRUE., kiolib = jprstlib ) 
    187          CALL iom_open( TRIM(clnam4), inum4, ldwrt = .TRUE., kiolib = jprstlib ) 
    188          ! 
    189       END SELECT 
    190        
    191       !                                                         ! masks (inum2)  
    192       CALL iom_rstput( 0, 0, inum2, 'tmask', tmask, ktype = jp_i1 )     !    ! land-sea mask 
    193       CALL iom_rstput( 0, 0, inum2, 'umask', umask, ktype = jp_i1 ) 
    194       CALL iom_rstput( 0, 0, inum2, 'vmask', vmask, ktype = jp_i1 ) 
    195       CALL iom_rstput( 0, 0, inum2, 'fmask', fmask, ktype = jp_i1 ) 
     90      clnam = 'mesh_mask'  ! filename (mesh and mask informations) 
     91       
     92      !                                  ! ============================ 
     93      !                                  !  create 'mesh_mask.nc' file 
     94      !                                  ! ============================ 
     95      CALL iom_open( TRIM(clnam), inum, ldwrt = .TRUE., kiolib = jprstlib ) 
     96      ! 
     97      !                                                         ! global domain size 
     98      CALL iom_rstput( 0, 0, inum, 'jpiglo', REAL( jpiglo, wp), ktype = jp_i4 ) 
     99      CALL iom_rstput( 0, 0, inum, 'jpjglo', REAL( jpjglo, wp), ktype = jp_i4 ) 
     100      CALL iom_rstput( 0, 0, inum, 'jpkglo', REAL( jpkglo, wp), ktype = jp_i4 ) 
     101 
     102      !                                                         ! domain characteristics 
     103      CALL iom_rstput( 0, 0, inum, 'jperio', REAL( jperio, wp), ktype = jp_i4 ) 
     104      !                                                         ! type of vertical coordinate 
     105      IF( ln_zco    ) THEN   ;   izco = 1   ;   ELSE   ;   izco = 0   ;   ENDIF 
     106      IF( ln_zps    ) THEN   ;   izps = 1   ;   ELSE   ;   izps = 0   ;   ENDIF 
     107      IF( ln_sco    ) THEN   ;   isco = 1   ;   ELSE   ;   isco = 0   ;   ENDIF 
     108      CALL iom_rstput( 0, 0, inum, 'ln_zco'   , REAL( izco, wp), ktype = jp_i4 ) 
     109      CALL iom_rstput( 0, 0, inum, 'ln_zps'   , REAL( izps, wp), ktype = jp_i4 ) 
     110      CALL iom_rstput( 0, 0, inum, 'ln_sco'   , REAL( isco, wp), ktype = jp_i4 ) 
     111      !                                                         ! ocean cavities under iceshelves 
     112      IF( ln_isfcav ) THEN   ;   icav = 1   ;   ELSE   ;   icav = 0   ;   ENDIF 
     113      CALL iom_rstput( 0, 0, inum, 'ln_isfcav', REAL( icav, wp), ktype = jp_i4 ) 
     114   
     115      !                                                         ! masks 
     116      CALL iom_rstput( 0, 0, inum, 'tmask', tmask, ktype = jp_i1 )     !    ! land-sea mask 
     117      CALL iom_rstput( 0, 0, inum, 'umask', umask, ktype = jp_i1 ) 
     118      CALL iom_rstput( 0, 0, inum, 'vmask', vmask, ktype = jp_i1 ) 
     119      CALL iom_rstput( 0, 0, inum, 'fmask', fmask, ktype = jp_i1 ) 
    196120       
    197121      CALL dom_uniq( zprw, 'T' ) 
    198122      DO jj = 1, jpj 
    199123         DO ji = 1, jpi 
    200             jk=mikt(ji,jj)  
    201             zprt(ji,jj) = tmask(ji,jj,jk) * zprw(ji,jj)                        !    ! unique point mask 
     124            zprt(ji,jj) = ssmask(ji,jj) * zprw(ji,jj)                        !    ! unique point mask 
    202125         END DO 
    203126      END DO                             !    ! unique point mask 
    204       CALL iom_rstput( 0, 0, inum2, 'tmaskutil', zprt, ktype = jp_i1 )   
     127      CALL iom_rstput( 0, 0, inum, 'tmaskutil', zprt, ktype = jp_i1 )   
    205128      CALL dom_uniq( zprw, 'U' ) 
    206129      DO jj = 1, jpj 
    207130         DO ji = 1, jpi 
    208             jk=miku(ji,jj)  
    209             zprt(ji,jj) = umask(ji,jj,jk) * zprw(ji,jj)                        !    ! unique point mask 
     131            zprt(ji,jj) = ssumask(ji,jj) * zprw(ji,jj)                        !    ! unique point mask 
    210132         END DO 
    211133      END DO 
    212       CALL iom_rstput( 0, 0, inum2, 'umaskutil', zprt, ktype = jp_i1 )   
     134      CALL iom_rstput( 0, 0, inum, 'umaskutil', zprt, ktype = jp_i1 )   
    213135      CALL dom_uniq( zprw, 'V' ) 
    214136      DO jj = 1, jpj 
    215137         DO ji = 1, jpi 
    216             jk=mikv(ji,jj)  
    217             zprt(ji,jj) = vmask(ji,jj,jk) * zprw(ji,jj)                        !    ! unique point mask 
     138            zprt(ji,jj) = ssvmask(ji,jj) * zprw(ji,jj)                        !    ! unique point mask 
    218139         END DO 
    219140      END DO 
    220       CALL iom_rstput( 0, 0, inum2, 'vmaskutil', zprt, ktype = jp_i1 )   
    221       CALL dom_uniq( zprw, 'F' ) 
    222       DO jj = 1, jpj 
    223          DO ji = 1, jpi 
    224             jk=mikf(ji,jj)  
    225             zprt(ji,jj) = fmask(ji,jj,jk) * zprw(ji,jj)                        !    ! unique point mask 
    226          END DO 
    227       END DO 
    228       CALL iom_rstput( 0, 0, inum2, 'fmaskutil', zprt, ktype = jp_i1 )   
     141      CALL iom_rstput( 0, 0, inum, 'vmaskutil', zprt, ktype = jp_i1 )   
     142!!gm  ssfmask has been removed  ==>> find another solution to defined fmaskutil 
     143!!    Here we just remove the output of fmaskutil. 
     144!      CALL dom_uniq( zprw, 'F' ) 
     145!      DO jj = 1, jpj 
     146!         DO ji = 1, jpi 
     147!            zprt(ji,jj) = ssfmask(ji,jj) * zprw(ji,jj)                        !    ! unique point mask 
     148!         END DO 
     149!      END DO 
     150!      CALL iom_rstput( 0, 0, inum, 'fmaskutil', zprt, ktype = jp_i1 )   
     151!!gm 
    229152 
    230153      !                                                         ! horizontal mesh (inum3) 
    231       CALL iom_rstput( 0, 0, inum3, 'glamt', glamt, ktype = jp_r8 )     !    ! latitude 
    232       CALL iom_rstput( 0, 0, inum3, 'glamu', glamu, ktype = jp_r8 ) 
    233       CALL iom_rstput( 0, 0, inum3, 'glamv', glamv, ktype = jp_r8 ) 
    234       CALL iom_rstput( 0, 0, inum3, 'glamf', glamf, ktype = jp_r8 ) 
    235        
    236       CALL iom_rstput( 0, 0, inum3, 'gphit', gphit, ktype = jp_r8 )     !    ! longitude 
    237       CALL iom_rstput( 0, 0, inum3, 'gphiu', gphiu, ktype = jp_r8 ) 
    238       CALL iom_rstput( 0, 0, inum3, 'gphiv', gphiv, ktype = jp_r8 ) 
    239       CALL iom_rstput( 0, 0, inum3, 'gphif', gphif, ktype = jp_r8 ) 
    240        
    241       CALL iom_rstput( 0, 0, inum3, 'e1t', e1t, ktype = jp_r8 )         !    ! e1 scale factors 
    242       CALL iom_rstput( 0, 0, inum3, 'e1u', e1u, ktype = jp_r8 ) 
    243       CALL iom_rstput( 0, 0, inum3, 'e1v', e1v, ktype = jp_r8 ) 
    244       CALL iom_rstput( 0, 0, inum3, 'e1f', e1f, ktype = jp_r8 ) 
    245        
    246       CALL iom_rstput( 0, 0, inum3, 'e2t', e2t, ktype = jp_r8 )         !    ! e2 scale factors 
    247       CALL iom_rstput( 0, 0, inum3, 'e2u', e2u, ktype = jp_r8 ) 
    248       CALL iom_rstput( 0, 0, inum3, 'e2v', e2v, ktype = jp_r8 ) 
    249       CALL iom_rstput( 0, 0, inum3, 'e2f', e2f, ktype = jp_r8 ) 
    250        
    251       CALL iom_rstput( 0, 0, inum3, 'ff', ff, ktype = jp_r8 )           !    ! coriolis factor 
     154      CALL iom_rstput( 0, 0, inum, 'glamt', glamt, ktype = jp_r8 )     !    ! latitude 
     155      CALL iom_rstput( 0, 0, inum, 'glamu', glamu, ktype = jp_r8 ) 
     156      CALL iom_rstput( 0, 0, inum, 'glamv', glamv, ktype = jp_r8 ) 
     157      CALL iom_rstput( 0, 0, inum, 'glamf', glamf, ktype = jp_r8 ) 
     158       
     159      CALL iom_rstput( 0, 0, inum, 'gphit', gphit, ktype = jp_r8 )     !    ! longitude 
     160      CALL iom_rstput( 0, 0, inum, 'gphiu', gphiu, ktype = jp_r8 ) 
     161      CALL iom_rstput( 0, 0, inum, 'gphiv', gphiv, ktype = jp_r8 ) 
     162      CALL iom_rstput( 0, 0, inum, 'gphif', gphif, ktype = jp_r8 ) 
     163       
     164      CALL iom_rstput( 0, 0, inum, 'e1t', e1t, ktype = jp_r8 )         !    ! e1 scale factors 
     165      CALL iom_rstput( 0, 0, inum, 'e1u', e1u, ktype = jp_r8 ) 
     166      CALL iom_rstput( 0, 0, inum, 'e1v', e1v, ktype = jp_r8 ) 
     167      CALL iom_rstput( 0, 0, inum, 'e1f', e1f, ktype = jp_r8 ) 
     168       
     169      CALL iom_rstput( 0, 0, inum, 'e2t', e2t, ktype = jp_r8 )         !    ! e2 scale factors 
     170      CALL iom_rstput( 0, 0, inum, 'e2u', e2u, ktype = jp_r8 ) 
     171      CALL iom_rstput( 0, 0, inum, 'e2v', e2v, ktype = jp_r8 ) 
     172      CALL iom_rstput( 0, 0, inum, 'e2f', e2f, ktype = jp_r8 ) 
     173       
     174      CALL iom_rstput( 0, 0, inum, 'ff_f', ff_f, ktype = jp_r8 )       !    ! coriolis factor 
     175      CALL iom_rstput( 0, 0, inum, 'ff_t', ff_t, ktype = jp_r8 ) 
    252176       
    253177      ! note that mbkt is set to 1 over land ==> use surface tmask 
    254178      zprt(:,:) = ssmask(:,:) * REAL( mbkt(:,:) , wp ) 
    255       CALL iom_rstput( 0, 0, inum4, 'mbathy', zprt, ktype = jp_i2 )     !    ! nb of ocean T-points 
     179      CALL iom_rstput( 0, 0, inum, 'mbathy', zprt, ktype = jp_i4 )     !    ! nb of ocean T-points 
    256180      zprt(:,:) = ssmask(:,:) * REAL( mikt(:,:) , wp ) 
    257       CALL iom_rstput( 0, 0, inum4, 'misf', zprt, ktype = jp_i2 )       !    ! nb of ocean T-points 
     181      CALL iom_rstput( 0, 0, inum, 'misf', zprt, ktype = jp_i4 )       !    ! nb of ocean T-points 
    258182      zprt(:,:) = ssmask(:,:) * REAL( risfdep(:,:) , wp ) 
    259       CALL iom_rstput( 0, 0, inum4, 'isfdraft', zprt, ktype = jp_r8 )       !    ! nb of ocean T-points 
    260              
    261       IF( ln_sco ) THEN                                         ! s-coordinate 
    262          CALL iom_rstput( 0, 0, inum4, 'hbatt', hbatt ) 
    263          CALL iom_rstput( 0, 0, inum4, 'hbatu', hbatu ) 
    264          CALL iom_rstput( 0, 0, inum4, 'hbatv', hbatv ) 
    265          CALL iom_rstput( 0, 0, inum4, 'hbatf', hbatf ) 
    266          ! 
    267          CALL iom_rstput( 0, 0, inum4, 'gsigt', gsigt )         !    ! scaling coef. 
    268          CALL iom_rstput( 0, 0, inum4, 'gsigw', gsigw )   
    269          CALL iom_rstput( 0, 0, inum4, 'gsi3w', gsi3w ) 
    270          CALL iom_rstput( 0, 0, inum4, 'esigt', esigt ) 
    271          CALL iom_rstput( 0, 0, inum4, 'esigw', esigw ) 
    272          ! 
    273          CALL iom_rstput( 0, 0, inum4, 'e3t_0', e3t_0 )         !    ! scale factors 
    274          CALL iom_rstput( 0, 0, inum4, 'e3u_0', e3u_0 ) 
    275          CALL iom_rstput( 0, 0, inum4, 'e3v_0', e3v_0 ) 
    276          CALL iom_rstput( 0, 0, inum4, 'e3w_0', e3w_0 ) 
    277          CALL iom_rstput( 0, 0, inum4, 'rx1', rx1 )             !    ! Max. grid stiffness ratio 
    278          ! 
    279          CALL iom_rstput( 0, 0, inum4, 'gdept_1d' , gdept_1d )  !    ! stretched system 
    280          CALL iom_rstput( 0, 0, inum4, 'gdepw_1d' , gdepw_1d ) 
    281          CALL iom_rstput( 0, 0, inum4, 'gdept_0', gdept_0, ktype = jp_r8 )      
    282          CALL iom_rstput( 0, 0, inum4, 'gdepw_0', gdepw_0, ktype = jp_r8 ) 
     183      CALL iom_rstput( 0, 0, inum, 'isfdraft', zprt, ktype = jp_r8 )   !    ! nb of ocean T-points 
     184      !                                                         ! vertical mesh 
     185      CALL iom_rstput( 0, 0, inum, 'e3t_0', e3t_0, ktype = jp_r8  )    !    ! scale factors 
     186      CALL iom_rstput( 0, 0, inum, 'e3u_0', e3u_0, ktype = jp_r8  ) 
     187      CALL iom_rstput( 0, 0, inum, 'e3v_0', e3v_0, ktype = jp_r8  ) 
     188      CALL iom_rstput( 0, 0, inum, 'e3w_0', e3w_0, ktype = jp_r8  ) 
     189      ! 
     190      CALL iom_rstput( 0, 0, inum, 'gdept_1d' , gdept_1d , ktype = jp_r8 )  ! stretched system 
     191      CALL iom_rstput( 0, 0, inum, 'gdepw_1d' , gdepw_1d , ktype = jp_r8 ) 
     192      CALL iom_rstput( 0, 0, inum, 'gdept_0'  , gdept_0  , ktype = jp_r8 ) 
     193      CALL iom_rstput( 0, 0, inum, 'gdepw_0'  , gdepw_0  , ktype = jp_r8 ) 
     194      ! 
     195      IF( ln_sco ) THEN                                         ! s-coordinate stiffness 
     196         CALL dom_stiff( zprt ) 
     197         CALL iom_rstput( 0, 0, inum, 'stiffness', zprt )       ! Max. grid stiffness ratio 
    283198      ENDIF 
    284        
    285       IF( ln_zps ) THEN                                         ! z-coordinate - partial steps 
    286          ! 
    287          IF( nmsh <= 6 ) THEN                                   !    ! 3D vertical scale factors 
    288             CALL iom_rstput( 0, 0, inum4, 'e3t_0', e3t_0 )          
    289             CALL iom_rstput( 0, 0, inum4, 'e3u_0', e3u_0 ) 
    290             CALL iom_rstput( 0, 0, inum4, 'e3v_0', e3v_0 ) 
    291             CALL iom_rstput( 0, 0, inum4, 'e3w_0', e3w_0 ) 
    292          ELSE                                                   !    ! 2D masked bottom ocean scale factors 
    293             DO jj = 1,jpj    
    294                DO ji = 1,jpi 
    295                   e3tp(ji,jj) = e3t_0(ji,jj,mbkt(ji,jj)) * ssmask(ji,jj) 
    296                   e3wp(ji,jj) = e3w_0(ji,jj,mbkt(ji,jj)) * ssmask(ji,jj) 
    297                END DO 
    298             END DO 
    299             CALL iom_rstput( 0, 0, inum4, 'e3t_ps', e3tp )       
    300             CALL iom_rstput( 0, 0, inum4, 'e3w_ps', e3wp ) 
    301          END IF 
    302          ! 
    303          IF( nmsh <= 3 ) THEN                                   !    ! 3D depth 
    304             CALL iom_rstput( 0, 0, inum4, 'gdept_0', gdept_0, ktype = jp_r8 )      
    305             DO jk = 1,jpk    
    306                DO jj = 1, jpjm1    
    307                   DO ji = 1, fs_jpim1   ! vector opt. 
    308                      zdepu(ji,jj,jk) = MIN( gdept_0(ji,jj,jk) , gdept_0(ji+1,jj  ,jk) ) 
    309                      zdepv(ji,jj,jk) = MIN( gdept_0(ji,jj,jk) , gdept_0(ji  ,jj+1,jk) ) 
    310                   END DO    
    311                END DO    
    312             END DO 
    313             CALL lbc_lnk( zdepu, 'U', 1. )   ;   CALL lbc_lnk( zdepv, 'V', 1. )  
    314             CALL iom_rstput( 0, 0, inum4, 'gdepu', zdepu, ktype = jp_r8 ) 
    315             CALL iom_rstput( 0, 0, inum4, 'gdepv', zdepv, ktype = jp_r8 ) 
    316             CALL iom_rstput( 0, 0, inum4, 'gdepw_0', gdepw_0, ktype = jp_r8 ) 
    317          ELSE                                                   !    ! 2D bottom depth 
    318             DO jj = 1,jpj    
    319                DO ji = 1,jpi 
    320                   zprt(ji,jj) = gdept_0(ji,jj,mbkt(ji,jj)  ) * ssmask(ji,jj) 
    321                   zprw(ji,jj) = gdepw_0(ji,jj,mbkt(ji,jj)+1) * ssmask(ji,jj) 
    322                END DO 
    323             END DO 
    324             CALL iom_rstput( 0, 0, inum4, 'hdept', zprt, ktype = jp_r8 )      
    325             CALL iom_rstput( 0, 0, inum4, 'hdepw', zprw, ktype = jp_r8 )  
    326          ENDIF 
    327          ! 
    328          CALL iom_rstput( 0, 0, inum4, 'gdept_1d', gdept_1d )   !    ! reference z-coord. 
    329          CALL iom_rstput( 0, 0, inum4, 'gdepw_1d', gdepw_1d ) 
    330          CALL iom_rstput( 0, 0, inum4, 'e3t_1d'  , e3t_1d   ) 
    331          CALL iom_rstput( 0, 0, inum4, 'e3w_1d'  , e3w_1d   ) 
    332       ENDIF 
    333        
    334       IF( ln_zco ) THEN 
    335          !                                                      ! z-coordinate - full steps 
    336          CALL iom_rstput( 0, 0, inum4, 'gdept_1d', gdept_1d )   !    ! depth 
    337          CALL iom_rstput( 0, 0, inum4, 'gdepw_1d', gdepw_1d ) 
    338          CALL iom_rstput( 0, 0, inum4, 'e3t_1d'  , e3t_1d   )   !    ! scale factors 
    339          CALL iom_rstput( 0, 0, inum4, 'e3w_1d'  , e3w_1d   ) 
     199      ! 
     200      IF( ln_wd ) THEN                                          ! wetting and drying domain 
     201         CALL iom_rstput( 0, 0, inum, 'ht_0'   , ht_0   , ktype = jp_r8 ) 
     202         CALL iom_rstput( 0, 0, inum, 'ht_wd'  , ht_wd  , ktype = jp_r8 ) 
    340203      ENDIF 
    341204      !                                     ! ============================ 
    342       !                                     !        close the files  
     205      CALL iom_close( inum )                !        close the files  
    343206      !                                     ! ============================ 
    344       SELECT CASE ( MOD(nmsh, 3) ) 
    345       CASE ( 1 )                 
    346          CALL iom_close( inum0 ) 
    347       CASE ( 2 ) 
    348          CALL iom_close( inum1 ) 
    349          CALL iom_close( inum2 ) 
    350       CASE ( 0 ) 
    351          CALL iom_close( inum2 ) 
    352          CALL iom_close( inum3 ) 
    353          CALL iom_close( inum4 ) 
    354       END SELECT 
    355207      ! 
    356208      CALL wrk_dealloc( jpi, jpj, zprt, zprw ) 
    … …  
    371223      !!                2) check which elements have been changed 
    372224      !!---------------------------------------------------------------------- 
    373       ! 
    374225      CHARACTER(len=1)        , INTENT(in   ) ::   cdgrd   !  
    375226      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   puniq   !  
    … …  
    405256   END SUBROUTINE dom_uniq 
    406257 
     258 
     259   SUBROUTINE dom_stiff( px1 ) 
     260      !!---------------------------------------------------------------------- 
     261      !!                  ***  ROUTINE dom_stiff  *** 
     262      !!                      
     263      !! ** Purpose :   Diagnose maximum grid stiffness/hydrostatic consistency 
     264      !! 
     265      !! ** Method  :   Compute Haney (1991) hydrostatic condition ratio 
     266      !!                Save the maximum in the vertical direction 
     267      !!                (this number is only relevant in s-coordinates) 
     268      !! 
     269      !!                Haney, 1991, J. Phys. Oceanogr., 21, 610-619. 
     270      !!---------------------------------------------------------------------- 
     271      REAL(wp), DIMENSION(:,:), INTENT(out), OPTIONAL ::   px1   ! stiffness 
     272      ! 
     273      INTEGER  ::   ji, jj, jk  
     274      REAL(wp) ::   zrxmax 
     275      REAL(wp), DIMENSION(4) ::   zr1 
     276      REAL(wp), DIMENSION(jpi,jpj) ::   zx1 
     277      !!---------------------------------------------------------------------- 
     278      zx1(:,:) = 0._wp 
     279      zrxmax   = 0._wp 
     280      zr1(:)   = 0._wp 
     281      ! 
     282      DO ji = 2, jpim1 
     283         DO jj = 2, jpjm1 
     284            DO jk = 1, jpkm1 
     285!!gm   remark: dk(gdepw) = e3t   ===>>>  possible simplification of the following calculation.... 
     286!!             especially since it is gde3w which is used to compute the pressure gradient 
     287!!             furthermore, I think gdept_0 should be used below instead of w point in the numerator 
     288!!             so that the ratio is computed at the same point (i.e. uw and vw) .... 
     289               zr1(1) = ABS(  ( gdepw_0(ji  ,jj,jk  )-gdepw_0(ji-1,jj,jk  )               &  
     290                    &          +gdepw_0(ji  ,jj,jk+1)-gdepw_0(ji-1,jj,jk+1) )             & 
     291                    &       / ( gdepw_0(ji  ,jj,jk  )+gdepw_0(ji-1,jj,jk  )               & 
     292                    &          -gdepw_0(ji  ,jj,jk+1)-gdepw_0(ji-1,jj,jk+1) + rsmall )  ) * umask(ji-1,jj,jk) 
     293               zr1(2) = ABS(  ( gdepw_0(ji+1,jj,jk  )-gdepw_0(ji  ,jj,jk  )               & 
     294                    &          +gdepw_0(ji+1,jj,jk+1)-gdepw_0(ji  ,jj,jk+1) )             & 
     295                    &       / ( gdepw_0(ji+1,jj,jk  )+gdepw_0(ji  ,jj,jk  )               & 
     296                    &          -gdepw_0(ji+1,jj,jk+1)-gdepw_0(ji  ,jj,jk+1) + rsmall )  ) * umask(ji  ,jj,jk) 
     297               zr1(3) = ABS(  ( gdepw_0(ji,jj+1,jk  )-gdepw_0(ji,jj  ,jk  )               & 
     298                    &          +gdepw_0(ji,jj+1,jk+1)-gdepw_0(ji,jj  ,jk+1) )             & 
     299                    &       / ( gdepw_0(ji,jj+1,jk  )+gdepw_0(ji,jj  ,jk  )               & 
     300                    &          -gdepw_0(ji,jj+1,jk+1)-gdepw_0(ji,jj  ,jk+1) + rsmall )  ) * vmask(ji,jj  ,jk) 
     301               zr1(4) = ABS(  ( gdepw_0(ji,jj  ,jk  )-gdepw_0(ji,jj-1,jk  )               & 
     302                    &          +gdepw_0(ji,jj  ,jk+1)-gdepw_0(ji,jj-1,jk+1) )             & 
     303                    &       / ( gdepw_0(ji,jj  ,jk  )+gdepw_0(ji,jj-1,jk  )               & 
     304                    &          -gdepw_0(ji,jj  ,jk+1)-gdepw_0(ji,jj-1,jk+1) + rsmall )  ) * vmask(ji,jj-1,jk) 
     305               zrxmax = MAXVAL( zr1(1:4) ) 
     306               zx1(ji,jj) = MAX( zx1(ji,jj) , zrxmax ) 
     307            END DO 
     308         END DO 
     309      END DO 
     310      CALL lbc_lnk( zx1, 'T', 1. ) 
     311      ! 
     312      IF( PRESENT( px1 ) )    px1 = zx1 
     313      ! 
     314      zrxmax = MAXVAL( zx1 ) 
     315      ! 
     316      IF( lk_mpp )   CALL mpp_max( zrxmax ) ! max over the global domain 
     317      ! 
     318      IF(lwp) THEN 
     319         WRITE(numout,*) 
     320         WRITE(numout,*) 'dom_stiff : maximum grid stiffness ratio: ', zrxmax 
     321         WRITE(numout,*) '~~~~~~~~~' 
     322      ENDIF 
     323      ! 
     324   END SUBROUTINE dom_stiff 
     325 
    407326   !!====================================================================== 
    408327END MODULE domwri 
Note: See TracChangeset for help on using the changeset viewer.