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 7806 for branches/2015/dev_r5003_MERCATOR6_CRS/NEMOGCM/NEMO/OFF_SRC/domrea.F90 – NEMO

Ignore:
Timestamp:
2017-03-17T08:46:30+01:00 (7 years ago)
Author:
cbricaud
Message:

phaze dev_r5003_MERCATOR6_CRS branch with rev7805 of 3.6_stable branch

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/2015/dev_r5003_MERCATOR6_CRS/NEMOGCM/NEMO/OFF_SRC/domrea.F90

    r7217 r7806  
    1717   USE lib_mpp         ! distributed memory computing library 
    1818 
     19   USE iom 
    1920   USE domstp          ! domain: set the time-step 
    2021 
     
    7374 
    7475      CALL dom_nam      ! read namelist ( namrun, namdom, namcla ) 
     76      CALL dom_msk      ! Masks 
     77      CALL dom_hgr      ! Horizontal grid 
    7578      CALL dom_zgr      ! Vertical mesh and bathymetry option 
    76       CALL dom_grd      ! Create a domain file 
    77  
    78      ! 
    79       ! - ML - Used in dom_vvl_sf_nxt and lateral diffusion routines 
    80       !        but could be usefull in many other routines 
     79      ! 
    8180      e12t    (:,:) = e1t(:,:) * e2t(:,:) 
    8281      e1e2t   (:,:) = e1t(:,:) * e2t(:,:) 
     
    9190      re1v_e2v(:,:) = e1v(:,:) / e2v(:,:) 
    9291      ! 
    93       hu(:,:) = 0._wp                          ! Ocean depth at U- and V-points 
    94       hv(:,:) = 0._wp 
    95       DO jk = 1, jpk 
    96          hu(:,:) = hu(:,:) + fse3u_n(:,:,jk) * umask(:,:,jk) 
    97          hv(:,:) = hv(:,:) + fse3v_n(:,:,jk) * vmask(:,:,jk) 
    98       END DO 
    99       !                                        ! Inverse of the local depth 
    100       hur(:,:) = 1._wp / ( hu(:,:) + 1._wp - umask(:,:,1) ) * umask(:,:,1) 
    101       hvr(:,:) = 1._wp / ( hv(:,:) + 1._wp - vmask(:,:,1) ) * vmask(:,:,1) 
    102  
    10392      CALL dom_stp      ! Time step 
    104       CALL dom_msk      ! Masks 
    10593      CALL dom_ctl      ! Domain control 
    10694 
     
    178166      nstocklist = nn_stocklist 
    179167      nwrite = nn_write 
    180  
    181  
    182168      !                             ! control of output frequency 
    183169      IF ( nstock == 0 .OR. nstock > nitend ) THEN 
     
    222208904   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdom in configuration namelist', lwp ) 
    223209      IF(lwm) WRITE ( numond, namdom ) 
     210 
    224211 
    225212      IF(lwp) THEN 
     
    321308   END SUBROUTINE dom_nam 
    322309 
     310   SUBROUTINE dom_msk 
     311      !!--------------------------------------------------------------------- 
     312      !!                 ***  ROUTINE dom_msk  *** 
     313      !! ** Purpose :  Read the NetCDF file(s) which contain(s) all the 
     314      !!      ocean mask informations and defines the interior domain T-mask. 
     315      !! 
     316      !! ** Method  :  Read in a file all the arrays generated in routines 
     317      !!               dommsk:   'mask.nc' file 
     318      !!              The interior ocean/land mask is computed from tmask 
     319      !!              setting to zero the duplicated row and lines due to 
     320      !!              MPP exchange halos, est-west cyclic and north fold 
     321      !!              boundary conditions. 
     322      !! 
     323      !! ** Action :   tmask_i  : interiorland/ocean mask at t-point 
     324      !!               tpol     : ??? 
     325      !!---------------------------------------------------------------------- 
     326      ! 
     327      INTEGER  ::  inum   ! local integers 
     328      INTEGER  ::   ji, jj, jk                   ! dummy loop indices 
     329      INTEGER  ::   iif, iil, ijf, ijl       ! local integers 
     330      REAL(wp), POINTER, DIMENSION(:,:) :: zmbk 
     331      ! 
     332      !!--------------------------------------------------------------------- 
     333       
     334 
     335 
     336      IF(lwp) WRITE(numout,*) 
     337      IF(lwp) WRITE(numout,*) 'dom_rea : read NetCDF mesh and mask information file(s)' 
     338      IF(lwp) WRITE(numout,*) '~~~~~~~' 
     339 
     340      CALL wrk_alloc( jpi, jpj, zmbk ) 
     341      zmbk(:,:) = 0._wp 
     342 
     343      IF(lwp) WRITE(numout,*) '          one file in "mesh_mask.nc" ' 
     344      CALL iom_open( 'mask', inum ) 
     345 
     346         !                                                         ! masks (inum2)  
     347      CALL iom_get( inum, jpdom_data, 'tmask', tmask ) 
     348      CALL iom_get( inum, jpdom_data, 'umask', umask ) 
     349      CALL iom_get( inum, jpdom_data, 'vmask', vmask ) 
     350      CALL iom_get( inum, jpdom_data, 'fmask', fmask ) 
     351 
     352      CALL lbc_lnk( tmask, 'T', 1._wp )    ! Lateral boundary conditions 
     353      CALL lbc_lnk( umask, 'U', 1._wp )       
     354      CALL lbc_lnk( vmask, 'V', 1._wp ) 
     355      CALL lbc_lnk( fmask, 'F', 1._wp ) 
     356 
     357#if defined key_c1d 
     358      ! set umask and vmask equal tmask in 1D configuration 
     359      IF(lwp) WRITE(numout,*) 
     360      IF(lwp) WRITE(numout,*) '**********  1D configuration : set umask and vmask equal tmask ********' 
     361      IF(lwp) WRITE(numout,*) '**********                                                     ********' 
     362 
     363      umask(:,:,:) = tmask(:,:,:) 
     364      vmask(:,:,:) = tmask(:,:,:) 
     365#endif 
     366 
     367#if defined key_degrad 
     368      CALL iom_get( inum, jpdom_data, 'facvolt', facvol ) 
     369#endif 
     370 
     371      CALL iom_get( inum, jpdom_data, 'mbathy', zmbk )              ! number of ocean t-points 
     372      mbathy (:,:) = INT( zmbk(:,:) ) 
     373      misfdep(:,:) = 1                                               ! ice shelf case not yet done 
     374       
     375      CALL zgr_bot_level                                             ! mbk. arrays (deepest ocean t-, u- & v-points 
     376 
     377      !                                     ! ============================ 
     378      !                                     !        close the files  
     379      !                                     ! ============================ 
     380 
     381      ! 
     382      ! Interior domain mask (used for global sum) 
     383      ! -------------------- 
     384      ssmask(:,:)  = tmask(:,:,1) 
     385      tmask_i(:,:) = tmask(:,:,1) 
     386      iif = jpreci                        ! thickness of exchange halos in i-axis 
     387      iil = nlci - jpreci + 1 
     388      ijf = jprecj                        ! thickness of exchange halos in j-axis 
     389      ijl = nlcj - jprecj + 1 
     390      ! 
     391      tmask_i( 1 :iif,   :   ) = 0._wp    ! first columns 
     392      tmask_i(iil:jpi,   :   ) = 0._wp    ! last  columns (including mpp extra columns) 
     393      tmask_i(   :   , 1 :ijf) = 0._wp    ! first rows 
     394      tmask_i(   :   ,ijl:jpj) = 0._wp    ! last  rows (including mpp extra rows) 
     395      ! 
     396      !                                   ! north fold mask 
     397      tpol(1:jpiglo) = 1._wp 
     398      !                                 
     399      IF( jperio == 3 .OR. jperio == 4 )   tpol(jpiglo/2+1:jpiglo) = 0._wp    ! T-point pivot 
     400      IF( jperio == 5 .OR. jperio == 6 )   tpol(     1    :jpiglo) = 0._wp    ! F-point pivot 
     401      IF( jperio == 3 .OR. jperio == 4 ) THEN      ! T-point pivot: only half of the nlcj-1 row 
     402         IF( mjg(ijl-1) == jpjglo-1 ) THEN 
     403            DO ji = iif+1, iil-1 
     404               tmask_i(ji,ijl-1) = tmask_i(ji,ijl-1) * tpol(mig(ji)) 
     405            END DO 
     406         ENDIF 
     407      ENDIF  
     408      ! 
     409      ! (ISF) MIN(1,SUM(umask)) is here to check if you have effectively at 
     410      ! least 1 wet u point 
     411      DO jj = 1, jpjm1 
     412         DO ji = 1, fs_jpim1   ! vector loop 
     413            umask_i(ji,jj)  = ssmask(ji,jj) * ssmask(ji+1,jj  )  * MIN(1._wp,SUM(umask(ji,jj,:))) 
     414            vmask_i(ji,jj)  = ssmask(ji,jj) * ssmask(ji  ,jj+1)  * MIN(1._wp,SUM(vmask(ji,jj,:))) 
     415         END DO 
     416         DO ji = 1, jpim1      ! NO vector opt. 
     417            fmask_i(ji,jj) =  ssmask(ji,jj  ) * ssmask(ji+1,jj  )   & 
     418               &            * ssmask(ji,jj+1) * ssmask(ji+1,jj+1) * MIN(1._wp,SUM(fmask(ji,jj,:))) 
     419         END DO 
     420      END DO 
     421      CALL lbc_lnk( umask_i, 'U', 1._wp )      ! Lateral boundary conditions 
     422      CALL lbc_lnk( vmask_i, 'V', 1._wp ) 
     423      CALL lbc_lnk( fmask_i, 'F', 1._wp ) 
     424 
     425      ! 3. Ocean/land mask at wu-, wv- and w points  
     426      !---------------------------------------------- 
     427      wmask (:,:,1) = tmask(:,:,1) ! ???????? 
     428      wumask(:,:,1) = umask(:,:,1) ! ???????? 
     429      wvmask(:,:,1) = vmask(:,:,1) ! ???????? 
     430      DO jk = 2, jpk 
     431         wmask (:,:,jk) = tmask(:,:,jk) * tmask(:,:,jk-1) 
     432         wumask(:,:,jk) = umask(:,:,jk) * umask(:,:,jk-1)    
     433         wvmask(:,:,jk) = vmask(:,:,jk) * vmask(:,:,jk-1) 
     434      END DO 
     435      ! 
     436      CALL wrk_dealloc( jpi, jpj, zmbk ) 
     437      ! 
     438      CALL iom_close( inum ) 
     439      ! 
     440   END SUBROUTINE dom_msk 
     441 
     442   SUBROUTINE zgr_bot_level 
     443      !!---------------------------------------------------------------------- 
     444      !!                    ***  ROUTINE zgr_bot_level  *** 
     445      !! 
     446      !! ** Purpose :   defines the vertical index of ocean bottom (mbk. arrays) 
     447      !! 
     448      !! ** Method  :   computes from mbathy with a minimum value of 1 over land 
     449      !! 
     450      !! ** Action  :   mbkt, mbku, mbkv :   vertical indices of the deeptest 
     451      !!                                     ocean level at t-, u- & v-points 
     452      !!                                     (min value = 1 over land) 
     453      !!---------------------------------------------------------------------- 
     454      ! 
     455      INTEGER ::   ji, jj   ! dummy loop indices 
     456      REAL(wp), POINTER, DIMENSION(:,:) :: zmbk 
     457      !!---------------------------------------------------------------------- 
     458 
     459      ! 
     460      IF(lwp) WRITE(numout,*) 
     461      IF(lwp) WRITE(numout,*) '    zgr_bot_level : ocean bottom k-index of T-, U-, V- and W-levels ' 
     462      IF(lwp) WRITE(numout,*) '    ~~~~~~~~~~~~~' 
     463      ! 
     464      CALL wrk_alloc( jpi, jpj, zmbk ) 
     465      ! 
     466      mbkt(:,:) = MAX( mbathy(:,:) , 1 )    ! bottom k-index of T-level (=1 over land) 
     467      mikt(:,:) = 1 ; miku(:,:) = 1; mikv(:,:) = 1; ! top k-index of T-level (=1 over open ocean; >1 beneath ice shelf) 
     468      !                                     ! bottom k-index of W-level = mbkt+1 
     469      DO jj = 1, jpjm1                      ! bottom k-index of u- (v-) level 
     470         DO ji = 1, jpim1 
     471            mbku(ji,jj) = MIN(  mbkt(ji+1,jj  ) , mbkt(ji,jj)  ) 
     472            mbkv(ji,jj) = MIN(  mbkt(ji  ,jj+1) , mbkt(ji,jj)  ) 
     473         END DO 
     474      END DO 
     475      ! converte into REAL to use lbc_lnk ; impose a min value of 1 as a zero can be set in lbclnk  
     476      zmbk(:,:) = REAL( mbku(:,:), wp )   ;   CALL lbc_lnk(zmbk,'U',1.)   ;   mbku  (:,:) = MAX( INT( zmbk(:,:) ), 1 ) 
     477      zmbk(:,:) = REAL( mbkv(:,:), wp )   ;   CALL lbc_lnk(zmbk,'V',1.)   ;   mbkv  (:,:) = MAX( INT( zmbk(:,:) ), 1 ) 
     478      ! 
     479      CALL wrk_dealloc( jpi, jpj, zmbk ) 
     480      ! 
     481   END SUBROUTINE zgr_bot_level 
     482 
     483   SUBROUTINE dom_hgr 
     484      !!---------------------------------------------------------------------- 
     485      !!                  ***  ROUTINE dom_hgr  *** 
     486      !!                    
     487      !! ** Purpose :  Read the NetCDF file(s) which contain(s) all the 
     488      !!      ocean horizontal mesh informations  
     489      !! 
     490      !! ** Method  :   Read in a file all the arrays generated in routines 
     491      !!                domhgr:   'mesh_hgr.nc' file 
     492      !!---------------------------------------------------------------------- 
     493      !! 
     494      INTEGER ::   ji, jj   ! dummy loop indices 
     495      INTEGER  ::  inum    ! local integers 
     496      !!---------------------------------------------------------------------- 
     497 
     498      IF(lwp) WRITE(numout,*) 
     499      IF(lwp) WRITE(numout,*) 'dom_grd_hgr : read NetCDF mesh and mask information file(s)' 
     500      IF(lwp) WRITE(numout,*) '~~~~~~~' 
     501 
     502      IF(lwp) WRITE(numout,*) '          one file in "mesh_mask.nc" ' 
     503      CALL iom_open( 'mesh_hgr', inum ) 
     504 
     505      !                                                         ! horizontal mesh (inum3) 
     506      CALL iom_get( inum, jpdom_data, 'glamt', glamt ) 
     507      CALL iom_get( inum, jpdom_data, 'glamu', glamu ) 
     508      CALL iom_get( inum, jpdom_data, 'glamv', glamv ) 
     509      CALL iom_get( inum, jpdom_data, 'glamf', glamf ) 
     510 
     511      CALL iom_get( inum, jpdom_data, 'gphit', gphit ) 
     512      CALL iom_get( inum, jpdom_data, 'gphiu', gphiu ) 
     513      CALL iom_get( inum, jpdom_data, 'gphiv', gphiv ) 
     514      CALL iom_get( inum, jpdom_data, 'gphif', gphif ) 
     515 
     516      CALL iom_get( inum, jpdom_data, 'e1t', e1t ) 
     517      CALL iom_get( inum, jpdom_data, 'e1u', e1u ) 
     518      CALL iom_get( inum, jpdom_data, 'e1v', e1v ) 
     519       
     520      CALL iom_get( inum, jpdom_data, 'e2t', e2t ) 
     521      CALL iom_get( inum, jpdom_data, 'e2u', e2u ) 
     522      CALL iom_get( inum, jpdom_data, 'e2v', e2v ) 
     523 
     524      CALL iom_get( inum, jpdom_data, 'ff', ff ) 
     525 
     526 
     527      ! Control printing : Grid informations (if not restart) 
     528      ! ---------------- 
     529 
     530      IF(lwp .AND. .NOT.ln_rstart ) THEN 
     531         WRITE(numout,*) 
     532         WRITE(numout,*) '          longitude and e1 scale factors' 
     533         WRITE(numout,*) '          ------------------------------' 
     534         WRITE(numout,9300) ( ji, glamt(ji,1), glamu(ji,1),   & 
     535            glamv(ji,1), glamf(ji,1),   & 
     536            e1t(ji,1), e1u(ji,1),   & 
     537            e1v(ji,1), ji = 1, jpi,10) 
     538 
     539         WRITE(numout,*) 
     540         WRITE(numout,*) '          latitude and e2 scale factors' 
     541         WRITE(numout,*) '          -----------------------------' 
     542         WRITE(numout,9300) ( jj, gphit(1,jj), gphiu(1,jj),   & 
     543            &                     gphiv(1,jj), gphif(1,jj),   & 
     544            &                     e2t  (1,jj), e2u  (1,jj),   & 
     545            &                     e2v  (1,jj), jj = 1, jpj, 10 ) 
     546      ENDIF 
     547 
     548      !                                     ! ============================ 
     549      !                                     !        close the files  
     550      !                                     ! ============================ 
     551      CALL iom_close( inum ) 
     552      ! 
     5539300     FORMAT( 1x, i4, f8.2,1x, f8.2,1x, f8.2,1x, f8.2, 1x,    & 
     554            f19.10, 1x, f19.10, 1x, f19.10 ) 
     555   END SUBROUTINE dom_hgr 
     556 
     557 
    323558   SUBROUTINE dom_zgr 
    324559      !!---------------------------------------------------------------------- 
    325560      !!                ***  ROUTINE dom_zgr  *** 
    326561      !!                    
    327       !! ** Purpose :  set the depth of model levels and the resulting  
    328       !!      vertical scale factors. 
     562      !! ** Purpose :  Read the NetCDF file(s) which contain(s) all the 
     563      !!      ocean horizontal mesh informations and/or set the depth of model levels  
     564      !!      and the resulting vertical scale factors. 
    329565      !! 
    330566      !! ** Method  : - reference 1D vertical coordinate (gdep._1d, e3._1d) 
     
    338574      !! ** Action  :   define gdep., e3., mbathy and bathy 
    339575      !!---------------------------------------------------------------------- 
    340       INTEGER ::   ioptio = 0   ! temporary integer 
    341       INTEGER ::   ios 
     576      INTEGER  ::  ioptio = 0   ! temporary integer 
     577      INTEGER  ::  inum, ios 
     578      INTEGER  ::  ji, jj, jk, ik 
     579      REAL(wp) ::  zrefdep 
    342580      !! 
    343581      NAMELIST/namzgr/ ln_zco, ln_zps, ln_sco, ln_isfcav 
     582      REAL(wp), POINTER, DIMENSION(:,:) :: zprt, zprw 
    344583      !!---------------------------------------------------------------------- 
    345584 
     
    372611      IF ( ioptio == 33 )   CALL ctl_stop( ' isf cavity with off line module not yet done    ' ) 
    373612 
    374    END SUBROUTINE dom_zgr 
    375  
    376    SUBROUTINE dom_ctl 
    377       !!---------------------------------------------------------------------- 
    378       !!                     ***  ROUTINE dom_ctl  *** 
    379       !! 
    380       !! ** Purpose :   Domain control. 
    381       !! 
    382       !! ** Method  :   compute and print extrema of masked scale factors 
    383       !! 
    384       !! History : 
    385       !!   8.5  !  02-08  (G. Madec)    Original code 
    386       !!---------------------------------------------------------------------- 
    387       !! * Local declarations 
    388       INTEGER ::   iimi1, ijmi1, iimi2, ijmi2, iima1, ijma1, iima2, ijma2 
    389       INTEGER, DIMENSION(2) ::   iloc      !  
    390       REAL(wp) ::   ze1min, ze1max, ze2min, ze2max 
    391       !!---------------------------------------------------------------------- 
    392  
    393       ! Extrema of the scale factors 
    394  
    395       IF(lwp)WRITE(numout,*) 
    396       IF(lwp)WRITE(numout,*) 'dom_ctl : extrema of the masked scale factors' 
    397       IF(lwp)WRITE(numout,*) '~~~~~~~' 
    398  
    399       IF (lk_mpp) THEN 
    400          CALL mpp_minloc( e1t(:,:), tmask(:,:,1), ze1min, iimi1,ijmi1 ) 
    401          CALL mpp_minloc( e2t(:,:), tmask(:,:,1), ze2min, iimi2,ijmi2 ) 
    402          CALL mpp_maxloc( e1t(:,:), tmask(:,:,1), ze1max, iima1,ijma1 ) 
    403          CALL mpp_maxloc( e2t(:,:), tmask(:,:,1), ze2max, iima2,ijma2 ) 
     613      IF(lwp) WRITE(numout,*) '          one file in "mesh_mask.nc" ' 
     614      CALL iom_open( 'mesh_zgr', inum ) 
     615 
     616      CALL iom_get( inum, jpdom_unknown, 'gdept_1d', gdept_1d ) ! depth 
     617      CALL iom_get( inum, jpdom_unknown, 'gdepw_1d', gdepw_1d ) 
     618      IF( ln_zco .OR. ln_zps ) THEN 
     619         CALL iom_get( inum, jpdom_unknown, 'e3t_1d'  , e3t_1d   )    ! reference scale factors 
     620         CALL iom_get( inum, jpdom_unknown, 'e3w_1d'  , e3w_1d   ) 
     621      ENDIF 
     622 
     623!!gm BUG in s-coordinate this does not work! 
     624      ! deepest/shallowest W level Above/Below ~10m 
     625      zrefdep = 10._wp - ( 0.1_wp * MINVAL(e3w_1d) )                 ! ref. depth with tolerance (10% of minimum layer thickness) 
     626      nlb10 = MINLOC( gdepw_1d, mask = gdepw_1d > zrefdep, dim = 1 ) ! shallowest W level Below ~10m 
     627      nla10 = nlb10 - 1                                              ! deepest    W level Above ~10m 
     628!!gm end bug 
     629 
     630      IF(lwp) THEN 
     631         WRITE(numout,*) 
     632         WRITE(numout,*) '              Reference z-coordinate depth and scale factors:' 
     633         WRITE(numout, "(9x,' level   gdept    gdepw     e3t      e3w  ')" ) 
     634         WRITE(numout, "(10x, i4, 4f9.2)" ) ( jk, gdept_1d(jk), gdepw_1d(jk), e3t_1d(jk), e3w_1d(jk), jk = 1, jpk ) 
     635      ENDIF 
     636 
     637      DO jk = 1, jpk 
     638         IF( e3w_1d  (jk) <= 0._wp .OR. e3t_1d  (jk) <= 0._wp )   CALL ctl_stop( ' e3w_1d or e3t_1d =< 0 ' ) 
     639         IF( gdepw_1d(jk) <  0._wp .OR. gdept_1d(jk) <  0._wp )   CALL ctl_stop( ' gdepw_1d or gdept_1d < 0 ' ) 
     640      END DO 
     641 
     642      IF( lk_vvl ) THEN 
     643          CALL iom_get( inum, jpdom_data, 'e3t_0', e3t_0(:,:,:) ) 
     644          CALL iom_get( inum, jpdom_data, 'e3u_0', e3u_0(:,:,:) ) 
     645          CALL iom_get( inum, jpdom_data, 'e3v_0', e3v_0(:,:,:) ) 
     646          CALL iom_get( inum, jpdom_data, 'e3w_0', e3w_0(:,:,:) ) 
     647          CALL iom_get( inum, jpdom_data, 'gdept_0', gdept_0(:,:,:) ) 
     648          CALL iom_get( inum, jpdom_data, 'gdepw_0', gdepw_0(:,:,:) ) 
     649          ht_0(:,:) = 0.0_wp                       ! Reference ocean depth at  T-points 
     650          DO jk = 1, jpk 
     651             ht_0(:,:) = ht_0(:,:) + e3t_0(:,:,jk) * tmask(:,:,jk) 
     652          END DO 
    404653      ELSE 
    405          ze1min = MINVAL( e1t(:,:), mask = tmask(:,:,1) == 1.e0 )     
    406          ze2min = MINVAL( e2t(:,:), mask = tmask(:,:,1) == 1.e0 )     
    407          ze1max = MAXVAL( e1t(:,:), mask = tmask(:,:,1) == 1.e0 )     
    408          ze2max = MAXVAL( e2t(:,:), mask = tmask(:,:,1) == 1.e0 )     
    409  
    410          iloc  = MINLOC( e1t(:,:), mask = tmask(:,:,1) == 1.e0 ) 
    411          iimi1 = iloc(1) + nimpp - 1 
    412          ijmi1 = iloc(2) + njmpp - 1 
    413          iloc  = MINLOC( e2t(:,:), mask = tmask(:,:,1) == 1.e0 ) 
    414          iimi2 = iloc(1) + nimpp - 1 
    415          ijmi2 = iloc(2) + njmpp - 1 
    416          iloc  = MAXLOC( e1t(:,:), mask = tmask(:,:,1) == 1.e0 ) 
    417          iima1 = iloc(1) + nimpp - 1 
    418          ijma1 = iloc(2) + njmpp - 1 
    419          iloc  = MAXLOC( e2t(:,:), mask = tmask(:,:,1) == 1.e0 ) 
    420          iima2 = iloc(1) + nimpp - 1 
    421          ijma2 = iloc(2) + njmpp - 1 
    422       ENDIF 
    423  
    424       IF(lwp) THEN 
    425          WRITE(numout,"(14x,'e1t maxi: ',1f10.2,' at i = ',i5,' j= ',i5)") ze1max, iima1, ijma1 
    426          WRITE(numout,"(14x,'e1t mini: ',1f10.2,' at i = ',i5,' j= ',i5)") ze1min, iimi1, ijmi1 
    427          WRITE(numout,"(14x,'e2t maxi: ',1f10.2,' at i = ',i5,' j= ',i5)") ze2max, iima2, ijma2 
    428          WRITE(numout,"(14x,'e2t mini: ',1f10.2,' at i = ',i5,' j= ',i5)") ze2min, iimi2, ijmi2 
    429       ENDIF 
    430  
    431    END SUBROUTINE dom_ctl 
    432  
    433    SUBROUTINE dom_grd 
    434       !!---------------------------------------------------------------------- 
    435       !!                  ***  ROUTINE dom_grd  *** 
    436       !!                    
    437       !! ** Purpose :  Read the NetCDF file(s) which contain(s) all the 
    438       !!      ocean domain informations (mesh and mask arrays). This (these) 
    439       !!      file(s) is (are) used for visualisation (SAXO software) and 
    440       !!      diagnostic computation. 
    441       !! 
    442       !! ** Method  :   Read in a file all the arrays generated in routines 
    443       !!      domhgr, domzgr, and dommsk. Note: the file contain depends on 
    444       !!      the vertical coord. used (z-coord, partial steps, s-coord) 
    445       !!                    nmsh = 1  :   'mesh_mask.nc' file 
    446       !!                         = 2  :   'mesh.nc' and mask.nc' files 
    447       !!                         = 3  :   'mesh_hgr.nc', 'mesh_zgr.nc' and 
    448       !!                                  'mask.nc' files 
    449       !!      For huge size domain, use option 2 or 3 depending on your  
    450       !!      vertical coordinate. 
    451       !! 
    452       !! ** input file :  
    453       !!      meshmask.nc  : domain size, horizontal grid-point position, 
    454       !!                     masks, depth and vertical scale factors 
    455       !!---------------------------------------------------------------------- 
    456       USE iom 
    457       !! 
    458       INTEGER  ::   ji, jj, jk   ! dummy loop indices 
    459       INTEGER  ::   ik, inum0 , inum1 , inum2 , inum3 , inum4   ! local integers 
    460       REAL(wp) ::   zrefdep         ! local real 
    461       REAL(wp), POINTER, DIMENSION(:,:) :: zmbk, zprt, zprw 
    462       !!---------------------------------------------------------------------- 
    463  
    464       IF(lwp) WRITE(numout,*) 
    465       IF(lwp) WRITE(numout,*) 'dom_rea : read NetCDF mesh and mask information file(s)' 
    466       IF(lwp) WRITE(numout,*) '~~~~~~~' 
    467  
    468       CALL wrk_alloc( jpi, jpj, zmbk, zprt, zprw ) 
    469  
    470       zmbk(:,:) = 0._wp 
    471  
    472       SELECT CASE (nmsh) 
    473          !                                     ! ============================ 
    474          CASE ( 1 )                            !  create 'mesh_mask.nc' file 
    475             !                                  ! ============================ 
    476  
    477             IF(lwp) WRITE(numout,*) '          one file in "mesh_mask.nc" ' 
    478             CALL iom_open( 'mesh_mask', inum0 ) 
    479  
    480             inum2 = inum0                                            ! put all the informations 
    481             inum3 = inum0                                            ! in unit inum0 
    482             inum4 = inum0 
    483  
    484             !                                  ! ============================ 
    485          CASE ( 2 )                            !  create 'mesh.nc' and  
    486             !                                  !         'mask.nc' files 
    487             !                                  ! ============================ 
    488  
    489             IF(lwp) WRITE(numout,*) '          two files in "mesh.nc" and "mask.nc" ' 
    490             CALL iom_open( 'mesh', inum1 ) 
    491             CALL iom_open( 'mask', inum2 ) 
    492  
    493             inum3 = inum1                                            ! put mesh informations  
    494             inum4 = inum1                                            ! in unit inum1  
    495  
    496             !                                  ! ============================ 
    497          CASE ( 3 )                            !  create 'mesh_hgr.nc' 
    498             !                                  !         'mesh_zgr.nc' and 
    499             !                                  !         'mask.nc'     files 
    500             !                                  ! ============================ 
    501  
    502             IF(lwp) WRITE(numout,*) '          three files in "mesh_hgr.nc" , "mesh_zgr.nc" and "mask.nc" ' 
    503             CALL iom_open( 'mesh_hgr', inum3 ) ! create 'mesh_hgr.nc' 
    504             CALL iom_open( 'mesh_zgr', inum4 ) ! create 'mesh_zgr.nc' 
    505             CALL iom_open( 'mask'    , inum2 ) ! create 'mask.nc' 
    506  
    507             !                                  ! =========================== 
    508          CASE DEFAULT                          ! return error  
    509             !                                  ! mesh has to be provided 
    510             !                                  ! =========================== 
    511             CALL ctl_stop( ' OFFLINE mode requires the input mesh mask(s). ',   & 
    512             &                                 ' Invalid nn_msh value in the namelist (0 is not allowed)' ) 
    513  
    514          END SELECT 
    515  
    516          !                                                         ! masks (inum2)  
    517          CALL iom_get( inum2, jpdom_data, 'tmask', tmask ) 
    518          CALL iom_get( inum2, jpdom_data, 'umask', umask ) 
    519          CALL iom_get( inum2, jpdom_data, 'vmask', vmask ) 
    520          CALL iom_get( inum2, jpdom_data, 'fmask', fmask ) 
    521  
    522          CALL lbc_lnk( tmask, 'T', 1._wp )    ! Lateral boundary conditions 
    523          CALL lbc_lnk( umask, 'U', 1._wp )       
    524          CALL lbc_lnk( vmask, 'V', 1._wp ) 
    525          CALL lbc_lnk( fmask, 'F', 1._wp ) 
    526  
    527 #if defined key_c1d 
    528          ! set umask and vmask equal tmask in 1D configuration 
    529          IF(lwp) WRITE(numout,*) 
    530          IF(lwp) WRITE(numout,*) '**********  1D configuration : set umask and vmask equal tmask ********' 
    531          IF(lwp) WRITE(numout,*) '**********                                                     ********' 
    532  
    533          umask(:,:,:) = tmask(:,:,:) 
    534          vmask(:,:,:) = tmask(:,:,:) 
    535 #endif 
    536  
    537 #if defined key_degrad 
    538          CALL iom_get( inum2, jpdom_data, 'facvolt', facvol ) 
    539 #endif 
    540  
    541          !                                                         ! horizontal mesh (inum3) 
    542          CALL iom_get( inum3, jpdom_data, 'glamt', glamt ) 
    543          CALL iom_get( inum3, jpdom_data, 'glamu', glamu ) 
    544          CALL iom_get( inum3, jpdom_data, 'glamv', glamv ) 
    545          CALL iom_get( inum3, jpdom_data, 'glamf', glamf ) 
    546  
    547          CALL iom_get( inum3, jpdom_data, 'gphit', gphit ) 
    548          CALL iom_get( inum3, jpdom_data, 'gphiu', gphiu ) 
    549          CALL iom_get( inum3, jpdom_data, 'gphiv', gphiv ) 
    550          CALL iom_get( inum3, jpdom_data, 'gphif', gphif ) 
    551  
    552          CALL iom_get( inum3, jpdom_data, 'e1t', e1t ) 
    553          CALL iom_get( inum3, jpdom_data, 'e1u', e1u ) 
    554          CALL iom_get( inum3, jpdom_data, 'e1v', e1v ) 
    555           
    556          CALL iom_get( inum3, jpdom_data, 'e2t', e2t ) 
    557          CALL iom_get( inum3, jpdom_data, 'e2u', e2u ) 
    558          CALL iom_get( inum3, jpdom_data, 'e2v', e2v ) 
    559  
    560          CALL iom_get( inum3, jpdom_data, 'ff', ff ) 
    561  
    562          CALL iom_get( inum4, jpdom_data, 'mbathy', zmbk )              ! number of ocean t-points 
    563          mbathy (:,:) = INT( zmbk(:,:) ) 
    564          misfdep(:,:) = 1                                               ! ice shelf case not yet done 
    565           
    566          CALL zgr_bot_level                                             ! mbk. arrays (deepest ocean t-, u- & v-points 
    567          ! 
    568654         IF( ln_sco ) THEN                                         ! s-coordinate 
    569             CALL iom_get( inum4, jpdom_data, 'hbatt', hbatt ) 
    570             CALL iom_get( inum4, jpdom_data, 'hbatu', hbatu ) 
    571             CALL iom_get( inum4, jpdom_data, 'hbatv', hbatv ) 
    572             CALL iom_get( inum4, jpdom_data, 'hbatf', hbatf ) 
     655            CALL iom_get( inum, jpdom_data, 'hbatt', hbatt ) 
     656            CALL iom_get( inum, jpdom_data, 'hbatu', hbatu ) 
     657            CALL iom_get( inum, jpdom_data, 'hbatv', hbatv ) 
     658            CALL iom_get( inum, jpdom_data, 'hbatf', hbatf ) 
    573659             
    574             CALL iom_get( inum4, jpdom_unknown, 'gsigt', gsigt ) ! scaling coef. 
    575             CALL iom_get( inum4, jpdom_unknown, 'gsigw', gsigw ) 
    576             CALL iom_get( inum4, jpdom_unknown, 'gsi3w', gsi3w )  
    577             CALL iom_get( inum4, jpdom_unknown, 'esigt', esigt ) 
    578             CALL iom_get( inum4, jpdom_unknown, 'esigw', esigw ) 
    579  
    580             CALL iom_get( inum4, jpdom_data, 'e3t_0', fse3t_n(:,:,:) ) ! scale factors 
    581             CALL iom_get( inum4, jpdom_data, 'e3u_0', fse3u_n(:,:,:) ) 
    582             CALL iom_get( inum4, jpdom_data, 'e3v_0', fse3v_n(:,:,:) ) 
    583             CALL iom_get( inum4, jpdom_data, 'e3w_0', fse3w_n(:,:,:) ) 
    584  
    585             CALL iom_get( inum4, jpdom_unknown, 'gdept_1d', gdept_1d ) ! depth 
    586             CALL iom_get( inum4, jpdom_unknown, 'gdepw_1d', gdepw_1d ) 
     660            CALL iom_get( inum, jpdom_unknown, 'gsigt', gsigt ) ! scaling coef. 
     661            CALL iom_get( inum, jpdom_unknown, 'gsigw', gsigw ) 
     662            CALL iom_get( inum, jpdom_unknown, 'gsi3w', gsi3w )  
     663            CALL iom_get( inum, jpdom_unknown, 'esigt', esigt ) 
     664            CALL iom_get( inum, jpdom_unknown, 'esigw', esigw ) 
     665 
     666            CALL iom_get( inum, jpdom_data, 'e3t_0', fse3t_n(:,:,:) ) ! scale factors 
     667            CALL iom_get( inum, jpdom_data, 'e3u_0', fse3u_n(:,:,:) ) 
     668            CALL iom_get( inum, jpdom_data, 'e3v_0', fse3v_n(:,:,:) ) 
     669            CALL iom_get( inum, jpdom_data, 'e3w_0', fse3w_n(:,:,:) ) 
    587670         ENDIF 
    588671 
    589672  
    590673         IF( ln_zps ) THEN                                           ! z-coordinate - partial steps 
    591             CALL iom_get( inum4, jpdom_unknown, 'gdept_1d', gdept_1d )  ! reference depth 
    592             CALL iom_get( inum4, jpdom_unknown, 'gdepw_1d', gdepw_1d ) 
    593             CALL iom_get( inum4, jpdom_unknown, 'e3t_1d'  , e3t_1d   )    ! reference scale factors 
    594             CALL iom_get( inum4, jpdom_unknown, 'e3w_1d'  , e3w_1d   ) 
    595674            ! 
    596             IF( nmsh <= 6 ) THEN                                        ! 3D vertical scale factors 
    597                CALL iom_get( inum4, jpdom_data, 'e3t_0', fse3t_n(:,:,:) ) 
    598                CALL iom_get( inum4, jpdom_data, 'e3u_0', fse3u_n(:,:,:) ) 
    599                CALL iom_get( inum4, jpdom_data, 'e3v_0', fse3v_n(:,:,:) ) 
    600                CALL iom_get( inum4, jpdom_data, 'e3w_0', fse3w_n(:,:,:) ) 
     675            IF( iom_varid( inum, 'e3t_0', ldstop = .FALSE. ) > 0 ) THEN 
     676               CALL iom_get( inum, jpdom_data, 'e3t_0', fse3t_n(:,:,:) ) 
     677               CALL iom_get( inum, jpdom_data, 'e3u_0', fse3u_n(:,:,:) ) 
     678               CALL iom_get( inum, jpdom_data, 'e3v_0', fse3v_n(:,:,:) ) 
     679               CALL iom_get( inum, jpdom_data, 'e3w_0', fse3w_n(:,:,:) ) 
    601680            ELSE                                                        ! 2D bottom scale factors 
    602                CALL iom_get( inum4, jpdom_data, 'e3t_ps', e3tp ) 
    603                CALL iom_get( inum4, jpdom_data, 'e3w_ps', e3wp ) 
     681               CALL iom_get( inum, jpdom_data, 'e3t_ps', e3tp ) 
     682               CALL iom_get( inum, jpdom_data, 'e3w_ps', e3wp ) 
    604683               !                                                        ! deduces the 3D scale factors 
    605684               DO jk = 1, jpk 
     
    633712            END IF 
    634713 
    635             IF( iom_varid( inum4, 'gdept_0', ldstop = .FALSE. ) > 0 ) THEN   ! 3D depth of t- and w-level 
    636                CALL iom_get( inum4, jpdom_data, 'gdept_0', fsdept_n(:,:,:) ) 
    637                CALL iom_get( inum4, jpdom_data, 'gdepw_0', fsdepw_n(:,:,:) ) 
     714            IF( iom_varid( inum, 'gdept_0', ldstop = .FALSE. ) > 0 ) THEN   ! 3D depth of t- and w-level 
     715               CALL iom_get( inum, jpdom_data, 'gdept_0', fsdept_n(:,:,:) ) 
     716               CALL iom_get( inum, jpdom_data, 'gdepw_0', fsdepw_n(:,:,:) ) 
    638717            ELSE                                                           ! 2D bottom depth 
    639                CALL iom_get( inum4, jpdom_data, 'hdept', zprt ) 
    640                CALL iom_get( inum4, jpdom_data, 'hdepw', zprw ) 
     718               CALL wrk_alloc( jpi, jpj, zprt, zprw ) 
     719               ! 
     720               CALL iom_get( inum, jpdom_data, 'hdept', zprt ) 
     721               CALL iom_get( inum, jpdom_data, 'hdepw', zprw ) 
    641722               ! 
    642723               DO jk = 1, jpk                                              ! deduces the 3D depth 
     
    654735                  END DO 
    655736               END DO 
     737               CALL wrk_dealloc( jpi, jpj, zprt, zprw ) 
    656738            ENDIF 
    657739            ! 
     
    659741 
    660742         IF( ln_zco ) THEN           ! Vertical coordinates and scales factors 
    661             CALL iom_get( inum4, jpdom_unknown, 'gdept_1d', gdept_1d ) ! depth 
    662             CALL iom_get( inum4, jpdom_unknown, 'gdepw_1d', gdepw_1d ) 
    663             CALL iom_get( inum4, jpdom_unknown, 'e3t_1d'  , e3t_1d   ) 
    664             CALL iom_get( inum4, jpdom_unknown, 'e3w_1d'  , e3w_1d   ) 
    665743            DO jk = 1, jpk 
    666744               fse3t_n(:,:,jk) = e3t_1d(jk)                              ! set to the ref. factors 
     
    672750            END DO 
    673751         ENDIF 
    674  
    675 !!gm BUG in s-coordinate this does not work! 
    676       ! deepest/shallowest W level Above/Below ~10m 
    677       zrefdep = 10._wp - ( 0.1_wp * MINVAL(e3w_1d) )                 ! ref. depth with tolerance (10% of minimum layer thickness) 
    678       nlb10 = MINLOC( gdepw_1d, mask = gdepw_1d > zrefdep, dim = 1 ) ! shallowest W level Below ~10m 
    679       nla10 = nlb10 - 1                                              ! deepest    W level Above ~10m 
    680 !!gm end bug 
    681  
    682       ! Control printing : Grid informations (if not restart) 
    683       ! ---------------- 
    684  
    685       IF(lwp .AND. .NOT.ln_rstart ) THEN 
    686          WRITE(numout,*) 
    687          WRITE(numout,*) '          longitude and e1 scale factors' 
    688          WRITE(numout,*) '          ------------------------------' 
    689          WRITE(numout,9300) ( ji, glamt(ji,1), glamu(ji,1),   & 
    690             glamv(ji,1), glamf(ji,1),   & 
    691             e1t(ji,1), e1u(ji,1),   & 
    692             e1v(ji,1), ji = 1, jpi,10) 
    693 9300     FORMAT( 1x, i4, f8.2,1x, f8.2,1x, f8.2,1x, f8.2, 1x,    & 
    694             f19.10, 1x, f19.10, 1x, f19.10 ) 
    695  
    696          WRITE(numout,*) 
    697          WRITE(numout,*) '          latitude and e2 scale factors' 
    698          WRITE(numout,*) '          -----------------------------' 
    699          WRITE(numout,9300) ( jj, gphit(1,jj), gphiu(1,jj),   & 
    700             &                     gphiv(1,jj), gphif(1,jj),   & 
    701             &                     e2t  (1,jj), e2u  (1,jj),   & 
    702             &                     e2v  (1,jj), jj = 1, jpj, 10 ) 
    703       ENDIF 
    704  
    705  
    706       IF( nprint == 1 .AND. lwp ) THEN 
    707          WRITE(numout,*) '          e1u e2u ' 
    708          CALL prihre( e1u,jpi,jpj,jpi-5,jpi,1,jpj-5,jpj,1,0.,numout ) 
    709          CALL prihre( e2u,jpi,jpj,jpi-5,jpi,1,jpj-5,jpj,1,0.,numout ) 
    710          WRITE(numout,*) '          e1v e2v  ' 
    711          CALL prihre( e1v,jpi,jpj,jpi-5,jpi,1,jpj-5,jpj,1,0.,numout ) 
    712          CALL prihre( e2v,jpi,jpj,jpi-5,jpi,1,jpj-5,jpj,1,0.,numout ) 
    713       ENDIF 
    714  
    715       IF(lwp) THEN 
    716          WRITE(numout,*) 
    717          WRITE(numout,*) '              Reference z-coordinate depth and scale factors:' 
    718          WRITE(numout, "(9x,' level   gdept    gdepw     e3t      e3w  ')" ) 
    719          WRITE(numout, "(10x, i4, 4f9.2)" ) ( jk, gdept_1d(jk), gdepw_1d(jk), e3t_1d(jk), e3w_1d(jk), jk = 1, jpk ) 
    720       ENDIF 
    721  
    722       DO jk = 1, jpk 
    723          IF( e3w_1d  (jk) <= 0._wp .OR. e3t_1d  (jk) <= 0._wp )   CALL ctl_stop( ' e3w_1d or e3t_1d =< 0 ' ) 
    724          IF( gdepw_1d(jk) <  0._wp .OR. gdept_1d(jk) <  0._wp )   CALL ctl_stop( ' gdepw_1d or gdept_1d < 0 ' ) 
    725       END DO 
     752         ! 
     753      ENDIF 
    726754      !                                     ! ============================ 
    727755      !                                     !        close the files  
    728756      !                                     ! ============================ 
    729       SELECT CASE ( nmsh ) 
    730          CASE ( 1 )                 
    731             CALL iom_close( inum0 ) 
    732          CASE ( 2 ) 
    733             CALL iom_close( inum1 ) 
    734             CALL iom_close( inum2 ) 
    735          CASE ( 3 ) 
    736             CALL iom_close( inum2 ) 
    737             CALL iom_close( inum3 ) 
    738             CALL iom_close( inum4 ) 
    739       END SELECT 
    740       ! 
    741       CALL wrk_dealloc( jpi, jpj, zmbk, zprt, zprw ) 
    742       ! 
    743    END SUBROUTINE dom_grd 
    744  
    745  
    746    SUBROUTINE zgr_bot_level 
    747       !!---------------------------------------------------------------------- 
    748       !!                    ***  ROUTINE zgr_bot_level  *** 
    749       !! 
    750       !! ** Purpose :   defines the vertical index of ocean bottom (mbk. arrays) 
    751       !! 
    752       !! ** Method  :   computes from mbathy with a minimum value of 1 over land 
    753       !! 
    754       !! ** Action  :   mbkt, mbku, mbkv :   vertical indices of the deeptest 
    755       !!                                     ocean level at t-, u- & v-points 
    756       !!                                     (min value = 1 over land) 
    757       !!---------------------------------------------------------------------- 
    758       ! 
    759       INTEGER ::   ji, jj   ! dummy loop indices 
    760       REAL(wp), POINTER, DIMENSION(:,:) :: zmbk 
    761       !!---------------------------------------------------------------------- 
    762  
    763       ! 
    764       IF(lwp) WRITE(numout,*) 
    765       IF(lwp) WRITE(numout,*) '    zgr_bot_level : ocean bottom k-index of T-, U-, V- and W-levels ' 
    766       IF(lwp) WRITE(numout,*) '    ~~~~~~~~~~~~~' 
    767       ! 
    768       CALL wrk_alloc( jpi, jpj, zmbk ) 
    769       ! 
    770       mbkt(:,:) = MAX( mbathy(:,:) , 1 )    ! bottom k-index of T-level (=1 over land) 
    771       mikt(:,:) = 1 ; miku(:,:) = 1; mikv(:,:) = 1; ! top k-index of T-level (=1 over open ocean; >1 beneath ice shelf) 
    772       !                                     ! bottom k-index of W-level = mbkt+1 
    773       DO jj = 1, jpjm1                      ! bottom k-index of u- (v-) level 
    774          DO ji = 1, jpim1 
    775             mbku(ji,jj) = MIN(  mbkt(ji+1,jj  ) , mbkt(ji,jj)  ) 
    776             mbkv(ji,jj) = MIN(  mbkt(ji  ,jj+1) , mbkt(ji,jj)  ) 
    777          END DO 
    778       END DO 
    779       ! converte into REAL to use lbc_lnk ; impose a min value of 1 as a zero can be set in lbclnk  
    780       zmbk(:,:) = REAL( mbku(:,:), wp )   ;   CALL lbc_lnk(zmbk,'U',1.)   ;   mbku  (:,:) = MAX( INT( zmbk(:,:) ), 1 ) 
    781       zmbk(:,:) = REAL( mbkv(:,:), wp )   ;   CALL lbc_lnk(zmbk,'V',1.)   ;   mbkv  (:,:) = MAX( INT( zmbk(:,:) ), 1 ) 
    782       ! 
    783       CALL wrk_dealloc( jpi, jpj, zmbk ) 
    784       ! 
    785    END SUBROUTINE zgr_bot_level 
    786  
    787    SUBROUTINE dom_msk 
    788       !!--------------------------------------------------------------------- 
    789       !!                 ***  ROUTINE dom_msk  *** 
    790       !! 
    791       !! ** Purpose :   Off-line case: defines the interior domain T-mask. 
    792       !! 
    793       !! ** Method  :   The interior ocean/land mask is computed from tmask 
    794       !!              setting to zero the duplicated row and lines due to 
    795       !!              MPP exchange halos, est-west cyclic and north fold 
    796       !!              boundary conditions. 
    797       !! 
    798       !! ** Action :   tmask_i  : interiorland/ocean mask at t-point 
    799       !!               tpol     : ??? 
    800       !!---------------------------------------------------------------------- 
    801       ! 
    802       INTEGER  ::   ji, jj, jk                   ! dummy loop indices 
    803       INTEGER  ::   iif, iil, ijf, ijl       ! local integers 
    804       INTEGER, POINTER, DIMENSION(:,:) ::  imsk  
    805       ! 
    806       !!--------------------------------------------------------------------- 
    807        
    808       CALL wrk_alloc( jpi, jpj, imsk ) 
    809       ! 
    810       ! Interior domain mask (used for global sum) 
    811       ! -------------------- 
    812       ssmask(:,:)  = tmask(:,:,1) 
    813       tmask_i(:,:) = tmask(:,:,1) 
    814       iif = jpreci                        ! thickness of exchange halos in i-axis 
    815       iil = nlci - jpreci + 1 
    816       ijf = jprecj                        ! thickness of exchange halos in j-axis 
    817       ijl = nlcj - jprecj + 1 
    818       ! 
    819       tmask_i( 1 :iif,   :   ) = 0._wp    ! first columns 
    820       tmask_i(iil:jpi,   :   ) = 0._wp    ! last  columns (including mpp extra columns) 
    821       tmask_i(   :   , 1 :ijf) = 0._wp    ! first rows 
    822       tmask_i(   :   ,ijl:jpj) = 0._wp    ! last  rows (including mpp extra rows) 
    823       ! 
    824       !                                   ! north fold mask 
    825       tpol(1:jpiglo) = 1._wp 
    826       !                                 
    827       IF( jperio == 3 .OR. jperio == 4 )   tpol(jpiglo/2+1:jpiglo) = 0._wp    ! T-point pivot 
    828       IF( jperio == 5 .OR. jperio == 6 )   tpol(     1    :jpiglo) = 0._wp    ! F-point pivot 
    829       IF( jperio == 3 .OR. jperio == 4 ) THEN      ! T-point pivot: only half of the nlcj-1 row 
    830          IF( mjg(ijl-1) == jpjglo-1 ) THEN 
    831             DO ji = iif+1, iil-1 
    832                tmask_i(ji,ijl-1) = tmask_i(ji,ijl-1) * tpol(mig(ji)) 
    833             END DO 
    834          ENDIF 
    835       ENDIF  
    836       ! 
    837       ! (ISF) MIN(1,SUM(umask)) is here to check if you have effectively at 
    838       ! least 1 wet u point 
    839       DO jj = 1, jpjm1 
    840          DO ji = 1, fs_jpim1   ! vector loop 
    841             umask_i(ji,jj)  = ssmask(ji,jj) * ssmask(ji+1,jj  )  * MIN(1._wp,SUM(umask(ji,jj,:))) 
    842             vmask_i(ji,jj)  = ssmask(ji,jj) * ssmask(ji  ,jj+1)  * MIN(1._wp,SUM(vmask(ji,jj,:))) 
    843          END DO 
    844          DO ji = 1, jpim1      ! NO vector opt. 
    845             fmask_i(ji,jj) =  ssmask(ji,jj  ) * ssmask(ji+1,jj  )   & 
    846                &            * ssmask(ji,jj+1) * ssmask(ji+1,jj+1) * MIN(1._wp,SUM(fmask(ji,jj,:))) 
    847          END DO 
    848       END DO 
    849       CALL lbc_lnk( umask_i, 'U', 1._wp )      ! Lateral boundary conditions 
    850       CALL lbc_lnk( vmask_i, 'V', 1._wp ) 
    851       CALL lbc_lnk( fmask_i, 'F', 1._wp ) 
    852  
    853       ! 3. Ocean/land mask at wu-, wv- and w points  
    854       !---------------------------------------------- 
    855       wmask (:,:,1) = tmask(:,:,1) ! ???????? 
    856       wumask(:,:,1) = umask(:,:,1) ! ???????? 
    857       wvmask(:,:,1) = vmask(:,:,1) ! ???????? 
    858       DO jk=2,jpk 
    859          wmask (:,:,jk)=tmask(:,:,jk) * tmask(:,:,jk-1) 
    860          wumask(:,:,jk)=umask(:,:,jk) * umask(:,:,jk-1)    
    861          wvmask(:,:,jk)=vmask(:,:,jk) * vmask(:,:,jk-1) 
    862       END DO 
    863       ! 
    864       IF( nprint == 1 .AND. lwp ) THEN    ! Control print 
    865          imsk(:,:) = INT( tmask_i(:,:) ) 
    866          WRITE(numout,*) ' tmask_i : ' 
    867          CALL prihin( imsk(:,:), jpi, jpj, 1, jpi, 1, 1, jpj, 1, 1, numout) 
    868          WRITE (numout,*) 
    869          WRITE (numout,*) ' dommsk: tmask for each level' 
    870          WRITE (numout,*) ' ----------------------------' 
    871          DO jk = 1, jpk 
    872             imsk(:,:) = INT( tmask(:,:,jk) ) 
    873             WRITE(numout,*) 
    874             WRITE(numout,*) ' level = ',jk 
    875             CALL prihin( imsk(:,:), jpi, jpj, 1, jpi, 1, 1, jpj, 1, 1, numout) 
    876          END DO 
    877       ENDIF 
    878       ! 
    879       CALL wrk_dealloc( jpi, jpj, imsk ) 
    880       ! 
    881    END SUBROUTINE dom_msk 
     757      CALL iom_close( inum ) 
     758      ! 
     759      ! 
     760   END SUBROUTINE dom_zgr 
     761 
     762   SUBROUTINE dom_ctl 
     763      !!---------------------------------------------------------------------- 
     764      !!                     ***  ROUTINE dom_ctl  *** 
     765      !! 
     766      !! ** Purpose :   Domain control. 
     767      !! 
     768      !! ** Method  :   compute and print extrema of masked scale factors 
     769      !! 
     770      !! History : 
     771      !!   8.5  !  02-08  (G. Madec)    Original code 
     772      !!---------------------------------------------------------------------- 
     773      !! * Local declarations 
     774      INTEGER ::   iimi1, ijmi1, iimi2, ijmi2, iima1, ijma1, iima2, ijma2 
     775      INTEGER, DIMENSION(2) ::   iloc      !  
     776      REAL(wp) ::   ze1min, ze1max, ze2min, ze2max 
     777      !!---------------------------------------------------------------------- 
     778 
     779      ! Extrema of the scale factors 
     780 
     781      IF(lwp)WRITE(numout,*) 
     782      IF(lwp)WRITE(numout,*) 'dom_ctl : extrema of the masked scale factors' 
     783      IF(lwp)WRITE(numout,*) '~~~~~~~' 
     784 
     785      IF (lk_mpp) THEN 
     786         CALL mpp_minloc( e1t(:,:), tmask(:,:,1), ze1min, iimi1,ijmi1 ) 
     787         CALL mpp_minloc( e2t(:,:), tmask(:,:,1), ze2min, iimi2,ijmi2 ) 
     788         CALL mpp_maxloc( e1t(:,:), tmask(:,:,1), ze1max, iima1,ijma1 ) 
     789         CALL mpp_maxloc( e2t(:,:), tmask(:,:,1), ze2max, iima2,ijma2 ) 
     790      ELSE 
     791         ze1min = MINVAL( e1t(:,:), mask = tmask(:,:,1) == 1.e0 )     
     792         ze2min = MINVAL( e2t(:,:), mask = tmask(:,:,1) == 1.e0 )     
     793         ze1max = MAXVAL( e1t(:,:), mask = tmask(:,:,1) == 1.e0 )     
     794         ze2max = MAXVAL( e2t(:,:), mask = tmask(:,:,1) == 1.e0 )     
     795 
     796         iloc  = MINLOC( e1t(:,:), mask = tmask(:,:,1) == 1.e0 ) 
     797         iimi1 = iloc(1) + nimpp - 1 
     798         ijmi1 = iloc(2) + njmpp - 1 
     799         iloc  = MINLOC( e2t(:,:), mask = tmask(:,:,1) == 1.e0 ) 
     800         iimi2 = iloc(1) + nimpp - 1 
     801         ijmi2 = iloc(2) + njmpp - 1 
     802         iloc  = MAXLOC( e1t(:,:), mask = tmask(:,:,1) == 1.e0 ) 
     803         iima1 = iloc(1) + nimpp - 1 
     804         ijma1 = iloc(2) + njmpp - 1 
     805         iloc  = MAXLOC( e2t(:,:), mask = tmask(:,:,1) == 1.e0 ) 
     806         iima2 = iloc(1) + nimpp - 1 
     807         ijma2 = iloc(2) + njmpp - 1 
     808      ENDIF 
     809 
     810      IF(lwp) THEN 
     811         WRITE(numout,"(14x,'e1t maxi: ',1f10.2,' at i = ',i5,' j= ',i5)") ze1max, iima1, ijma1 
     812         WRITE(numout,"(14x,'e1t mini: ',1f10.2,' at i = ',i5,' j= ',i5)") ze1min, iimi1, ijmi1 
     813         WRITE(numout,"(14x,'e2t maxi: ',1f10.2,' at i = ',i5,' j= ',i5)") ze2max, iima2, ijma2 
     814         WRITE(numout,"(14x,'e2t mini: ',1f10.2,' at i = ',i5,' j= ',i5)") ze2min, iimi2, ijmi2 
     815      ENDIF 
     816 
     817   END SUBROUTINE dom_ctl 
    882818 
    883819   !!====================================================================== 
Note: See TracChangeset for help on using the changeset viewer.