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 2528 for trunk/NEMOGCM/NEMO/OPA_SRC/DOM – NEMO

Ignore:
Timestamp:
2010-12-27T18:33:53+01:00 (13 years ago)
Author:
rblod
Message:

Update NEMOGCM from branch nemo_v3_3_beta

Location:
trunk/NEMOGCM/NEMO/OPA_SRC/DOM
Files:
14 edited
1 copied

Legend:

Unmodified
Added
Removed

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

    • Property svn:eol-style deleted
    r1601 r2528  
    4444#  include "vectopt_loop_substitute.h90" 
    4545   !!---------------------------------------------------------------------- 
    46    !!  OPA 9.0 , LOCEAN-IPSL (2006)  
     46   !! NEMO/OPA 3.3 , NEMO Consortium (2010) 
    4747   !! $Id$ 
    48    !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 
     48   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 
    4949   !!---------------------------------------------------------------------- 
    5050 
    … …  
    180180      INTEGER                     ::   ji, jj, jc, jn   ! dummy loop indices 
    181181      REAL(wp)                    ::   zze2 
    182       REAL(wp), DIMENSION (jpncs) ::   zemp 
     182      REAL(wp), DIMENSION (jpncs) ::   zfwf  
     183  
    183184      !!---------------------------------------------------------------------- 
    184185      ! 
    … …  
    216217      !                                                   !--------------------! 
    217218      !                                                   !  update emp, emps  ! 
    218       zemp = 0.e0                                         !--------------------! 
     219      zfwf = 0.e0                                         !--------------------! 
    219220      DO jc = 1, jpncs 
    220221         DO jj = ncsj1(jc), ncsj2(jc) 
    221222            DO ji = ncsi1(jc), ncsi2(jc) 
    222                zemp(jc) = zemp(jc) + e1t(ji,jj) * e2t(ji,jj) * emp(ji,jj) * tmask_i(ji,jj) 
     223               zfwf(jc) = zfwf(jc) + e1t(ji,jj) * e2t(ji,jj) * ( emp(ji,jj)-rnf(ji,jj) ) * tmask_i(ji,jj)  
    223224            END DO   
    224225         END DO  
    225226      END DO 
    226       IF( lk_mpp )   CALL mpp_sum ( zemp(:) , jpncs )       ! mpp: sum over all the global domain 
     227      IF( lk_mpp )   CALL mpp_sum ( zfwf(:) , jpncs )       ! mpp: sum over all the global domain 
    227228 
    228229      IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN      ! Black Sea case for ORCA_R2 configuration 
    229          zze2    = ( zemp(3) + zemp(4) ) / 2. 
    230          zemp(3) = zze2 
    231          zemp(4) = zze2 
     230         zze2    = ( zfwf(3) + zfwf(4) ) / 2. 
     231         zfwf(3) = zze2 
     232         zfwf(4) = zze2 
    232233      ENDIF 
    233234 
    … …  
    236237         IF( ncstt(jc) == 0 ) THEN  
    237238            ! water/evap excess is shared by all open ocean 
    238             emp (:,:) = emp (:,:) + zemp(jc) / surf(jpncs+1) 
    239             emps(:,:) = emps(:,:) + zemp(jc) / surf(jpncs+1) 
     239            emp (:,:) = emp (:,:) + zfwf(jc) / surf(jpncs+1) 
     240            emps(:,:) = emps(:,:) + zfwf(jc) / surf(jpncs+1) 
    240241         ELSEIF( ncstt(jc) == 1 ) THEN  
    241242            ! Excess water in open sea, at outflow location, excess evap shared 
    242             IF ( zemp(jc) <= 0.e0 ) THEN  
     243            IF ( zfwf(jc) <= 0.e0 ) THEN  
    243244                DO jn = 1, ncsnr(jc) 
    244245                  ji = mi0(ncsir(jc,jn)) 
    … …  
    246247                  IF (      ji > 1 .AND. ji < jpi   & 
    247248                      .AND. jj > 1 .AND. jj < jpj ) THEN  
    248                       emp (ji,jj) = emp (ji,jj) + zemp(jc) /   & 
     249                      emp (ji,jj) = emp (ji,jj) + zfwf(jc) /   & 
    249250                         (FLOAT(ncsnr(jc)) * e1t(ji,jj) * e2t(ji,jj)) 
    250                       emps(ji,jj) = emps(ji,jj) + zemp(jc) /   & 
     251                      emps(ji,jj) = emps(ji,jj) + zfwf(jc) /   & 
    251252                          (FLOAT(ncsnr(jc)) * e1t(ji,jj) * e2t(ji,jj)) 
    252253                  END IF  
    253254                END DO  
    254255            ELSE  
    255                 emp (:,:) = emp (:,:) + zemp(jc) / surf(jpncs+1) 
    256                 emps(:,:) = emps(:,:) + zemp(jc) / surf(jpncs+1) 
     256                emp (:,:) = emp (:,:) + zfwf(jc) / surf(jpncs+1) 
     257                emps(:,:) = emps(:,:) + zfwf(jc) / surf(jpncs+1) 
    257258            ENDIF 
    258259         ELSEIF( ncstt(jc) == 2 ) THEN  
    … …  
    263264                  ji = mi0(ncsir(jc,jn)) 
    264265                  jj = mj0(ncsjr(jc,jn)) ! Location of outflow in open ocean 
    265                   emp (ji,jj) = emp (ji,jj) + zemp(jc)   & 
     266                  emp (ji,jj) = emp (ji,jj) + zfwf(jc)   & 
    266267                      / (FLOAT(ncsnr(jc)) *  e1t(ji,jj) * e2t(ji,jj) ) 
    267                   emps(ji,jj) = emps(ji,jj) + zemp(jc)   & 
     268                  emps(ji,jj) = emps(ji,jj) + zfwf(jc)   & 
    268269                      / (FLOAT(ncsnr(jc)) *  e1t(ji,jj) * e2t(ji,jj) ) 
    269270                END DO  
    … …  
    273274         DO jj = ncsj1(jc), ncsj2(jc) 
    274275            DO ji = ncsi1(jc), ncsi2(jc) 
    275                emp (ji,jj) = emp (ji,jj) - zemp(jc) / surf(jc) 
    276                emps(ji,jj) = emps(ji,jj) - zemp(jc) / surf(jc) 
     276               emp (ji,jj) = emp (ji,jj) - zfwf(jc) / surf(jc) 
     277               emps(ji,jj) = emps(ji,jj) - zfwf(jc) / surf(jc) 
    277278            END DO   
    278279         END DO  

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

    • Property svn:eol-style deleted
    r2191 r2528  
    3333   USE prtctl          ! Print control 
    3434   USE restart         !  
     35   USE trc_oce, ONLY : lk_offline ! offline flag 
    3536 
    3637   IMPLICIT NONE 
    … …  
    4344 
    4445   !!---------------------------------------------------------------------- 
    45    !! NEMO/OPA 3.2 , LOCEAN-IPSL (2009)  
     46   !! NEMO/OPA 3.3 , NEMO Consortium (2010) 
    4647   !! $Id$ 
    47    !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 
     48   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 
    4849   !!---------------------------------------------------------------------- 
    4950 
    … …  
    6768      !!              - nmonth_len, nyear_len, nmonth_half, nmonth_end through day_mth 
    6869      !!---------------------------------------------------------------------- 
     70      INTEGER :: inbday, idweek 
     71      REAL(wp) :: zjul 
     72      !!---------------------------------------------------------------------- 
    6973 
    7074      ! all calendar staff is based on the fact that MOD( rday, rdttra(1) ) == 0 
    … …  
    7781      ndt05   = NINT(0.5 * rdttra(1)) 
    7882 
    79       CALL day_rst( nit000, 'READ' )  
     83      IF( .NOT. lk_offline ) CALL day_rst( nit000, 'READ' )  
    8084 
    8185      ! set the calandar from ndastp (read in restart file and namelist) 
    … …  
    105109      ! day since january 1st 
    106110      nday_year = nday + SUM( nmonth_len(1:nmonth - 1) ) 
    107        
    108       ! number of seconds since the beginning of current year/month at the middle of the time-step 
     111 
     112      !compute number of days between last monday and today       
     113      CALL ymds2ju( 1900, 01, 01, 0.0, zjul )  ! compute julian day value of 01.01.1900 (our reference that was a Monday) 
     114      inbday = NINT(fjulday - zjul)            ! compute nb day between  01.01.1900 and current day   
     115      idweek = MOD(inbday, 7)                  ! compute nb day between last monday and current day   
     116 
     117      ! number of seconds since the beginning of current year/month/week/day at the middle of the time-step 
    109118      nsec_year  = nday_year * nsecd - ndt05   ! 1 time step before the middle of the first time step 
    110119      nsec_month = nday      * nsecd - ndt05   ! because day will be called at the beginning of step 
     120      nsec_week  = idweek    * nsecd - ndt05 
    111121      nsec_day   =             nsecd - ndt05 
    112122 
    113123      ! control print 
    114124      IF(lwp) WRITE(numout,'(a,i6,a,i2,a,i2,a,i6)')' ==============>> 1/2 time step before the start of the run DATE Y/M/D = ',   & 
    115            &                   nyear, '/', nmonth, '/', nday, '  nsec_day:', nsec_day 
     125           &                   nyear, '/', nmonth, '/', nday, '  nsec_day:', nsec_day, '  nsec_week:', nsec_week 
    116126 
    117127      ! Up to now, calendar parameters are related to the end of previous run (nit000-1) 
    … …  
    200210      nsec_year  = nsec_year  + ndt  
    201211      nsec_month = nsec_month + ndt                  
     212      nsec_week  = nsec_week  + ndt 
    202213      nsec_day   = nsec_day   + ndt                 
    203214      adatrj  = adatrj  + rdttra(1) / rday 
    … …  
    206217      IF( ABS(adatrj  - REAL(NINT(adatrj ),wp)) < zprec )   adatrj  = REAL(NINT(adatrj ),wp)   ! avoid truncation error 
    207218       
    208       IF( nsec_day > nsecd ) THEN                        ! NEW day 
     219      IF( nsec_day > nsecd ) THEN                       ! New day 
    209220         ! 
    210221         nday      = nday + 1 
    … …  
    212223         nsec_day  = ndt05 
    213224         ! 
    214          IF( nday == nmonth_len(nmonth) + 1 ) THEN      ! NEW month 
     225         IF( nday == nmonth_len(nmonth) + 1 ) THEN      ! New month 
    215226            nday   = 1 
    216227            nmonth = nmonth + 1 
    217228            nsec_month = ndt05 
    218             IF( nmonth == 13 ) THEN                     ! NEW year 
     229            IF( nmonth == 13 ) THEN                     ! New year 
    219230               nyear     = nyear + 1 
    220231               nmonth    = 1 
    … …  
    226237         ENDIF 
    227238         ! 
    228          ndastp = nyear * 10000 + nmonth * 100 + nday   ! NEW date 
     239         ndastp = nyear * 10000 + nmonth * 100 + nday   ! New date 
     240         ! 
     241         !compute first day of the year in julian days 
     242         CALL ymds2ju( nyear, 01, 01, 0.0, fjulstartyear ) 
    229243         ! 
    230244         IF(lwp) WRITE(numout,'(a,i8,a,i4.4,a,i2.2,a,i2.2,a,i3.3)') '======>> time-step =', kt,   & 
    231245              &   '      New day, DATE Y/M/D = ', nyear, '/', nmonth, '/', nday, '      nday_year = ', nday_year 
    232246         IF(lwp) WRITE(numout,'(a,i8,a,i7,a,i5)') '         nsec_year = ', nsec_year,   & 
    233               &   '   nsec_month = ', nsec_month, '   nsec_day = ', nsec_day 
    234       ENDIF 
     247              &   '   nsec_month = ', nsec_month, '   nsec_day = ', nsec_day, '   nsec_week = ', nsec_week 
     248      ENDIF 
     249 
     250      IF( nsec_week > 7*nsecd )   nsec_week = ndt05     ! New week 
    235251       
    236252      IF(ln_ctl) THEN 
    … …  
    239255      ENDIF 
    240256 
    241       CALL rst_opn( kt )                                ! Open the restart file if needed and control lrst_oce 
    242       IF( lrst_oce )  CALL day_rst( kt, 'WRITE' )      ! write day restart information 
     257      IF( .NOT. lk_offline ) CALL rst_opn( kt )               ! Open the restart file if needed and control lrst_oce 
     258      IF( lrst_oce         ) CALL day_rst( kt, 'WRITE' )      ! write day restart information 
    243259      ! 
    244260   END SUBROUTINE day 

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

    • Property svn:eol-style deleted
    r1886 r2528  
    66   !!====================================================================== 
    77   !! History :  1.0  ! 2005-10  (A. Beckmann, G. Madec)  reactivate s-coordinate  
    8    !!---------------------------------------------------------------------- 
    9    !! NEMO/OPA 3.2 , LOCEAN-IPSL (2009)  
    10    !! $Id$  
    11    !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 
     8   !!            3.3  ! 2010-11  (G. Madec) add mbk. arrays associated to the deepest ocean level 
    129   !!---------------------------------------------------------------------- 
    1310   USE par_oce      ! ocean parameters 
    … …  
    2118   !! ---------------------------- 
    2219   !                                              !!* Namelist namdom : time & space domain * 
    23    INTEGER , PUBLIC ::   nn_bathy     =  0         !: = 0/1 ,compute/read the bathymetry file 
    24    REAL(wp), PUBLIC ::   rn_e3zps_min = 5.0_wp     !: miminum thickness for partial steps (meters) 
    25    REAL(wp), PUBLIC ::   rn_e3zps_rat = 0.1_wp     !: minimum thickness ration for partial steps 
    26    INTEGER , PUBLIC ::   nn_msh       = 0          !: = 1 create a mesh-mask file 
    27    INTEGER , PUBLIC ::   nn_acc       = 0          !: = 0/1 use of the acceleration of convergence technique 
    28    REAL(wp), PUBLIC ::   rn_atfp      = 0.1_wp     !: asselin time filter parameter 
     20   INTEGER , PUBLIC ::   nn_bathy     =    0       !: = 0/1 ,compute/read the bathymetry file 
     21   REAL(wp), PUBLIC ::   rn_hmin      =   -3.0_wp  !: minimum ocean depth (>0) or minimum number of ocean levels (<0) 
     22   REAL(wp), PUBLIC ::   rn_e3zps_min =    5.0_wp  !: miminum thickness for partial steps (meters) 
     23   REAL(wp), PUBLIC ::   rn_e3zps_rat =    0.1_wp  !: minimum thickness ration for partial steps 
     24   INTEGER , PUBLIC ::   nn_msh       =    0       !: = 1 create a mesh-mask file 
     25   INTEGER , PUBLIC ::   nn_acc       =    0       !: = 0/1 use of the acceleration of convergence technique 
     26   REAL(wp), PUBLIC ::   rn_atfp      =    0.1_wp  !: asselin time filter parameter 
    2927   REAL(wp), PUBLIC ::   rn_rdt       = 3600._wp   !: time step for the dynamics (and tracer if nacc=0) 
    3028   REAL(wp), PUBLIC ::   rn_rdtmin    = 3600._wp   !: minimum time step on tracers 
    3129   REAL(wp), PUBLIC ::   rn_rdtmax    = 3600._wp   !: maximum time step on tracers 
    3230   REAL(wp), PUBLIC ::   rn_rdth      =  800._wp   !: depth variation of tracer step 
    33    INTEGER , PUBLIC ::   nn_baro      = 64         !: number of barotropic time steps (key_dynspg_ts) 
    34    INTEGER , PUBLIC ::   nn_closea    =  0         !: =0 suppress closed sea/lake from the ORCA domain or not (=1) 
    35  
    36    !                                          ! old non-DOCTOR names still used in the model 
    37    INTEGER , PUBLIC ::   ntopo                !: = 0/1 ,compute/read the bathymetry file 
    38    REAL(wp), PUBLIC ::   e3zps_min            !: miminum thickness for partial steps (meters) 
    39    REAL(wp), PUBLIC ::   e3zps_rat            !: minimum thickness ration for partial steps 
    40    INTEGER , PUBLIC ::   nmsh                 !: = 1 create a mesh-mask file 
    41    INTEGER , PUBLIC ::   nacc                 !: = 0/1 use of the acceleration of convergence technique 
    42    REAL(wp), PUBLIC ::   atfp                 !: asselin time filter parameter 
    43    REAL(wp), PUBLIC ::   rdt                  !: time step for the dynamics (and tracer if nacc=0) 
    44    REAL(wp), PUBLIC ::   rdtmin               !: minimum time step on tracers 
    45    REAL(wp), PUBLIC ::   rdtmax               !: maximum time step on tracers 
    46    REAL(wp), PUBLIC ::   rdth                 !: depth variation of tracer step 
    47    INTEGER , PUBLIC ::   nclosea              !: =0 suppress closed sea/lake from the ORCA domain or not (=1) 
    48  
    49  
    50    !                                         !!! associated variables 
     31   INTEGER , PUBLIC ::   nn_baro      =   64       !: number of barotropic time steps (key_dynspg_ts) 
     32   INTEGER , PUBLIC ::   nn_closea    =    0       !: =0 suppress closed sea/lake from the ORCA domain or not (=1) 
     33 
     34   !                                    !! old non-DOCTOR names still used in the model 
     35   INTEGER , PUBLIC ::   ntopo           !: = 0/1 ,compute/read the bathymetry file 
     36   REAL(wp), PUBLIC ::   e3zps_min       !: miminum thickness for partial steps (meters) 
     37   REAL(wp), PUBLIC ::   e3zps_rat       !: minimum thickness ration for partial steps 
     38   INTEGER , PUBLIC ::   nmsh            !: = 1 create a mesh-mask file 
     39   INTEGER , PUBLIC ::   nacc            !: = 0/1 use of the acceleration of convergence technique 
     40   REAL(wp), PUBLIC ::   atfp            !: asselin time filter parameter 
     41   REAL(wp), PUBLIC ::   rdt             !: time step for the dynamics (and tracer if nacc=0) 
     42   REAL(wp), PUBLIC ::   rdtmin          !: minimum time step on tracers 
     43   REAL(wp), PUBLIC ::   rdtmax          !: maximum time step on tracers 
     44   REAL(wp), PUBLIC ::   rdth            !: depth variation of tracer step 
     45   INTEGER , PUBLIC ::   nclosea         !: =0 suppress closed sea/lake from the ORCA domain or not (=1) 
     46 
     47 
     48   !                                                  !!! associated variables 
    5149   INTEGER , PUBLIC                 ::   neuler  = 0   !: restart euler forward option (0=Euler) 
    5250   REAL(wp), PUBLIC                 ::   atfp1         !: asselin time filter coeff. (atfp1= 1-2*atfp) 
    … …  
    5553   !                                         !!* Namelist namcla : cross land advection 
    5654   INTEGER, PUBLIC ::   nn_cla = 0            !: =1 cross land advection for exchanges through some straits (ORCA2) 
    57  
    58    !                                          ! old non-DOCTOR names still used in the model 
    59    INTEGER, PUBLIC ::   n_cla = 0             !: =1 cross land advection for exchanges through some straits (ORCA2) 
    6055 
    6156   !!---------------------------------------------------------------------- 
    … …  
    120115   LOGICAL, PUBLIC ::   ln_zps     =  .FALSE.   !: z-coordinate - partial step 
    121116   LOGICAL, PUBLIC ::   ln_sco     =  .FALSE.   !: s-coordinate or hybrid z-s coordinate 
    122 #if defined key_zco 
    123    LOGICAL, PUBLIC, PARAMETER ::   lk_zco = .TRUE.    !: z-coordinate flag (1D arrays) 
    124 #else 
    125    LOGICAL, PUBLIC, PARAMETER ::   lk_zco = .FALSE.   !: z-coordinate flag (3D arrays) 
    126117 
    127118   !! All coordinates 
    … …  
    133124   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   e3vw            !: analytical vertical scale factors at  VW-- 
    134125   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   e3w   , e3uw    !:                                        W--UW  points (m) 
    135 #endif 
    136126#if defined key_vvl 
    137127   LOGICAL, PUBLIC, PARAMETER ::   lk_vvl = .TRUE.    !: variable grid flag 
    … …  
    145135   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   e3vw_1             !: analytical vertical scale factors at  VW-- 
    146136   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   e3w_1  , e3uw_1    !:                                       W--UW  points (m) 
     137   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   e3t_b              !: before         -      -      -    -   T      points (m) 
     138   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   e3u_b  , e3v_b     !:   -            -      -      -    -   U--V   points (m) 
    147139#else 
    148140   LOGICAL, PUBLIC, PARAMETER ::   lk_vvl = .FALSE.   !: fixed grid flag 
    … …  
    159151   REAL(wp), PUBLIC, DIMENSION(jpk)     ::   gdept_0, gdepw_0   !: reference depth of t- and w-points (m) 
    160152   REAL(wp), PUBLIC, DIMENSION(jpk)     ::   e3t_0  , e3w_0     !: reference vertical scale factors at T- and W-pts (m) 
    161    REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   hdept  , hdepw     !: ocean bottom depth at T and W points 
    162153   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   e3tp   , e3wp      !: ocean bottom level thickness at T and W points 
    163154 
    … …  
    172163   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   scosrf, scobot   !: ocean surface and bottom topographies  
    173164   !                                                          !  (if deviating from coordinate surfaces in HYBRID) 
    174    REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   hifv  , hiff     !: interface depth between stretching    at  V--F 
    175    REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   hift  , hifu     !: and quasi-uniform spacing                 T--U  points (m) 
     165   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   hifv  , hiff     !: interface depth between stretching at  V--F 
     166   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   hift  , hifu     !: and quasi-uniform spacing              T--U  points (m) 
    176167 
    177168   !!---------------------------------------------------------------------- 
    178169   !! masks, bathymetry 
    179170   !! --------------------------------------------------------------------- 
    180    INTEGER , PUBLIC, DIMENSION(jpi,jpj) ::   mbathy    !: number of ocean level (=0, 1, ... , jpk-1) 
    181    REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   bathy     !: ocean depth (meters) 
    182    REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   tmask_i   !: interior domain T-point mask 
    183    REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   bmask     !: land/ocean mask of barotropic stream function 
    184  
    185    REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   tmask, umask, vmask, fmask   !: land/ocean mask at T-, U-, V- and F-points 
    186  
    187    REAL(wp), PUBLIC, DIMENSION(jpiglo) ::   tpol, fpol          !: north fold mask (nperio= 3 or 4) 
     171   INTEGER , PUBLIC, DIMENSION(jpi,jpj) ::   mbathy       !: number of ocean level (=0, 1, ... , jpk-1) 
     172   INTEGER , PUBLIC, DIMENSION(jpi,jpj) ::   mbkt         !: vertical index of the bottom last T- ocean level 
     173   INTEGER , PUBLIC, DIMENSION(jpi,jpj) ::   mbku, mbkv   !: vertical index of the bottom last U- and W- ocean level 
     174   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   bathy        !: ocean depth (meters) 
     175   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   tmask_i      !: interior domain T-point mask 
     176   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   bmask        !: land/ocean mask of barotropic stream function 
     177 
     178   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   tmask, umask, vmask, fmask   !: land/ocean mask at T-, U-, V- and F-pts 
     179 
     180   REAL(wp), PUBLIC, DIMENSION(jpiglo) ::   tpol, fpol          !: north fold mask (jperio= 3 or 4) 
    188181 
    189182#if defined key_noslip_accurate 
    … …  
    195188   !! calendar variables 
    196189   !! --------------------------------------------------------------------- 
    197    INTEGER , PUBLIC ::   nyear       !: current year 
    198    INTEGER , PUBLIC ::   nmonth      !: current month 
    199    INTEGER , PUBLIC ::   nday        !: current day of the month 
    200    INTEGER , PUBLIC ::   ndastp      !: time step date in yyyymmdd format 
    201    INTEGER , PUBLIC ::   nday_year   !: current day counted from jan 1st of the current year 
    202    INTEGER , PUBLIC ::   nsec_year   !: current time step counted in second since 00h jan 1st of the current year 
    203    INTEGER , PUBLIC ::   nsec_month  !: current time step counted in second since 00h 1st day of the current month 
    204    INTEGER , PUBLIC ::   nsec_day    !: current time step counted in second since 00h of the current day 
    205    REAL(wp), PUBLIC ::   fjulday     !: julian day  
    206    REAL(wp), PUBLIC ::   adatrj      !: number of elapsed days since the begining of the whole simulation 
    207    !                                 !: (cumulative duration of previous runs that may have used different time-step size) 
     190   INTEGER , PUBLIC ::   nyear         !: current year 
     191   INTEGER , PUBLIC ::   nmonth        !: current month 
     192   INTEGER , PUBLIC ::   nday          !: current day of the month 
     193   INTEGER , PUBLIC ::   ndastp        !: time step date in yyyymmdd format 
     194   INTEGER , PUBLIC ::   nday_year     !: current day counted from jan 1st of the current year 
     195   INTEGER , PUBLIC ::   nsec_year     !: current time step counted in second since 00h jan 1st of the current year 
     196   INTEGER , PUBLIC ::   nsec_month    !: current time step counted in second since 00h 1st day of the current month 
     197   INTEGER , PUBLIC ::   nsec_week     !: current time step counted in second since 00h of last monday 
     198   INTEGER , PUBLIC ::   nsec_day      !: current time step counted in second since 00h of the current day 
     199   REAL(wp), PUBLIC ::   fjulday       !: current julian day  
     200   REAL(wp), PUBLIC ::   fjulstartyear !: first day of the current year in julian days 
     201   REAL(wp), PUBLIC ::   adatrj        !: number of elapsed days since the begining of the whole simulation 
     202   !                                   !: (cumulative duration of previous runs that may have used different time-step size) 
    208203   INTEGER , PUBLIC, DIMENSION(0: 1) ::   nyear_len     !: length in days of the previous/current year 
    209204   INTEGER , PUBLIC, DIMENSION(0:13) ::   nmonth_len    !: length in days of the months of the current year 
    … …  
    213208 
    214209   !!---------------------------------------------------------------------- 
     210   !! mpp reproducibility 
     211   !!---------------------------------------------------------------------- 
     212#if defined key_mpp_rep 
     213   LOGICAL, PUBLIC, PARAMETER ::   lk_mpp_rep = .TRUE.    !: agrif flag 
     214#else 
     215   LOGICAL, PUBLIC, PARAMETER ::   lk_mpp_rep = .FALSE.   !: agrif flag 
     216#endif 
     217   !!---------------------------------------------------------------------- 
    215218   !! agrif domain 
    216219   !!---------------------------------------------------------------------- 
    … …  
    229232   END FUNCTION Agrif_CFixed 
    230233#endif 
    231  
     234   !!---------------------------------------------------------------------- 
     235   !! NEMO/OPA 3.3 , NEMO Consortium (2010) 
     236   !! $Id$  
     237   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
    232238   !!====================================================================== 
    233239END MODULE dom_oce 

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

    • Property svn:eol-style deleted
    r1792 r2528  
    1111   !!   NEMO     1.0  !  2002-08  (G. Madec)  F90: Free form and module 
    1212   !!            2.0  !  2005-11  (V. Garnier) Surface pressure gradient organization 
     13   !!            3.3  !  2010-11  (G. Madec)  initialisation in C1D configuration 
    1314   !!---------------------------------------------------------------------- 
    1415    
    … …  
    1819   !!   dom_ctl        : control print for the ocean domain 
    1920   !!---------------------------------------------------------------------- 
    20    USE oce             !  
    21    USE dom_oce         ! ocean space and time domain 
     21   USE oce             ! ocean variables 
     22   USE dom_oce         ! domain: ocean 
    2223   USE sbc_oce         ! surface boundary condition: ocean 
    2324   USE phycst          ! physical constants 
    … …  
    3233   USE domwri          ! domain: write the meshmask file 
    3334   USE domvvl          ! variable volume 
     35   USE c1d             ! 1D vertical configuration 
     36   USE dyncor_c1d      ! Coriolis term (c1d case)         (cor_c1d routine) 
    3437 
    3538   IMPLICIT NONE 
    … …  
    4144#  include "domzgr_substitute.h90" 
    4245   !!------------------------------------------------------------------------- 
    43    !! NEMO/OPA 3.2 , LOCEAN-IPSL (2009)  
     46   !! NEMO/OPA 3.3 , NEMO Consortium (2010) 
    4447   !! $Id$ 
    45    !! Software is governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 
     48   !! Software governed by the CeCILL licence        (NEMOGCM/NEMO_CeCILL.txt) 
    4649   !!------------------------------------------------------------------------- 
    47  
    4850CONTAINS 
    4951 
    … …  
    6264      !!              - dom_stp: defined the model time step 
    6365      !!              - dom_wri: create the meshmask file if nmsh=1 
     66      !!              - 1D configuration, move Coriolis, u and v at T-point 
    6467      !!---------------------------------------------------------------------- 
    6568      INTEGER ::   jk                ! dummy loop argument 
    … …  
    7982                             CALL dom_msk      ! Masks 
    8083      IF( lk_vvl         )   CALL dom_vvl      ! Vertical variable mesh 
     84      ! 
     85      IF( lk_c1d ) THEN                        ! 1D configuration  
     86         CALL cor_c1d                          ! Coriolis set at T-point 
     87         umask(:,:,:) = tmask(:,:,:)           ! U, V moved at T-point 
     88         vmask(:,:,:) = tmask(:,:,:) 
     89      END IF 
    8190      ! 
    8291      hu(:,:) = 0.e0                           ! Ocean depth at U- and V-points 
    … …  
    106115      !!              - namdom namelist 
    107116      !!              - namcla namelist 
     117      !!              - namnc4 namelist   ! "key_netcdf4" only 
    108118      !!---------------------------------------------------------------------- 
    109119      USE ioipsl 
    … …  
    111121         &             nn_it000, nn_itend  , nn_date0    , nn_leapy     , nn_istate , nn_stock ,   & 
    112122         &             nn_write, ln_dimgnnn, ln_mskland  , ln_clobber   , nn_chunksz 
    113       NAMELIST/namdom/ nn_bathy , rn_e3zps_min, rn_e3zps_rat, nn_msh   ,   & 
    114          &             nn_acc   , rn_atfp     , rn_rdt      , rn_rdtmin,   & 
     123      NAMELIST/namdom/ nn_bathy , rn_e3zps_min, rn_e3zps_rat, nn_msh    , rn_hmin,   & 
     124         &             nn_acc   , rn_atfp     , rn_rdt      , rn_rdtmin ,            & 
    115125         &             rn_rdtmax, rn_rdth     , nn_baro     , nn_closea 
    116126      NAMELIST/namcla/ nn_cla 
     127#if defined key_netcdf4 
     128      NAMELIST/namnc4/ nn_nchunks_i, nn_nchunks_j, nn_nchunks_k, ln_nc4zip 
     129#endif 
    117130      !!---------------------------------------------------------------------- 
    118131 
    … …  
    166179      ENDIF 
    167180 
     181#if defined key_agrif 
    168182      IF( Agrif_Root() ) THEN 
    169          SELECT CASE ( nleapy )        ! Choose calendar for IOIPSL 
    170          CASE (  1 )  
    171             CALL ioconf_calendar('gregorian') 
    172             IF(lwp) WRITE(numout,*) '   The IOIPSL calendar is "gregorian", i.e. leap year' 
    173          CASE (  0 ) 
    174             CALL ioconf_calendar('noleap') 
    175             IF(lwp) WRITE(numout,*) '   The IOIPSL calendar is "noleap", i.e. no leap year' 
    176          CASE ( 30 ) 
    177             CALL ioconf_calendar('360d') 
    178             IF(lwp) WRITE(numout,*) '   The IOIPSL calendar is "360d", i.e. 360 days in a year' 
    179          END SELECT 
    180       ENDIF 
    181  
    182       REWIND( numnam )             ! Namelist namdom : space & time domain (bathymetry, mesh, timestep) 
     183#endif 
     184      SELECT CASE ( nleapy )        ! Choose calendar for IOIPSL 
     185      CASE (  1 )  
     186         CALL ioconf_calendar('gregorian') 
     187         IF(lwp) WRITE(numout,*) '   The IOIPSL calendar is "gregorian", i.e. leap year' 
     188      CASE (  0 ) 
     189         CALL ioconf_calendar('noleap') 
     190         IF(lwp) WRITE(numout,*) '   The IOIPSL calendar is "noleap", i.e. no leap year' 
     191      CASE ( 30 ) 
     192         CALL ioconf_calendar('360d') 
     193         IF(lwp) WRITE(numout,*) '   The IOIPSL calendar is "360d", i.e. 360 days in a year' 
     194      END SELECT 
     195#if defined key_agrif 
     196      ENDIF 
     197#endif 
     198 
     199      REWIND( numnam )              ! Namelist namdom : space & time domain (bathymetry, mesh, timestep) 
    183200      READ  ( numnam, namdom ) 
    184201 
    … …  
    187204         WRITE(numout,*) '   Namelist namdom : space & time domain' 
    188205         WRITE(numout,*) '      flag read/compute bathymetry      nn_bathy     = ', nn_bathy 
     206         WRITE(numout,*) '      min depth of the ocean    (>0) or    rn_hmin   = ', rn_hmin 
     207         WRITE(numout,*) '      min number of ocean level (<0)       ' 
    189208         WRITE(numout,*) '      minimum thickness of partial      rn_e3zps_min = ', rn_e3zps_min, ' (m)' 
    190209         WRITE(numout,*) '         step level                     rn_e3zps_rat = ', rn_e3zps_rat 
    191210         WRITE(numout,*) '      create mesh/mask file(s)          nn_msh       = ', nn_msh 
    192          WRITE(numout,*) '           = 0   no file created                 ' 
    193          WRITE(numout,*) '           = 1   mesh_mask                       ' 
    194          WRITE(numout,*) '           = 2   mesh and mask                   ' 
    195          WRITE(numout,*) '           = 3   mesh_hgr, msh_zgr and mask      ' 
     211         WRITE(numout,*) '           = 0   no file created           ' 
     212         WRITE(numout,*) '           = 1   mesh_mask                 ' 
     213         WRITE(numout,*) '           = 2   mesh and mask             ' 
     214         WRITE(numout,*) '           = 3   mesh_hgr, msh_zgr and mask' 
    196215         WRITE(numout,*) '      ocean time step                      rn_rdt    = ', rn_rdt 
    197216         WRITE(numout,*) '      asselin time filter parameter        rn_atfp   = ', rn_atfp 
    … …  
    216235      nclosea   = nn_closea 
    217236 
    218       REWIND( numnam )             ! Namelist cross land advection 
     237      REWIND( numnam )              ! Namelist cross land advection 
    219238      READ  ( numnam, namcla ) 
    220239      IF(lwp) THEN 
    … …  
    224243      ENDIF 
    225244 
    226       n_cla = nn_cla                ! conversion DOCTOR names into model names (this should disappear soon) 
    227  
    228       IF( nbit_cmp == 1 .AND. n_cla /= 0 )   CALL ctl_stop( ' Reproductibility tests (nbit_cmp=1) require n_cla = 0' ) 
     245#if defined key_netcdf4 
     246      !                             ! NetCDF 4 case   ("key_netcdf4" defined) 
     247      REWIND( numnam )                    ! Namelist namnc4 : netcdf4 chunking parameters 
     248      READ  ( numnam, namnc4 ) 
     249      IF(lwp) THEN                        ! control print 
     250         WRITE(numout,*) 
     251         WRITE(numout,*) '   Namelist namnc4 - Netcdf4 chunking parameters' 
     252         WRITE(numout,*) '      number of chunks in i-dimension      nn_nchunks_i   = ', nn_nchunks_i 
     253         WRITE(numout,*) '      number of chunks in j-dimension      nn_nchunks_j   = ', nn_nchunks_j 
     254         WRITE(numout,*) '      number of chunks in k-dimension      nn_nchunks_k   = ', nn_nchunks_k 
     255         WRITE(numout,*) '      apply netcdf4/hdf5 chunking & compression ln_nc4zip = ', ln_nc4zip 
     256      ENDIF 
     257 
     258      ! Put the netcdf4 settings into a simple structure (snc4set, defined in in_out_manager module) 
     259      ! Note the chunk size in the unlimited (time) dimension will be fixed at 1 
     260      snc4set%ni   = nn_nchunks_i 
     261      snc4set%nj   = nn_nchunks_j 
     262      snc4set%nk   = nn_nchunks_k 
     263      snc4set%luse = ln_nc4zip 
     264#else 
     265      snc4set%luse = .FALSE.        ! No NetCDF 4 case 
     266#endif 
    229267      ! 
    230268   END SUBROUTINE dom_nam 

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

    • Property svn:eol-style deleted
    r1566 r2528  
    2424   !! NEMO/OPA 3.2 , LODYC-IPSL  (2009) 
    2525   !! $Id$  
    26    !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)  
     26   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 
    2727   !!---------------------------------------------------------------------- 
    2828 

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

    • Property svn:eol-style deleted
    r1792 r2528  
    44   !! Ocean initialization : domain initialization 
    55   !!============================================================================== 
    6    !! History :       !  88-03  (G. Madec) 
    7    !!                 !  91-11  (G. Madec) 
    8    !!                 !  92-06  (M. Imbard) 
    9    !!                 !  96-01  (G. Madec)  terrain following coordinates 
    10    !!                 !  97-02  (G. Madec)  print mesh informations 
    11    !!                 !  99-11  (M. Imbard) NetCDF format with IO-IPSL 
    12    !!                 !  00-08  (D. Ludicone) Reduced section at Bab el Mandeb 
    13    !!                 !  01-09  (M. Levy)  eel config: grid in km, beta-plane 
    14    !!            8.5  !  02-08  (G. Madec)  F90: Free form and module, namelist 
    15    !!            9.0  !  04-01  (A.M. Treguier, J.M. Molines) Case 4 (Mercator mesh) 
    16    !!                           use of parameters in par_CONFIG-Rxx.h90, not in namelist 
    17    !!                 !  04-05  (A. Koch-Larrouy) Add Gyre configuration  
     6   !! History :  OPA  ! 1988-03  (G. Madec) Original code 
     7   !!            7.0  ! 1996-01  (G. Madec)  terrain following coordinates 
     8   !!            8.0  ! 1997-02  (G. Madec)  print mesh informations 
     9   !!            8.1  ! 1999-11  (M. Imbard) NetCDF format with IO-IPSL 
     10   !!            8.2  ! 2000-08  (D. Ludicone) Reduced section at Bab el Mandeb 
     11   !!             -   ! 2001-09  (M. Levy)  eel config: grid in km, beta-plane 
     12   !!  NEMO      1.0  ! 2002-08  (G. Madec)  F90: Free form and module, namelist 
     13   !!             -   ! 2004-01  (A.M. Treguier, J.M. Molines) Case 4 (Mercator mesh) 
     14   !!                            use of parameters in par_CONFIG-Rxx.h90, not in namelist 
     15   !!             -   ! 2004-05  (A. Koch-Larrouy) Add Gyre configuration  
    1816   !!---------------------------------------------------------------------- 
    1917 
    2018   !!---------------------------------------------------------------------- 
    21    !!   dom_hgr        : initialize the horizontal mesh  
    22    !!   hgr_read       : read "coordinate" NetCDF file  
     19   !!   dom_hgr       : initialize the horizontal mesh  
     20   !!   hgr_read      : read "coordinate" NetCDF file  
    2321   !!---------------------------------------------------------------------- 
    24    !! * Modules used 
    25    USE dom_oce         ! ocean space and time domain 
    26    USE phycst          ! physical constants 
    27    USE in_out_manager  ! I/O manager 
    28    USE lib_mpp 
     22   USE dom_oce        ! ocean space and time domain 
     23   USE phycst         ! physical constants 
     24   USE in_out_manager ! I/O manager 
     25   USE lib_mpp        ! MPP library 
    2926 
    3027   IMPLICIT NONE 
    3128   PRIVATE 
    3229 
    33    !! * Module variables 
    34    REAL(wp) ::   glam0, gphi0           ! variables corresponding to parameters 
    35       !                                 ! ppglam0 ppgphi0 set in par_oce 
    36  
    37    !! * Routine accessibility 
    38    PUBLIC dom_hgr        ! called by domain.F90 
     30   REAL(wp) ::   glam0, gphi0   ! variables corresponding to parameters ppglam0 ppgphi0 set in par_oce 
     31 
     32   PUBLIC   dom_hgr   ! called by domain.F90 
     33 
    3934   !!---------------------------------------------------------------------- 
    40    !!   OPA 9.0 , LOCEAN-IPSL (2005)  
     35   !! NEMO/OPA 3.3 , NEMO Consortium (2010) 
    4136   !! $Id$  
    42    !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 
     37   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
    4338   !!---------------------------------------------------------------------- 
    44  
    4539CONTAINS 
    4640 
    … …  
    10094      !!                Madec, Imbard, 1996, Clim. Dyn. 
    10195      !!---------------------------------------------------------------------- 
    102       INTEGER  ::   ji, jj              ! dummy loop indices 
    103       INTEGER  ::   ii0, ii1, ij0, ij1  ! temporary integers 
    104       INTEGER  ::   ijeq                ! index of equator T point (used in case 4) 
    105       REAL(wp) ::   & 
    106          zti, zui, zvi, zfi,         &  ! temporary scalars 
    107          ztj, zuj, zvj, zfj,         &  ! 
    108          zphi0, zbeta, znorme,       &  ! 
    109          zarg, zf0, zminff, zmaxff 
    110       REAL(wp) ::   & 
    111          zlam1, zcos_alpha, zim1 , zjm1 , ze1, ze1deg,   & 
    112          zphi1, zsin_alpha, zim05, zjm05 
     96      INTEGER  ::   ji, jj               ! dummy loop indices 
     97      INTEGER  ::   ii0, ii1, ij0, ij1   ! temporary integers 
     98      INTEGER  ::   ijeq                 ! index of equator T point (used in case 4) 
     99      REAL(wp) ::   zti, zui, zvi, zfi   ! local scalars 
     100      REAL(wp) ::   ztj, zuj, zvj, zfj   !   -      - 
     101      REAL(wp) ::   zphi0, zbeta, znorme ! 
     102      REAL(wp) ::   zarg, zf0, zminff, zmaxff 
     103      REAL(wp) ::   zlam1, zcos_alpha, zim1 , zjm1 , ze1, ze1deg 
     104      REAL(wp) ::   zphi1, zsin_alpha, zim05, zjm05 
    113105      !!---------------------------------------------------------------------- 
    114106 
    … …  
    138130         IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN    ! ORCA R2 configuration 
    139131            !                                             ! ===================== 
    140             IF( n_cla == 0 ) THEN 
     132            IF( nn_cla == 0 ) THEN 
    141133               ! 
    142134               ii0 = 139   ;   ii1 = 140        ! Gibraltar Strait (e2u = 20 km) 
    … …  
    157149            IF(lwp) WRITE(numout,*) 
    158150            IF(lwp) WRITE(numout,*) '             orca_r2: Danish Straits : e2u reduced to 10 km' 
     151            ! 
     152         ENDIF 
     153 
     154            !                                             ! ===================== 
     155         IF( cp_cfg == "orca" .AND. jp_cfg == 1 ) THEN    ! ORCA R1 configuration 
     156            !                                             ! ===================== 
     157 
     158            ii0 = 281   ;   ii1 = 282        ! Gibraltar Strait (e2u = 20 km) 
     159            ij0 = 200   ;   ij1 = 200   ;   e2u( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) =  20.e3 
     160            IF(lwp) WRITE(numout,*) 
     161            IF(lwp) WRITE(numout,*) '             orca_r1: Gibraltar : e2u reduced to 20 km' 
     162 
     163            ii0 = 314   ;   ii1 = 315        ! Bhosporus Strait (e2u = 10 km) 
     164            ij0 = 208   ;   ij1 = 208   ;   e2u( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) =  10.e3 
     165            IF(lwp) WRITE(numout,*) 
     166            IF(lwp) WRITE(numout,*) '             orca_r1: Bhosporus : e2u reduced to 10 km' 
     167 
     168            ii0 =  44   ;   ii1 =  44        ! Lombok Strait (e1v = 13 km) 
     169            ij0 = 124   ;   ij1 = 125   ;   e1v( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) =  13.e3 
     170            IF(lwp) WRITE(numout,*) 
     171            IF(lwp) WRITE(numout,*) '             orca_r1: Lombok : e1v reduced to 10 km' 
     172 
     173            ii0 =  48   ;   ii1 =  48        ! Sumba Strait (e1v = 8 km) [closed from bathy_11 on] 
     174            ij0 = 124   ;   ij1 = 125   ;   e1v( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) =  8.e3 
     175            IF(lwp) WRITE(numout,*) 
     176            IF(lwp) WRITE(numout,*) '             orca_r1: Sumba : e1v reduced to 8 km' 
     177 
     178            ii0 =  53   ;   ii1 =  53        ! Ombai Strait (e1v = 13 km) 
     179            ij0 = 124   ;   ij1 = 125   ;   e1v( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 13.e3 
     180            IF(lwp) WRITE(numout,*) 
     181            IF(lwp) WRITE(numout,*) '             orca_r1: Ombai : e1v reduced to 13 km' 
     182 
     183            ii0 =  56   ;   ii1 =  56        ! Timor Passage (e1v = 20 km) 
     184            ij0 = 124   ;   ij1 = 125   ;   e1v( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 20.e3 
     185            IF(lwp) WRITE(numout,*) 
     186            IF(lwp) WRITE(numout,*) '             orca_r1: Timor Passage : e1v reduced to 20 km' 
     187 
     188            ii0 =  55   ;   ii1 =  55        ! West Halmahera Strait (e1v = 30 km) 
     189            ij0 = 141   ;   ij1 = 142   ;   e1v( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 30.e3 
     190            IF(lwp) WRITE(numout,*) 
     191            IF(lwp) WRITE(numout,*) '             orca_r1: W Halmahera : e1v reduced to 30 km' 
     192 
     193            ii0 =  58   ;   ii1 =  58        ! East Halmahera Strait (e1v = 50 km) 
     194            ij0 = 141   ;   ij1 = 142   ;   e1v( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 50.e3 
     195            IF(lwp) WRITE(numout,*) 
     196            IF(lwp) WRITE(numout,*) '             orca_r1: E Halmahera : e1v reduced to 50 km' 
     197 
     198            ! 
     199 
     200            ! 
     201            ! 
     202            ! 
    159203            ! 
    160204         ENDIF 

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

    • Property svn:eol-style deleted
    r1707 r2528  
    55   !!====================================================================== 
    66   !! History :  OPA  ! 1987-07  (G. Madec)  Original code 
    7    !!             -   ! 1993-03  (M. Guyon)  symetrical conditions (M. Guyon) 
    8    !!             -   ! 1996-01  (G. Madec)  suppression of common work arrays 
     7   !!            6.0  ! 1993-03  (M. Guyon)  symetrical conditions (M. Guyon) 
     8   !!            7.0  ! 1996-01  (G. Madec)  suppression of common work arrays 
    99   !!             -   ! 1996-05  (G. Madec)  mask computed from tmask and sup- 
    1010   !!                 !                      pression of the double computation of bmask 
    11    !!             -   ! 1997-02  (G. Madec)  mesh information put in domhgr.F 
    12    !!             -   ! 1997-07  (G. Madec)  modification of mbathy and fmask 
     11   !!            8.0  ! 1997-02  (G. Madec)  mesh information put in domhgr.F 
     12   !!            8.1  ! 1997-07  (G. Madec)  modification of mbathy and fmask 
    1313   !!             -   ! 1998-05  (G. Roullet)  free surface 
    14    !!             -   ! 2000-03  (G. Madec)  no slip accurate 
     14   !!            8.2  ! 2000-03  (G. Madec)  no slip accurate 
    1515   !!             -   ! 2001-09  (J.-M. Molines)  Open boundaries 
    1616   !!   NEMO     1.0  ! 2002-08  (G. Madec)  F90: Free form and module 
    … …  
    4444   !! NEMO/OPA 3.2 , LODYC-IPSL  (2009) 
    4545   !! $Id$  
    46    !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)  
    47    !!---------------------------------------------------------------------- 
    48  
     46   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
     47   !!---------------------------------------------------------------------- 
    4948CONTAINS 
    5049    
    … …  
    5655      !!      zontal velocity points (u & v), vorticity points (f) and baro- 
    5756      !!      tropic stream function  points (b). 
    58       !!        Set mbathy to the number of non-zero w-levels of a water column 
    5957      !! 
    6058      !! ** Method  :   The ocean/land mask is computed from the basin bathy- 
    … …  
    7371      !!                   or mbathy(ji+1,jj)  or mbathy(ji+1,jj+1) =< 0 
    7472      !!                1. IF mbathy( ji ,jj) and mbathy( ji ,jj+1) 
    75       !!                and mbathy(ji+1,jj) and mbathy(ji+1,jj+1) >= jk. 
     73      !!                  and mbathy(ji+1,jj) and mbathy(ji+1,jj+1) >= jk. 
    7674      !!      b-point : the same definition as for f-point of the first ocean 
    7775      !!                level (surface level) but with 0 along coastlines. 
     76      !!      tmask_i : interior ocean mask at t-point, i.e. excluding duplicated 
     77      !!                rows/lines due to cyclic or North Fold boundaries as well 
     78      !!                as MPP halos. 
    7879      !! 
    7980      !!        The lateral friction is set through the value of fmask along 
    … …  
    99100      !!        - bmask is  set to 0 on the open boundaries. 
    100101      !! 
    101       !!      Set mbathy to the number of non-zero w-levels of a water column 
    102       !!                  mbathy = min( mbathy, 1 ) + 1 
    103       !!      (note that the minimum value of mbathy is 2). 
    104       !! 
    105102      !! ** Action :   tmask    : land/ocean mask at t-point (=0. or 1.) 
    106103      !!               umask    : land/ocean mask at u-point (=0. or 1.) 
    … …  
    110107      !!               bmask    : land/ocean mask at barotropic stream 
    111108      !!                          function point (=0. or 1.) and set to 0 along lateral boundaries 
    112       !!               mbathy   : number of non-zero w-levels  
     109      !!               tmask_i  : interior ocean mask 
    113110      !!---------------------------------------------------------------------- 
    114111      INTEGER  ::   ji, jj, jk      ! dummy loop indices 
    … …  
    132129      ENDIF 
    133130 
    134       IF     (      rn_shlat == 0.               ) THEN   ;   IF(lwp) WRITE(numout,*) '   ocean lateral free-slip ' 
     131      IF     (      rn_shlat == 0.               ) THEN   ;   IF(lwp) WRITE(numout,*) '   ocean lateral  free-slip ' 
    135132      ELSEIF (      rn_shlat == 2.               ) THEN   ;   IF(lwp) WRITE(numout,*) '   ocean lateral  no-slip ' 
    136133      ELSEIF ( 0. < rn_shlat .AND. rn_shlat < 2. ) THEN   ;   IF(lwp) WRITE(numout,*) '   ocean lateral  partial-slip ' 
    … …  
    145142      ! N.B. tmask has already the right boundary conditions since mbathy is ok 
    146143      ! 
    147       tmask(:,:,:) = 0.e0 
     144      tmask(:,:,:) = 0._wp 
    148145      DO jk = 1, jpk 
    149146         DO jj = 1, jpj 
    150147            DO ji = 1, jpi 
    151                IF( REAL( mbathy(ji,jj) - jk ) +.1 >= 0.e0 )   tmask(ji,jj,jk) = 1.e0 
     148               IF( REAL( mbathy(ji,jj) - jk, wp ) + 0.1_wp >= 0._wp )   tmask(ji,jj,jk) = 1._wp 
    152149            END DO   
    153150         END DO   
    … …  
    160157            ij0 =  87   ;   ij1 =  88 
    161158            ii0 = 160   ;   ii1 = 161 
    162             tmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) = 0.e0 
     159            tmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) = 0._wp 
    163160         ELSE 
    164161            IF(lwp) WRITE(numout,*) 
    … …  
    182179      ijl = nlcj - jprecj + 1 
    183180 
    184       tmask_i( 1 :iif,   :   ) = 0.e0      ! first columns 
    185       tmask_i(iil:jpi,   :   ) = 0.e0      ! last  columns (including mpp extra columns) 
    186       tmask_i(   :   , 1 :ijf) = 0.e0      ! first rows 
    187       tmask_i(   :   ,ijl:jpj) = 0.e0      ! last  rows (including mpp extra rows) 
     181      tmask_i( 1 :iif,   :   ) = 0._wp      ! first columns 
     182      tmask_i(iil:jpi,   :   ) = 0._wp      ! last  columns (including mpp extra columns) 
     183      tmask_i(   :   , 1 :ijf) = 0._wp      ! first rows 
     184      tmask_i(   :   ,ijl:jpj) = 0._wp      ! last  rows (including mpp extra rows) 
    188185 
    189186      ! north fold mask 
    190187      ! --------------- 
    191       tpol(1:jpiglo) = 1.e0  
    192       fpol(1:jpiglo) = 1.e0 
     188      tpol(1:jpiglo) = 1._wp  
     189      fpol(1:jpiglo) = 1._wp 
    193190      IF( jperio == 3 .OR. jperio == 4 ) THEN      ! T-point pivot 
    194          tpol(jpiglo/2+1:jpiglo) = 0.e0 
    195          fpol(     1    :jpiglo) = 0.e0 
     191         tpol(jpiglo/2+1:jpiglo) = 0._wp 
     192         fpol(     1    :jpiglo) = 0._wp 
    196193         IF( mjg(nlej) == jpjglo ) THEN                  ! only half of the nlcj-1 row 
    197194            DO ji = iif+1, iil-1 
    … …  
    201198      ENDIF 
    202199      IF( jperio == 5 .OR. jperio == 6 ) THEN      ! F-point pivot 
    203          tpol(     1    :jpiglo) = 0.e0 
    204          fpol(jpiglo/2+1:jpiglo) = 0.e0 
     200         tpol(     1    :jpiglo) = 0._wp 
     201         fpol(jpiglo/2+1:jpiglo) = 0._wp 
    205202      ENDIF 
    206203 
    … …  
    219216         END DO 
    220217      END DO 
    221       CALL lbc_lnk( umask, 'U', 1. )      ! Lateral boundary conditions 
    222       CALL lbc_lnk( vmask, 'V', 1. ) 
    223       CALL lbc_lnk( fmask, 'F', 1. ) 
     218      CALL lbc_lnk( umask, 'U', 1._wp )      ! Lateral boundary conditions 
     219      CALL lbc_lnk( vmask, 'V', 1._wp ) 
     220      CALL lbc_lnk( fmask, 'F', 1._wp ) 
    224221 
    225222 
    … …  
    231228      !                                    ! cyclic east-west : bmask must be set to 0. on rows 1 and jpi 
    232229      IF( nperio == 1 .OR. nperio == 4 .OR. nperio == 6 ) THEN 
    233          bmask( 1 ,:) = 0.e0 
    234          bmask(jpi,:) = 0.e0 
     230         bmask( 1 ,:) = 0._wp 
     231         bmask(jpi,:) = 0._wp 
    235232      ENDIF 
    236233      IF( nperio == 2 ) THEN               ! south symmetric :  bmask must be set to 0. on row 1 
    237          bmask(:, 1 ) = 0.e0 
     234         bmask(:, 1 ) = 0._wp 
    238235      ENDIF 
    239236      !                                    ! north fold :  
    … …  
    242239            ii = ji + nimpp - 1 
    243240            bmask(ji,jpj-1) = bmask(ji,jpj-1) * tpol(ii) 
    244             bmask(ji,jpj  ) = 0.e0 
     241            bmask(ji,jpj  ) = 0._wp 
    245242         END DO 
    246243      ENDIF 
    247244      IF( nperio == 5 .OR. nperio == 6 ) THEN   ! F-pt pivot and T-pt elliptic eq. : bmask set to 0. on row jpj 
    248          bmask(:,jpj) = 0.e0 
     245         bmask(:,jpj) = 0._wp 
    249246      ENDIF 
    250247      ! 
    251248      IF( lk_mpp ) THEN                    ! mpp specificities 
    252249         !                                      ! bmask is set to zero on the overlap region 
    253          IF( nbondi /= -1 .AND. nbondi /= 2 )   bmask(  1 :jpreci,:) = 0.e0 
    254          IF( nbondi /=  1 .AND. nbondi /= 2 )   bmask(nlci:jpi   ,:) = 0.e0 
    255          IF( nbondj /= -1 .AND. nbondj /= 2 )   bmask(:,  1 :jprecj) = 0.e0 
    256          IF( nbondj /=  1 .AND. nbondj /= 2 )   bmask(:,nlcj:jpj   ) = 0.e0 
     250         IF( nbondi /= -1 .AND. nbondi /= 2 )   bmask(  1 :jpreci,:) = 0._wp 
     251         IF( nbondi /=  1 .AND. nbondi /= 2 )   bmask(nlci:jpi   ,:) = 0._wp 
     252         IF( nbondj /= -1 .AND. nbondj /= 2 )   bmask(:,  1 :jprecj) = 0._wp 
     253         IF( nbondj /=  1 .AND. nbondj /= 2 )   bmask(:,nlcj:jpj   ) = 0._wp 
    257254         ! 
    258255         IF( npolj == 3 .OR. npolj == 4 ) THEN  ! north fold : bmask must be set to 0. on rows jpj-1 and jpj 
    … …  
    260257               ii = ji + nimpp - 1 
    261258               bmask(ji,nlcj-1) = bmask(ji,nlcj-1) * tpol(ii) 
    262                bmask(ji,nlcj  ) = 0.e0 
     259               bmask(ji,nlcj  ) = 0._wp 
    263260            END DO 
    264261         ENDIF 
    265262         IF( npolj == 5 .OR. npolj == 6 ) THEN  ! F-pt pivot and T-pt elliptic eq. : bmask set to 0. on row jpj 
    266263            DO ji = 1, nlci 
    267                bmask(ji,nlcj  ) = 0.e0 
     264               bmask(ji,nlcj  ) = 0._wp 
    268265            END DO 
    269266         ENDIF 
    … …  
    282279         DO jj = 2, jpjm1 
    283280            DO ji = fs_2, fs_jpim1   ! vector opt. 
    284                IF( fmask(ji,jj,jk) == 0. ) THEN 
    285                   fmask(ji,jj,jk) = rn_shlat * MIN( 1., MAX( zwf(ji+1,jj), zwf(ji,jj+1),   & 
    286                      &                                       zwf(ji-1,jj), zwf(ji,jj-1)  )  ) 
     281               IF( fmask(ji,jj,jk) == 0._wp ) THEN 
     282                  fmask(ji,jj,jk) = rn_shlat * MIN( 1._wp , MAX( zwf(ji+1,jj), zwf(ji,jj+1),   & 
     283                     &                                           zwf(ji-1,jj), zwf(ji,jj-1)  )  ) 
    287284               ENDIF 
    288285            END DO 
    289286         END DO 
    290287         DO jj = 2, jpjm1 
    291             IF( fmask(1,jj,jk) == 0. ) THEN 
    292                fmask(1  ,jj,jk) = rn_shlat * MIN( 1., MAX( zwf(2,jj), zwf(1,jj+1), zwf(1,jj-1) ) ) 
    293             ENDIF 
    294             IF( fmask(jpi,jj,jk) == 0. ) THEN 
    295                fmask(jpi,jj,jk) = rn_shlat * MIN( 1., MAX( zwf(jpi,jj+1), zwf(jpim1,jj), zwf(jpi,jj-1) ) ) 
     288            IF( fmask(1,jj,jk) == 0._wp ) THEN 
     289               fmask(1  ,jj,jk) = rn_shlat * MIN( 1._wp , MAX( zwf(2,jj), zwf(1,jj+1), zwf(1,jj-1) ) ) 
     290            ENDIF 
     291            IF( fmask(jpi,jj,jk) == 0._wp ) THEN 
     292               fmask(jpi,jj,jk) = rn_shlat * MIN( 1._wp , MAX( zwf(jpi,jj+1), zwf(jpim1,jj), zwf(jpi,jj-1) ) ) 
    296293            ENDIF 
    297294         END DO          
    298295         DO ji = 2, jpim1 
    299             IF( fmask(ji,1,jk) == 0. ) THEN 
    300                fmask(ji, 1 ,jk) = rn_shlat * MIN( 1., MAX( zwf(ji+1,1), zwf(ji,2), zwf(ji-1,1) ) ) 
    301             ENDIF 
    302             IF( fmask(ji,jpj,jk) == 0. ) THEN 
    303                fmask(ji,jpj,jk) = rn_shlat * MIN( 1., MAX( zwf(ji+1,jpj), zwf(ji-1,jpj), zwf(ji,jpjm1) ) ) 
     296            IF( fmask(ji,1,jk) == 0._wp ) THEN 
     297               fmask(ji, 1 ,jk) = rn_shlat * MIN( 1._wp , MAX( zwf(ji+1,1), zwf(ji,2), zwf(ji-1,1) ) ) 
     298            ENDIF 
     299            IF( fmask(ji,jpj,jk) == 0._wp ) THEN 
     300               fmask(ji,jpj,jk) = rn_shlat * MIN( 1._wp , MAX( zwf(ji+1,jpj), zwf(ji-1,jpj), zwf(ji,jpjm1) ) ) 
    304301            ENDIF 
    305302         END DO 
    … …  
    308305      IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN   ! ORCA_R2 configuration 
    309306         !                                                 ! Increased lateral friction near of some straits 
    310          IF( n_cla == 0 ) THEN 
     307         IF( nn_cla == 0 ) THEN 
    311308            !                                ! Gibraltar strait  : partial slip (fmask=0.5) 
    312309            ij0 = 101   ;   ij1 = 101 
    313             ii0 = 139   ;   ii1 = 140   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) =  0.5e0 
     310            ii0 = 139   ;   ii1 = 140   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) =  0.5_wp 
    314311            ij0 = 102   ;   ij1 = 102 
    315             ii0 = 139   ;   ii1 = 140   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) =  0.5e0 
     312            ii0 = 139   ;   ii1 = 140   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) =  0.5_wp 
    316313            ! 
    317314            !                                ! Bab el Mandeb : partial slip (fmask=1) 
    318315            ij0 =  87   ;   ij1 =  88 
    319             ii0 = 160   ;   ii1 = 160   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) =  1.e0 
     316            ii0 = 160   ;   ii1 = 160   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) =  1._wp 
    320317            ij0 =  88   ;   ij1 =  88 
    321             ii0 = 159   ;   ii1 = 159   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) =  1.e0 
     318            ii0 = 159   ;   ii1 = 159   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) =  1._wp 
    322319            ! 
    323320         ENDIF 
    324  
    325321         !                                ! Danish straits  : strong slip (fmask > 2) 
    326322! We keep this as an example but it is instable in this case  
    327323!         ij0 = 115   ;   ij1 = 115 
    328 !         ii0 = 145   ;   ii1 = 146   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) = 4.0e0 
     324!         ii0 = 145   ;   ii1 = 146   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) = 4._wp 
    329325!         ij0 = 116   ;   ij1 = 116 
    330 !         ii0 = 145   ;   ii1 = 146   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) = 4.0e0 
     326!         ii0 = 145   ;   ii1 = 146   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) = 4._wp 
    331327         ! 
    332328      ENDIF 
    333329      ! 
    334       CALL lbc_lnk( fmask, 'F', 1. )      ! Lateral boundary conditions on fmask 
    335  
    336        
    337       ! Mbathy set to the number of w-level (minimum value 2) 
    338       ! ----------------------------------- 
    339       DO jj = 1, jpj 
    340          DO ji = 1, jpi 
    341             mbathy(ji,jj) = MAX( 1, mbathy(ji,jj) ) + 1 
    342          END DO 
    343       END DO 
    344        
     330      IF( cp_cfg == "orca" .AND. jp_cfg == 1 ) THEN   ! ORCA R1 configuration 
     331         !                                                 ! Increased lateral friction near of some straits 
     332         IF(lwp) WRITE(numout,*) 
     333         IF(lwp) WRITE(numout,*) '   orca_r1: increase friction near the following straits : ' 
     334         IF(lwp) WRITE(numout,*) '      Gibraltar ' 
     335         ii0 = 283   ;   ii1 = 284        ! Gibraltar Strait  
     336         ij0 = 200   ;   ij1 = 200   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1), 1:jpk ) =  2._wp   
     337 
     338         IF(lwp) WRITE(numout,*) '      Bhosporus ' 
     339         ii0 = 314   ;   ii1 = 315        ! Bhosporus Strait  
     340         ij0 = 208   ;   ij1 = 208   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1), 1:jpk ) =  2._wp   
     341 
     342         IF(lwp) WRITE(numout,*) '      Makassar (Top) ' 
     343         ii0 =  48   ;   ii1 =  48        ! Makassar Strait (Top)  
     344         ij0 = 149   ;   ij1 = 150   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1), 1:jpk ) =  3._wp   
     345 
     346         IF(lwp) WRITE(numout,*) '      Lombok ' 
     347         ii0 =  44   ;   ii1 =  44        ! Lombok Strait  
     348         ij0 = 124   ;   ij1 = 125   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1), 1:jpk ) =  2._wp   
     349 
     350         IF(lwp) WRITE(numout,*) '      Ombai ' 
     351         ii0 =  53   ;   ii1 =  53        ! Ombai Strait  
     352         ij0 = 124   ;   ij1 = 125   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1), 1:jpk ) = 2._wp   
     353 
     354         IF(lwp) WRITE(numout,*) '      Timor Passage ' 
     355         ii0 =  56   ;   ii1 =  56        ! Timor Passage  
     356         ij0 = 124   ;   ij1 = 125   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1), 1:jpk ) = 2._wp   
     357 
     358         IF(lwp) WRITE(numout,*) '      West Halmahera ' 
     359         ii0 =  58   ;   ii1 =  58        ! West Halmahera Strait  
     360         ij0 = 141   ;   ij1 = 142   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1), 1:jpk ) = 3._wp   
     361 
     362         IF(lwp) WRITE(numout,*) '      East Halmahera ' 
     363         ii0 =  55   ;   ii1 =  55        ! East Halmahera Strait  
     364         ij0 = 141   ;   ij1 = 142   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1), 1:jpk ) = 3._wp   
     365         ! 
     366      ENDIF 
     367      ! 
     368      CALL lbc_lnk( fmask, 'F', 1._wp )      ! Lateral boundary conditions on fmask 
     369 
     370             
    345371      IF( nprint == 1 .AND. lwp ) THEN      ! Control print 
    346372         imsk(:,:) = INT( tmask_i(:,:) ) 
    … …  
    385411         imsk(:,:) = INT( bmask(:,:) ) 
    386412         CALL prihin( imsk(:,:), jpi, jpj, 1, jpi, 1,   & 
    387                &                           1, jpj, 1, 1, numout ) 
     413            &                              1, jpj, 1, 1, numout ) 
    388414      ENDIF 
    389415      ! 
    … …  
    404430      !! 
    405431      !! ** Action : 
    406       !! 
    407432      !!---------------------------------------------------------------------- 
    408433      INTEGER  :: ji, jj, jk, jl      ! dummy loop indices 
    … …  
    448473               zaa = tmask(ji  ,jj,jk) + tmask(ji  ,jj+1,jk)   & 
    449474                  &+ tmask(ji+1,jj,jk) + tmask(ji+1,jj+1,jk) 
    450                IF( ABS(zaa-3.) <= 0.1 )   fmask(ji,jj,jk) = 1. 
     475               IF( ABS(zaa-3._wp) <= 0.1_wp )   fmask(ji,jj,jk) = 1._wp 
    451476            END DO 
    452477         END DO 
    … …  
    461486            DO ji = 2, jpim1 
    462487               zaa = tmask(ji+1,jj,jk) + tmask(ji+1,jj+1,jk) 
    463                IF( ABS(zaa-2.) <= 0.1 .AND. fmask(ji,jj,jk) == 0 ) THEN 
     488               IF( ABS(zaa-2._wp) <= 0.1_wp .AND. fmask(ji,jj,jk) == 0._wp ) THEN 
    464489                  inw = inw + 1 
    465490                  nicoa(inw,1,jk) = ji 
    … …  
    468493               ENDIF 
    469494               zaa = tmask(ji,jj,jk) + tmask(ji,jj+1,jk) 
    470                IF( ABS(zaa-2.) <= 0.1 .AND. fmask(ji,jj,jk) == 0 ) THEN 
     495               IF( ABS(zaa-2._wp) <= 0.1_wp .AND. fmask(ji,jj,jk) == 0._wp ) THEN 
    471496                  ine = ine + 1 
    472497                  nicoa(ine,2,jk) = ji 
    … …  
    488513            DO ji =2, jpim1 
    489514               zaa = tmask(ji,jj+1,jk) + tmask(ji+1,jj+1,jk) 
    490                IF( ABS(zaa-2.) <= 0.1 .AND. fmask(ji,jj,jk) == 0 ) THEN 
     515               IF( ABS(zaa-2._wp) <= 0.1_wp .AND. fmask(ji,jj,jk) == 0._wp ) THEN 
    491516                  ins = ins + 1 
    492517                  nicoa(ins,3,jk) = ji 
    … …  
    495520               ENDIF 
    496521               zaa = tmask(ji+1,jj,jk) + tmask(ji,jj,jk) 
    497                IF( ABS(zaa-2.) <= 0.1 .AND. fmask(ji,jj,jk) == 0 ) THEN 
     522               IF( ABS(zaa-2._wp) <= 0.1_wp .AND. fmask(ji,jj,jk) == 0._wp ) THEN 
    498523                  inn = inn + 1 
    499524                  nicoa(inn,4,jk) = ji 
    … …  
    524549      iind = 0 
    525550      ijnd = 0 
    526       IF( nperio == 1 .OR. nperio == 4 .OR. nperio == 6 ) iind = 2 
    527       IF( nperio == 3 .OR. nperio == 4 .OR. nperio == 5 .OR. nperio == 6 ) ijnd = 2 
     551      IF( nperio == 1 .OR. nperio == 4 .OR. nperio == 6 )   iind = 2 
     552      IF( nperio == 3 .OR. nperio == 4 .OR. nperio == 5 .OR. nperio == 6 )   ijnd = 2 
    528553      DO jk = 1, jpk 
    529554         DO jl = 1, npcoa(1,jk) 
    … …  
    551576            ENDIF 
    552577         END DO 
    553          DO jl=1,npcoa(4,jk) 
     578         DO jl = 1, npcoa(4,jk) 
    554579            IF( njcoa(jl,4,jk)-2 < 1) THEN 
    555                ierror=ierror+1 
    556                icoord(ierror,1)=nicoa(jl,4,jk) 
    557                icoord(ierror,2)=njcoa(jl,4,jk) 
    558                icoord(ierror,3)=jk 
     580               ierror=ierror + 1 
     581               icoord(ierror,1) = nicoa(jl,4,jk) 
     582               icoord(ierror,2) = njcoa(jl,4,jk) 
     583               icoord(ierror,3) = jk 
    559584            ENDIF 
    560585         END DO 

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

    r1725 r2528  
    1919 
    2020   !!---------------------------------------------------------------------- 
    21    !! NEMO/OPA 3.2 , LOCEAN-IPSL (2008)  
     21   !! NEMO/OPA 3.3 , NEMO Consortium (2010) 
    2222   !! $Id$  
    23    !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 
     23   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 
    2424   !!---------------------------------------------------------------------- 
    2525 

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

    • Property svn:eol-style deleted
    r1152 r2528  
    2727#  include "domzgr_substitute.h90" 
    2828   !!---------------------------------------------------------------------- 
    29    !!   OPA 9.0 , LOCEAN-IPSL (2005)  
     29   !! NEMO/OPA 3.3 , NEMO Consortium (2010) 
    3030   !! $Id$  
    31    !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 
     31   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 
    3232   !!---------------------------------------------------------------------- 
    3333 

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

    • Property svn:eol-style deleted
    r1983 r2528  
    3636#  include "vectopt_loop_substitute.h90" 
    3737   !!---------------------------------------------------------------------- 
    38    !! NEMO/OPA 3.2 , LOCEAN-IPSL (2009)  
     38   !! NEMO/OPA 3.3 , NEMO Consortium (2010) 
    3939   !! $Id$ 
    40    !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 
     40   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 
    4141   !!---------------------------------------------------------------------- 
    4242 
    … …  
    5151      !!---------------------------------------------------------------------- 
    5252      INTEGER  ::   ji, jj, jk 
    53       REAL(wp) ::   zcoefu, zcoefv, zcoeff 
     53      REAL(wp) ::   zcoefu , zcoefv   , zcoeff                   ! temporary scalars 
     54      REAL(wp) ::   zv_t_ij, zv_t_ip1j, zv_t_ijp1, zv_t_ip1jp1   !     -        - 
     55      REAL(wp), DIMENSION(jpi,jpj) ::  zs_t, zs_u_1, zs_v_1      !     -     2D workspace 
    5456      !!---------------------------------------------------------------------- 
    5557 
    … …  
    6062      ENDIF 
    6163 
    62       IF( lk_zco )   CALL ctl_stop( 'dom_vvl : key_zco is incompatible with variable volume option key_vvl') 
    6364 
    64       IF( ln_zco) THEN 
    65          DO jk = 1, jpk 
    66             gdept(:,:,jk) = gdept_0(jk) 
    67             gdepw(:,:,jk) = gdepw_0(jk) 
    68             gdep3w(:,:,jk) = gdepw_0(jk) 
    69             e3t (:,:,jk) = e3t_0(jk) 
    70             e3u (:,:,jk) = e3t_0(jk) 
    71             e3v (:,:,jk) = e3t_0(jk) 
    72             e3f (:,:,jk) = e3t_0(jk) 
    73             e3w (:,:,jk) = e3w_0(jk) 
    74             e3uw(:,:,jk) = e3w_0(jk) 
    75             e3vw(:,:,jk) = e3w_0(jk) 
    76          END DO 
    77       ELSE 
    78          fsdept(:,:,:) = gdept (:,:,:) 
    79          fsdepw(:,:,:) = gdepw (:,:,:) 
    80          fsde3w(:,:,:) = gdep3w(:,:,:) 
    81          fse3t (:,:,:) = e3t   (:,:,:) 
    82          fse3u (:,:,:) = e3u   (:,:,:) 
    83          fse3v (:,:,:) = e3v   (:,:,:) 
    84          fse3f (:,:,:) = e3f   (:,:,:) 
    85          fse3w (:,:,:) = e3w   (:,:,:) 
    86          fse3uw(:,:,:) = e3uw  (:,:,:) 
    87          fse3vw(:,:,:) = e3vw  (:,:,:) 
    88       ENDIF 
     65      fsdept(:,:,:) = gdept (:,:,:) 
     66      fsdepw(:,:,:) = gdepw (:,:,:) 
     67      fsde3w(:,:,:) = gdep3w(:,:,:) 
     68      fse3t (:,:,:) = e3t   (:,:,:) 
     69      fse3u (:,:,:) = e3u   (:,:,:) 
     70      fse3v (:,:,:) = e3v   (:,:,:) 
     71      fse3f (:,:,:) = e3f   (:,:,:) 
     72      fse3w (:,:,:) = e3w   (:,:,:) 
     73      fse3uw(:,:,:) = e3uw  (:,:,:) 
     74      fse3vw(:,:,:) = e3vw  (:,:,:) 
    8975 
    9076      !                                 !==  mu computation  ==! 
    … …  
    130116         hv_0(:,:) = hv_0(:,:) + fse3v_0(:,:,jk) * vmask(:,:,jk) 
    131117      END DO 
     118       
     119      ! surface at t-points and inverse surface at (u/v)-points used in surface averaging computations 
     120      ! for ssh and scale factors 
     121      zs_t  (:,:) =       e1t(:,:) * e2t(:,:) 
     122      zs_u_1(:,:) = 0.5 / e1u(:,:) * e2u(:,:) 
     123      zs_v_1(:,:) = 0.5 / e1v(:,:) * e2v(:,:) 
    132124 
    133125      DO jj = 1, jpjm1                          ! initialise before and now Sea Surface Height at u-, v-, f-points 
    134126         DO ji = 1, jpim1   ! NO vector opt. 
    135             zcoefu = 0.5  * umask(ji,jj,1) / ( e1u(ji,jj) * e2u(ji,jj) ) 
    136             zcoefv = 0.5  * vmask(ji,jj,1) / ( e1v(ji,jj) * e2v(ji,jj) ) 
    137             zcoeff = 0.25 * umask(ji,jj,1) * umask(ji,jj+1,1) 
    138             sshu_b(ji,jj) = zcoefu * ( e1t(ji  ,jj) * e2t(ji  ,jj) * sshb(ji  ,jj)     & 
    139                &                     + e1t(ji+1,jj) * e2t(ji+1,jj) * sshb(ji+1,jj) )    
    140             sshv_b(ji,jj) = zcoefv * ( e1t(ji,jj  ) * e2t(ji,jj  ) * sshb(ji,jj  )     & 
    141                &                     + e1t(ji,jj+1) * e2t(ji,jj+1) * sshb(ji,jj+1) )    
    142             sshf_b(ji,jj) = zcoeff * ( sshb(ji  ,jj) + sshb(ji  ,jj+1)                 & 
    143                &                     + sshb(ji+1,jj) + sshb(ji+1,jj+1) )                
    144             sshu_n(ji,jj) = zcoefu * ( e1t(ji  ,jj) * e2t(ji  ,jj) * sshn(ji  ,jj)     & 
    145                &                     + e1t(ji+1,jj) * e2t(ji+1,jj) * sshn(ji+1,jj) )    
    146             sshv_n(ji,jj) = zcoefv * ( e1t(ji,jj  ) * e2t(ji,jj  ) * sshn(ji,jj  )     &  
    147                &                     + e1t(ji,jj+1) * e2t(ji,jj+1) * sshn(ji,jj+1) )      
    148             sshf_n(ji,jj) = zcoeff * ( sshn(ji  ,jj) + sshn(ji  ,jj+1)                 & 
    149                &                     + sshn(ji+1,jj) + sshn(ji+1,jj+1) )                
     127            zcoefu = umask(ji,jj,1) * zs_u_1(ji,jj) 
     128            zcoefv = vmask(ji,jj,1) * zs_v_1(ji,jj) 
     129            zcoeff = 0.5 * umask(ji,jj,1) * umask(ji,jj+1,1) / ( e1f(ji,jj) * e2f(ji,jj) ) 
     130            ! before fields 
     131            zv_t_ij       = zs_t(ji  ,jj  ) * sshb(ji  ,jj  ) 
     132            zv_t_ip1j     = zs_t(ji+1,jj  ) * sshb(ji+1,jj  ) 
     133            zv_t_ijp1     = zs_t(ji  ,jj+1) * sshb(ji  ,jj+1) 
     134            sshu_b(ji,jj) = zcoefu * ( zv_t_ij + zv_t_ip1j ) 
     135            sshv_b(ji,jj) = zcoefv * ( zv_t_ij + zv_t_ijp1 ) 
     136            ! now fields 
     137            zv_t_ij       = zs_t(ji  ,jj  ) * sshn(ji  ,jj  ) 
     138            zv_t_ip1j     = zs_t(ji+1,jj  ) * sshn(ji+1,jj  ) 
     139            zv_t_ijp1     = zs_t(ji  ,jj+1) * sshn(ji  ,jj+1) 
     140            zv_t_ip1jp1   = zs_t(ji  ,jj+1) * sshn(ji  ,jj+1) 
     141            sshu_n(ji,jj) = zcoefu * ( zv_t_ij + zv_t_ip1j ) 
     142            sshv_n(ji,jj) = zcoefv * ( zv_t_ij + zv_t_ijp1 ) 
     143            sshf_n(ji,jj) = zcoeff * ( zv_t_ij + zv_t_ip1j + zv_t_ijp1 + zv_t_ip1jp1 ) 
    150144         END DO 
    151145      END DO 
    152       CALL lbc_lnk( sshu_b, 'U', 1. )   ;   CALL lbc_lnk( sshu_n, 'U', 1. )      ! lateral boundary conditions 
    153       CALL lbc_lnk( sshv_b, 'V', 1. )   ;   CALL lbc_lnk( sshv_n, 'V', 1. ) 
    154       CALL lbc_lnk( sshf_b, 'F', 1. )   ;   CALL lbc_lnk( sshf_n, 'F', 1. ) 
     146      CALL lbc_lnk( sshu_n, 'U', 1. )   ;   CALL lbc_lnk( sshu_b, 'U', 1. )      ! lateral boundary conditions 
     147      CALL lbc_lnk( sshv_n, 'V', 1. )   ;   CALL lbc_lnk( sshv_b, 'V', 1. ) 
     148      CALL lbc_lnk( sshf_n, 'F', 1. ) 
     149 
     150                                                ! initialise before scale factors at (u/v)-points 
     151      ! Scale factor anomaly at (u/v)-points: surface averaging of scale factor at t-points 
     152      DO jk = 1, jpkm1 
     153         DO jj = 1, jpjm1 
     154            DO ji = 1, jpim1 
     155               zv_t_ij           = zs_t(ji  ,jj  ) * fse3t_b(ji  ,jj  ,jk) 
     156               zv_t_ip1j         = zs_t(ji+1,jj  ) * fse3t_b(ji+1,jj  ,jk) 
     157               zv_t_ijp1         = zs_t(ji  ,jj+1) * fse3t_b(ji  ,jj+1,jk) 
     158               fse3u_b(ji,jj,jk) = umask(ji,jj,jk) * ( zs_u_1(ji,jj) * ( zv_t_ij + zv_t_ip1j ) - fse3u_0(ji,jj,jk) ) 
     159               fse3v_b(ji,jj,jk) = vmask(ji,jj,jk) * ( zs_v_1(ji,jj) * ( zv_t_ij + zv_t_ijp1 ) - fse3v_0(ji,jj,jk) ) 
     160            END DO 
     161         END DO 
     162      END DO 
     163      CALL lbc_lnk( fse3u_b(:,:,:), 'U', 1. )               ! lateral boundary conditions 
     164      CALL lbc_lnk( fse3v_b(:,:,:), 'V', 1. ) 
     165      ! Add initial scale factor to scale factor anomaly 
     166      fse3u_b(:,:,:) = fse3u_b(:,:,:) + fse3u_0(:,:,:) 
     167      fse3v_b(:,:,:) = fse3v_b(:,:,:) + fse3v_0(:,:,:) 
    155168      ! 
    156          DO jk = 1, jpkm1 
    157             fsdept(:,:,jk) = fsdept_n(:,:,jk)          ! now local depths stored in fsdep. arrays 
    158             fsdepw(:,:,jk) = fsdepw_n(:,:,jk) 
    159             fsde3w(:,:,jk) = fsde3w_n(:,:,jk) 
    160             ! 
    161             fse3t (:,:,jk) = fse3t_n (:,:,jk)          ! vertical scale factors stored in fse3. arrays 
    162             fse3u (:,:,jk) = fse3u_n (:,:,jk) 
    163             fse3v (:,:,jk) = fse3v_n (:,:,jk) 
    164             fse3f (:,:,jk) = fse3f_n (:,:,jk) 
    165             fse3w (:,:,jk) = fse3w_n (:,:,jk) 
    166             fse3uw(:,:,jk) = fse3uw_n(:,:,jk) 
    167             fse3vw(:,:,jk) = fse3vw_n(:,:,jk) 
    168          END DO 
    169  
    170  
    171  
    172169   END SUBROUTINE dom_vvl 
    173170 

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

    • Property svn:eol-style deleted
    r1929 r2528  
    44   !! Ocean initialization : write the ocean domain mesh ask file(s) 
    55   !!====================================================================== 
     6   !! History :  OPA  ! 1997-02  (G. Madec)  Original code 
     7   !!            8.1  ! 1999-11  (M. Imbard)  NetCDF FORMAT with IOIPSL 
     8   !!   NEMO     1.0  ! 2002-08  (G. Madec)  F90 and several file 
     9   !!---------------------------------------------------------------------- 
    610 
    711   !!---------------------------------------------------------------------- 
    … …  
    1115   !!                         = 3  :   mesh_hgr, mesh_zgr and mask 
    1216   !!---------------------------------------------------------------------- 
    13    !! * Modules used 
    1417   USE dom_oce         ! ocean space and time domain 
    15    USE in_out_manager 
    16    USE iom 
    17    USE lbclnk          ! ocean lateral boundary conditions (or mpp link) 
    18    USE lib_mpp 
     18   USE in_out_manager  ! I/O manager 
     19   USE iom             ! I/O library 
     20   USE lbclnk          ! lateral boundary conditions - mpp exchanges 
     21   USE lib_mpp         ! MPP library 
    1922 
    2023   IMPLICIT NONE 
    … …  
    2629#  include "vectopt_loop_substitute.h90" 
    2730   !!---------------------------------------------------------------------- 
    28    !!   OPA 9.0 , LOCEAN-IPSL (2005)  
     31   !! NEMO/OPA 3.3 , NEMO Consortium (2010) 
    2932   !! $Id$  
    30    !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt  
    31    !!---------------------------------------------------------------------- 
    32  
     33   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
     34   !!---------------------------------------------------------------------- 
    3335CONTAINS 
    3436 
    … …  
    5456      !!      if     nmsh <= 3: write full 3D arrays for e3[tuvw] and gdep[tuvw] 
    5557      !!      if 3 < nmsh <= 6: write full 3D arrays for e3[tuvw] and 2D arrays  
    56       !!                        corresponding to the depth of the bottom points hdep[tw] 
     58      !!                        corresponding to the depth of the bottom t- and w-points 
    5759      !!      if 6 < nmsh <= 9: write 2D arrays corresponding to the depth and the 
    58       !!                        thickness of the bottom points hdep[tw] and e3[tw]_ps 
    59       !! 
    60       !! ** output file :  
    61       !!      meshmask.nc  : domain size, horizontal grid-point position, 
    62       !!                     masks, depth and vertical scale factors 
    63       !! 
    64       !! History : 
    65       !!        !  97-02  (G. Madec)  Original code 
    66       !!        !  99-11  (M. Imbard)  NetCDF FORMAT with IOIPSL 
    67       !!   9.0  !  02-08  (G. Madec)  F90 and several file 
    68       !!---------------------------------------------------------------------- 
    69       INTEGER                          ::   inum0    ! temprary units for 'mesh_mask.nc' file 
    70       INTEGER                          ::   inum1    ! temprary units for 'mesh.nc'      file 
    71       INTEGER                          ::   inum2    ! temprary units for 'mask.nc'      file 
    72       INTEGER                          ::   inum3    ! temprary units for 'mesh_hgr.nc'  file 
    73       INTEGER                          ::   inum4    ! temprary units for 'mesh_zgr.nc'  file 
    74       INTEGER                          ::   ji, jj, jk, ik 
    75       REAL(wp), DIMENSION(jpi,jpj)     ::   zprt     ! temporary array for bathymetry  
    76       REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zdepu    ! 3D depth of U point 
    77       REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zdepv    ! 3D depth of V point 
    78       CHARACTER(len=21)                ::   clnam0   ! filename (mesh and mask informations) 
    79       CHARACTER(len=21)                ::   clnam1   ! filename (mesh informations) 
    80       CHARACTER(len=21)                ::   clnam2   ! filename (mask informations) 
    81       CHARACTER(len=21)                ::   clnam3   ! filename (horizontal mesh informations) 
    82       CHARACTER(len=21)                ::   clnam4   ! filename (vertical   mesh informations) 
    83       !!---------------------------------------------------------------------- 
     60      !!                        thickness (e3[tw]_ps) of the bottom points  
     61      !! 
     62      !! ** output file :   meshmask.nc  : domain size, horizontal grid-point position, 
     63      !!                                   masks, depth and vertical scale factors 
     64      !!---------------------------------------------------------------------- 
     65      INTEGER           ::   inum0    ! temprary units for 'mesh_mask.nc' file 
     66      INTEGER           ::   inum1    ! temprary units for 'mesh.nc'      file 
     67      INTEGER           ::   inum2    ! temprary units for 'mask.nc'      file 
     68      INTEGER           ::   inum3    ! temprary units for 'mesh_hgr.nc'  file 
     69      INTEGER           ::   inum4    ! temprary units for 'mesh_zgr.nc'  file 
     70      CHARACTER(len=21) ::   clnam0   ! filename (mesh and mask informations) 
     71      CHARACTER(len=21) ::   clnam1   ! filename (mesh informations) 
     72      CHARACTER(len=21) ::   clnam2   ! filename (mask informations) 
     73      CHARACTER(len=21) ::   clnam3   ! filename (horizontal mesh informations) 
     74      CHARACTER(len=21) ::   clnam4   ! filename (vertical   mesh informations) 
     75      INTEGER           ::   ji, jj, jk   ! dummy loop indices 
     76      REAL(wp), DIMENSION(jpi,jpj)     ::   zprt , zprw    ! 2D workspace 
     77      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zdepu, zdepv   ! 3D workspace 
     78     !!---------------------------------------------------------------------- 
    8479 
    8580      IF(lwp) WRITE(numout,*) 
    … …  
    118113         CALL iom_open( TRIM(clnam3), inum3, ldwrt = .TRUE., kiolib = jprstlib ) 
    119114         CALL iom_open( TRIM(clnam4), inum4, ldwrt = .TRUE., kiolib = jprstlib ) 
    120           
     115         ! 
    121116      END SELECT 
    122117       
    … …  
    160155      CALL iom_rstput( 0, 0, inum3, 'ff', ff, ktype = jp_r8 )           !    ! coriolis factor 
    161156       
    162 !       note that mbathy has been modified in dommsk or in solver. 
    163 !       it is the number of non-zero "w" levels in the water, and the minimum  
    164 !       value (on land) is 2. We define zprt as the number of "T" points in the ocean  
    165 !       at any location, and zero on land.  
    166 ! 
    167       zprt = tmask(:,:,1)*(mbathy-1) 
    168       CALL iom_rstput( 0, 0, inum4, 'mbathy', zprt, ktype = jp_i2 ) 
     157      ! note that mbkt is set to 1 over land ==> use surface tmask 
     158      zprt(:,:) = tmask(:,:,1) * REAL( mbkt(:,:) , wp ) 
     159      CALL iom_rstput( 0, 0, inum4, 'mbathy', zprt, ktype = jp_i2 )     !    ! nb of ocean T-points 
    169160             
    170 #if ! defined key_zco 
    171161      IF( ln_sco ) THEN                                         ! s-coordinate 
    172162         CALL iom_rstput( 0, 0, inum4, 'hbatt', hbatt )         !    ! depth 
    … …  
    174164         CALL iom_rstput( 0, 0, inum4, 'hbatv', hbatv ) 
    175165         CALL iom_rstput( 0, 0, inum4, 'hbatf', hbatf ) 
    176           
     166         ! 
    177167         CALL iom_rstput( 0, 0, inum4, 'gsigt', gsigt )         !    ! scaling coef. 
    178168         CALL iom_rstput( 0, 0, inum4, 'gsigw', gsigw )   
    … …  
    180170         CALL iom_rstput( 0, 0, inum4, 'esigt', esigt ) 
    181171         CALL iom_rstput( 0, 0, inum4, 'esigw', esigw ) 
    182           
     172         ! 
    183173         CALL iom_rstput( 0, 0, inum4, 'e3t', e3t )             !    ! scale factors 
    184174         CALL iom_rstput( 0, 0, inum4, 'e3u', e3u ) 
    185175         CALL iom_rstput( 0, 0, inum4, 'e3v', e3v ) 
    186176         CALL iom_rstput( 0, 0, inum4, 'e3w', e3w ) 
    187           
     177         ! 
    188178         CALL iom_rstput( 0, 0, inum4, 'gdept_0' , gdept_0 )    !    ! stretched system 
    189179         CALL iom_rstput( 0, 0, inum4, 'gdepw_0' , gdepw_0 ) 
    … …  
    191181       
    192182      IF( ln_zps ) THEN                                         ! z-coordinate - partial steps 
    193  
     183         ! 
    194184         IF( nmsh <= 6 ) THEN                                   !    ! 3D vertical scale factors 
    195185            CALL iom_rstput( 0, 0, inum4, 'e3t', e3t )          
    … …  
    197187            CALL iom_rstput( 0, 0, inum4, 'e3v', e3v ) 
    198188            CALL iom_rstput( 0, 0, inum4, 'e3w', e3w ) 
    199          ELSE                                                   !    ! 2D bottom scale factors 
    200             DO jj = 1,jpj   ;   DO ji = 1,jpi 
    201                ik = NINT( zprt(ji,jj) )   ! take care that mbathy is not what you think it is here ! 
    202                IF ( ik /= 0 ) THEN   ;   e3tp(ji,jj) = e3t(ji,jj,ik)   ;   e3wp(ji,jj) = e3w(ji,jj,ik) 
    203                ELSE                  ;   e3tp(ji,jj) = 0.              ;   e3wp(ji,jj) = 0. 
    204                ENDIF 
    205             END DO   ;   END DO 
     189         ELSE                                                   !    ! 2D masked bottom ocean scale factors 
     190            DO jj = 1,jpj    
     191               DO ji = 1,jpi 
     192                  e3tp(ji,jj) = e3t(ji,jj,mbkt(ji,jj)) * tmask(ji,jj,1) 
     193                  e3wp(ji,jj) = e3w(ji,jj,mbkt(ji,jj)) * tmask(ji,jj,1) 
     194               END DO 
     195            END DO 
    206196            CALL iom_rstput( 0, 0, inum4, 'e3t_ps', e3tp )       
    207197            CALL iom_rstput( 0, 0, inum4, 'e3w_ps', e3wp ) 
    208198         END IF 
    209  
     199         ! 
    210200         IF( nmsh <= 3 ) THEN                                   !    ! 3D depth 
    211201            CALL iom_rstput( 0, 0, inum4, 'gdept', gdept, ktype = jp_r4 )      
    212             DO jk = 1,jpk   ;   DO jj = 1, jpjm1   ;   DO ji = 1, fs_jpim1   ! vector opt. 
    213                zdepu(ji,jj,jk) = MIN( gdept(ji,jj,jk), gdept(ji+1,jj  ,jk) ) 
    214                zdepv(ji,jj,jk) = MIN( gdept(ji,jj,jk), gdept(ji  ,jj+1,jk) ) 
    215             END DO   ;   END DO   ;   END DO 
     202            DO jk = 1,jpk    
     203               DO jj = 1, jpjm1    
     204                  DO ji = 1, fs_jpim1   ! vector opt. 
     205                     zdepu(ji,jj,jk) = MIN( gdept(ji,jj,jk) , gdept(ji+1,jj  ,jk) ) 
     206                     zdepv(ji,jj,jk) = MIN( gdept(ji,jj,jk) , gdept(ji  ,jj+1,jk) ) 
     207                  END DO    
     208               END DO    
     209            END DO 
    216210            CALL lbc_lnk( zdepu, 'U', 1. )   ;   CALL lbc_lnk( zdepv, 'V', 1. )  
    217211            CALL iom_rstput( 0, 0, inum4, 'gdepu', zdepu, ktype = jp_r4 ) 
    … …  
    219213            CALL iom_rstput( 0, 0, inum4, 'gdepw', gdepw, ktype = jp_r4 ) 
    220214         ELSE                                                   !    ! 2D bottom depth 
    221             DO jj = 1,jpj   ;   DO ji = 1,jpi 
    222                ik = NINT( zprt(ji,jj) )   ! take care that mbathy is not what you think it is here ! 
    223                IF ( ik /= 0 ) THEN   ;   hdept(ji,jj) = gdept(ji,jj,ik)   ;   hdepw(ji,jj) = gdepw(ji,jj,ik+1) 
    224                ELSE                  ;   hdept(ji,jj) = 0.                ;   hdepw(ji,jj) = 0. 
    225                ENDIF 
    226             END DO   ;   END DO 
    227             CALL iom_rstput( 0, 0, inum4, 'hdept' , hdept, ktype = jp_r4 )      
    228             CALL iom_rstput( 0, 0, inum4, 'hdepw' , hdepw, ktype = jp_r4 )  
     215            DO jj = 1,jpj    
     216               DO ji = 1,jpi 
     217                  zprt(ji,jj) = gdept(ji,jj,mbkt(ji,jj)  ) * tmask(ji,jj,1) 
     218                  zprw(ji,jj) = gdepw(ji,jj,mbkt(ji,jj)+1) * tmask(ji,jj,1) 
     219               END DO 
     220            END DO 
     221            CALL iom_rstput( 0, 0, inum4, 'hdept', zprt, ktype = jp_r4 )      
     222            CALL iom_rstput( 0, 0, inum4, 'hdepw', zprw, ktype = jp_r4 )  
    229223         ENDIF 
    230  
     224         ! 
    231225         CALL iom_rstput( 0, 0, inum4, 'gdept_0', gdept_0 )     !    ! reference z-coord. 
    232226         CALL iom_rstput( 0, 0, inum4, 'gdepw_0', gdepw_0 ) 
    … …  
    234228         CALL iom_rstput( 0, 0, inum4, 'e3w_0'  , e3w_0   ) 
    235229      ENDIF 
    236        
    237 #endif 
    238230       
    239231      IF( ln_zco ) THEN 
    … …  
    258250         CALL iom_close( inum4 ) 
    259251      END SELECT 
    260        
     252      ! 
    261253   END SUBROUTINE dom_wri 
    262254 
    … …  
    270262      !! ** Method  :   1) aplly lbc_lnk on an array with different values for each element 
    271263      !!                2) check which elements have been changed 
    272       !! 
    273264      !!---------------------------------------------------------------------- 
    274265      CHARACTER(len=1)            , INTENT(in   ) ::  cdgrd   !  
    275266      REAL(wp), DIMENSION(jpi,jpj)                ::  puniq   !  
    276  
     267      ! 
    277268      REAL(wp), DIMENSION(jpi,jpj  ) ::  ztstref   ! array with different values for each element  
    278269      REAL(wp)                       ::  zshift    ! shift value link to the process number 
    … …  
    280271      INTEGER                        ::  ji        ! dummy loop indices 
    281272      !!---------------------------------------------------------------------- 
    282  
     273      ! 
    283274      ! build an array with different values for each element  
    284275      ! in mpp: make sure that these values are different even between process 
    285276      ! -> apply a shift value according to the process number 
    286277      zshift = jpi * jpj * ( narea - 1 ) 
    287       ztstref(:,:) = RESHAPE( (/ (zshift + REAL(ji, wp), ji = 1, jpi*jpj) /), (/ jpi, jpj /) ) 
    288     
     278      ztstref(:,:) = RESHAPE( (/ (zshift + REAL(ji,wp), ji = 1, jpi*jpj) /), (/ jpi, jpj /) ) 
     279      ! 
    289280      puniq(:,:) = ztstref(:,:)                   ! default definition 
    290281      CALL lbc_lnk( puniq, cdgrd, 1. )            ! apply boundary conditions 
    291282      lldbl(:,:,1) = puniq(:,:) == ztstref(:,:)   ! check which values have been changed  
    292        
     283      ! 
    293284      puniq(:,:) = 1.                             ! default definition 
    294285      ! fill only the inner part of the cpu with llbl converted into real  
    295       puniq(nldi:nlei,nldj:nlej) = REAL(COUNT( lldbl(nldi:nlei,nldj:nlej,:), dim = 3 ), wp) 
    296  
     286      puniq(nldi:nlei,nldj:nlej) = REAL( COUNT( lldbl(nldi:nlei,nldj:nlej,:), dim = 3 ) , wp ) 
     287      ! 
    297288   END FUNCTION dom_uniq 
    298  
    299289 
    300290   !!====================================================================== 

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

    • Property svn:eol-style deleted
    r2465 r2528  
    77   !!                 ! 1997-07  (G. Madec)  lbc_lnk call 
    88   !!                 ! 1997-04  (J.-O. Beismann)  
    9    !!            8.5  ! 2002-09 (A. Bozec, G. Madec)  F90: Free form and module 
    10    !!             -   ! 2002-09 (A. de Miranda)  rigid-lid + islands 
     9   !!            8.5  ! 2002-09  (A. Bozec, G. Madec)  F90: Free form and module 
     10   !!             -   ! 2002-09  (A. de Miranda)  rigid-lid + islands 
    1111   !!  NEMO      1.0  ! 2003-08  (G. Madec)  F90: Free form and module 
    1212   !!             -   ! 2005-10  (A. Beckmann)  modifications for hybrid s-ccordinates & new stretching function 
    … …  
    1414   !!            3.0  ! 2008-06  (G. Madec)  insertion of domzgr_zps.h90 & conding style 
    1515   !!            3.2  ! 2009-07  (R. Benshila) Suppression of rigid-lid option 
     16   !!            3.3  ! 2010-11  (G. Madec) add mbk. arrays associated to the deepest ocean level 
    1617   !!---------------------------------------------------------------------- 
    1718 
    … …  
    2122   !!       zgr_bat_zoom : modify the bathymetry field if zoom domain 
    2223   !!       zgr_bat_ctl  : check the bathymetry files 
     24   !!       zgr_bot_level: deepest ocean level for t-, u, and v-points 
    2325   !!       zgr_z        : reference z-coordinate  
    2426   !!       zgr_zco      : z-coordinate  
    … …  
    2830   !!       dfssig       : derivative of the sigma coordinate function    !!gm  (currently missing!) 
    2931   !!--------------------------------------------------------------------- 
    30    USE oce             ! ocean dynamics and tracers 
    31    USE dom_oce         ! ocean space and time domain 
    32    USE in_out_manager  ! I/O manager 
    33    USE lib_mpp         ! distributed memory computing library 
    34    USE lbclnk          ! ocean lateral boundary conditions (or mpp link) 
    35    USE closea          ! closed seas 
    36    USE c1d             ! 1D configuration 
     32   USE oce               ! ocean variables 
     33   USE dom_oce           ! ocean domain 
     34   USE closea            ! closed seas 
     35   USE c1d               ! 1D vertical configuration 
     36   USE in_out_manager    ! I/O manager 
     37   USE iom               ! I/O library 
     38   USE lbclnk            ! ocean lateral boundary conditions (or mpp link) 
     39   USE lib_mpp           ! distributed memory computing library 
    3740 
    3841   IMPLICIT NONE 
    3942   PRIVATE 
    4043 
    41    PUBLIC   dom_zgr    ! called by dom_init.F90 
    42  
    43 !!gm   DOCTOR name for the namelist parameter... 
    44    !                                    !!! ** Namelist namzgr_sco ** 
    45    REAL(wp) ::   rn_sbot_min =  300.     ! minimum depth of s-bottom surface (>0) (m) 
    46    REAL(wp) ::   rn_sbot_max = 5250.     ! maximum depth of s-bottom surface (= ocean depth) (>0) (m) 
    47    REAL(wp) ::   rn_theta    =    6.0    ! surface control parameter (0<=rn_theta<=20) 
    48    REAL(wp) ::   rn_thetb    =    0.75   ! bottom control parameter  (0<=rn_thetb<= 1) 
    49    REAL(wp) ::   rn_rmax     =    0.15   ! maximum cut-off r-value allowed (0<rn_rmax<1) 
    50    LOGICAL  ::   ln_s_sigma  = .false.   ! use hybrid s-sigma -coordinate & stretching function fssig1 (ln_sco=T) 
    51    REAL(wp) ::   rn_bb       =    0.8    ! stretching parameter for song and haidvogel stretching 
    52    !                                     ! ( rn_bb=0; top only, rn_bb =1; top and bottom) 
    53    REAL(wp) ::   rn_hc       = 150.      ! Critical depth for s-sigma coordinates 
     44   PUBLIC   dom_zgr      ! called by dom_init.F90 
     45 
     46   !                                       !!* Namelist namzgr_sco * 
     47   REAL(wp) ::   rn_sbot_min =  300._wp     ! minimum depth of s-bottom surface (>0) (m) 
     48   REAL(wp) ::   rn_sbot_max = 5250._wp     ! maximum depth of s-bottom surface (= ocean depth) (>0) (m) 
     49   REAL(wp) ::   rn_theta    =    6.00_wp   ! surface control parameter (0<=rn_theta<=20) 
     50   REAL(wp) ::   rn_thetb    =    0.75_wp   ! bottom control parameter  (0<=rn_thetb<= 1) 
     51   REAL(wp) ::   rn_rmax     =    0.15_wp   ! maximum cut-off r-value allowed (0<rn_rmax<1) 
     52   LOGICAL  ::   ln_s_sigma  = .false.      ! use hybrid s-sigma -coordinate & stretching function fssig1 (ln_sco=T) 
     53   REAL(wp) ::   rn_bb       =    0.80_wp   ! stretching parameter for song and haidvogel stretching 
     54   !                                        ! ( rn_bb=0; top only, rn_bb =1; top and bottom) 
     55   REAL(wp) ::   rn_hc       =  150._wp     ! Critical depth for s-sigma coordinates 
    5456  
    5557   !! * Substitutions 
    … …  
    5759#  include "vectopt_loop_substitute.h90" 
    5860   !!---------------------------------------------------------------------- 
    59    !! NEMO/OPA 3.0 , LOCEAN-IPSL (2008)  
     61   !! NEMO/OPA 3.3 , NEMO Consortium (2010) 
    6062   !! $Id$ 
    61    !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 
     63   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
    6264   !!---------------------------------------------------------------------- 
    63  
    6465CONTAINS        
    6566 
    … …  
    7576      !!              - vertical coordinate (gdep., e3.) depending on the  
    7677      !!                coordinate chosen : 
    77       !!                   ln_zco=T   z-coordinate   (forced if lk_zco) 
     78      !!                   ln_zco=T   z-coordinate    
    7879      !!                   ln_zps=T   z-coordinate with partial steps 
    7980      !!                   ln_zco=T   s-coordinate  
    … …  
    8283      !!---------------------------------------------------------------------- 
    8384      INTEGER ::   ioptio = 0   ! temporary integer 
    84       !! 
     85      ! 
    8586      NAMELIST/namzgr/ ln_zco, ln_zps, ln_sco 
    8687      !!---------------------------------------------------------------------- 
    8788 
    88       REWIND ( numnam )                ! Read Namelist namzgr : vertical coordinate' 
    89       READ   ( numnam, namzgr ) 
     89      REWIND( numnam )                 ! Read Namelist namzgr : vertical coordinate' 
     90      READ  ( numnam, namzgr ) 
    9091 
    9192      IF(lwp) THEN                     ! Control print 
    … …  
    103104      IF( ln_zps ) ioptio = ioptio + 1 
    104105      IF( ln_sco ) ioptio = ioptio + 1 
    105       IF ( ioptio /= 1 )   CALL ctl_stop( ' none or several vertical coordinate options used' ) 
    106       IF( lk_zco ) THEN 
    107           IF(lwp) WRITE(numout,*) '          z-coordinate with reduced incore memory requirement' 
    108           IF( ln_zps .OR. ln_sco )   CALL ctl_stop( ' reduced memory with zps or sco option is impossible' ) 
    109       ENDIF 
    110  
     106      IF( ioptio /= 1 )   CALL ctl_stop( ' none or several vertical coordinate options used' ) 
     107      ! 
    111108      ! Build the vertical coordinate system 
    112109      ! ------------------------------------ 
    113                      CALL zgr_z        ! Reference z-coordinate system (always called) 
    114                      CALL zgr_bat      ! Bathymetry fields (levels and meters) 
    115       IF( ln_zco )   CALL zgr_zco      ! z-coordinate 
    116       IF( ln_zps )   CALL zgr_zps      ! Partial step z-coordinate 
    117       IF( ln_sco )   CALL zgr_sco      ! s-coordinate or hybrid z-s coordinate 
    118       ! 
    119      ! final adjustment of mbathy & check 
    120      ! ----------------------------------- 
    121      IF( lzoom       )   CALL zgr_bat_zoom     ! correct mbathy in case of zoom subdomain 
    122      IF( .NOT.lk_c1d )   CALL zgr_bat_ctl      ! check bathymetry (mbathy) and suppress isoated ocean points 
    123  
    124  
    125 !!bug gm control print: 
     110                          CALL zgr_z            ! Reference z-coordinate system (always called) 
     111                          CALL zgr_bat          ! Bathymetry fields (levels and meters) 
     112      IF( ln_zco      )   CALL zgr_zco          ! z-coordinate 
     113      IF( ln_zps      )   CALL zgr_zps          ! Partial step z-coordinate 
     114      IF( ln_sco      )   CALL zgr_sco          ! s-coordinate or hybrid z-s coordinate 
     115      ! 
     116      ! final adjustment of mbathy & check  
     117      ! ----------------------------------- 
     118      IF( lzoom       )   CALL zgr_bat_zoom     ! correct mbathy in case of zoom subdomain 
     119      IF( .NOT.lk_c1d )   CALL zgr_bat_ctl      ! check bathymetry (mbathy) and suppress isoated ocean points 
     120                          CALL zgr_bot_level    ! deepest ocean level for t-, u- and v-points 
     121      ! 
     122      ! 
    126123      IF( nprint == 1 .AND. lwp )   THEN 
    127124         WRITE(numout,*) ' MIN val mbathy ', MINVAL( mbathy(:,:) ), ' MAX ', MAXVAL( mbathy(:,:) ) 
    … …  
    140137            &                   ' w ',   MAXVAL( fse3w(:,:,:) ) 
    141138      ENDIF 
    142 !!!bug gm 
    143  
     139      ! 
    144140   END SUBROUTINE dom_zgr 
    145141 
    … …  
    171167      REAL(wp) ::   zacr, zdzmin, zhmax    ! par_CONFIG_Rxx.h90 
    172168      REAL(wp) ::   zrefdep                ! depth of the reference level (~10m) 
     169      REAL(wp) ::   za2, zkth2, zacr2      ! Values for optional double tanh function set from parameters  
    173170      !!---------------------------------------------------------------------- 
    174171 
    … …  
    177174       zkth = ppkth       ;   zacr = ppacr 
    178175       zdzmin = ppdzmin   ;   zhmax = pphmax 
     176       zkth2 = ppkth2     ;   zacr2 = ppacr2   ! optional (ldbletanh=T) double tanh parameters 
    179177 
    180178      ! If ppa1 and ppa0 and ppsur are et to pp_to_be_computed 
    181179      !  za0, za1, zsur are computed from ppdzmin , pphmax, ppkth, ppacr 
    182       ! 
    183180      IF(   ppa1  == pp_to_be_computed  .AND.  & 
    184181         &  ppa0  == pp_to_be_computed  .AND.  & 
    … …  
    192189      ELSE 
    193190         za1 = ppa1 ;       za0 = ppa0 ;          zsur = ppsur 
     191         za2 = ppa2                            ! optional (ldbletanh=T) double tanh parameter 
    194192      ENDIF 
    195193 
    … …  
    198196         WRITE(numout,*) '    zgr_z   : Reference vertical z-coordinates' 
    199197         WRITE(numout,*) '    ~~~~~~~' 
    200          IF(  ppkth == 0. ) THEN               
     198         IF(  ppkth == 0._wp ) THEN               
    201199              WRITE(numout,*) '            Uniform grid with ',jpk-1,' layers' 
    202200              WRITE(numout,*) '            Total depth    :', zhmax 
    203201              WRITE(numout,*) '            Layer thickness:', zhmax/(jpk-1) 
    204202         ELSE 
    205             IF( ppa1 == 0. .AND. ppa0 == 0. .AND. ppsur == 0. ) THEN 
     203            IF( ppa1 == 0._wp .AND. ppa0 == 0._wp .AND. ppsur == 0._wp ) THEN 
    206204               WRITE(numout,*) '         zsur, za0, za1 computed from ' 
    207205               WRITE(numout,*) '                 zdzmin = ', zdzmin 
    … …  
    214212            WRITE(numout,*) '                 zkth = ', zkth 
    215213            WRITE(numout,*) '                 zacr = ', zacr 
     214            IF( ldbletanh ) THEN 
     215               WRITE(numout,*) ' (Double tanh    za2  = ', za2 
     216               WRITE(numout,*) '  parameters)    zkth2= ', zkth2 
     217               WRITE(numout,*) '                 zacr2= ', zacr2 
     218            ENDIF 
    216219         ENDIF 
    217220      ENDIF 
    … …  
    220223      ! Reference z-coordinate (depth - scale factor at T- and W-points) 
    221224      ! ====================== 
    222       IF( ppkth == 0.e0 ) THEN            !  uniform vertical grid        
     225      IF( ppkth == 0._wp ) THEN            !  uniform vertical grid        
    223226         za1 = zhmax / FLOAT(jpk-1)  
    224227         DO jk = 1, jpk 
    225228            zw = FLOAT( jk ) 
    226             zt = FLOAT( jk ) + 0.5 
     229            zt = FLOAT( jk ) + 0.5_wp 
    227230            gdepw_0(jk) = ( zw - 1 ) * za1 
    228231            gdept_0(jk) = ( zt - 1 ) * za1 
    … …  
    231234         END DO 
    232235      ELSE                                ! Madec & Imbard 1996 function 
    233          DO jk = 1, jpk 
    234             zw = FLOAT( jk ) 
    235             zt = FLOAT( jk ) + 0.5 
    236             gdepw_0(jk) = ( zsur + za0 * zw + za1 * zacr * LOG ( COSH( (zw-zkth) / zacr ) )  ) 
    237             gdept_0(jk) = ( zsur + za0 * zt + za1 * zacr * LOG ( COSH( (zt-zkth) / zacr ) )  ) 
    238             e3w_0  (jk) =          za0      + za1        * TANH(       (zw-zkth) / zacr   ) 
    239             e3t_0  (jk) =          za0      + za1        * TANH(       (zt-zkth) / zacr   ) 
    240          END DO 
    241          gdepw_0(1) = 0.e0                     ! force first w-level to be exactly at zero 
     236         IF( .NOT. ldbletanh ) THEN 
     237            DO jk = 1, jpk 
     238               zw = REAL( jk , wp ) 
     239               zt = REAL( jk , wp ) + 0.5_wp 
     240               gdepw_0(jk) = ( zsur + za0 * zw + za1 * zacr * LOG ( COSH( (zw-zkth) / zacr ) )  ) 
     241               gdept_0(jk) = ( zsur + za0 * zt + za1 * zacr * LOG ( COSH( (zt-zkth) / zacr ) )  ) 
     242               e3w_0  (jk) =          za0      + za1        * TANH(       (zw-zkth) / zacr   ) 
     243               e3t_0  (jk) =          za0      + za1        * TANH(       (zt-zkth) / zacr   ) 
     244            END DO 
     245         ELSE 
     246            DO jk = 1, jpk 
     247               zw = FLOAT( jk ) 
     248               zt = FLOAT( jk ) + 0.5_wp 
     249               ! Double tanh function 
     250               gdepw_0(jk) = ( zsur + za0 * zw + za1 * zacr * LOG ( COSH( (zw-zkth ) / zacr  ) )    & 
     251                  &                            + za2 * zacr2* LOG ( COSH( (zw-zkth2) / zacr2 ) )  ) 
     252               gdept_0(jk) = ( zsur + za0 * zt + za1 * zacr * LOG ( COSH( (zt-zkth ) / zacr  ) )    & 
     253                  &                            + za2 * zacr2* LOG ( COSH( (zt-zkth2) / zacr2 ) )  ) 
     254               e3w_0  (jk) =          za0      + za1        * TANH(       (zw-zkth ) / zacr  )    & 
     255                  &                            + za2        * TANH(       (zw-zkth2) / zacr2 ) 
     256               e3t_0  (jk) =          za0      + za1        * TANH(       (zt-zkth ) / zacr  )    & 
     257                  &                            + za2        * TANH(       (zt-zkth2) / zacr2 ) 
     258            END DO 
     259         ENDIF 
     260         gdepw_0(1) = 0._wp                    ! force first w-level to be exactly at zero 
    242261      ENDIF 
    243262 
    244263!!gm BUG in s-coordinate this does not work! 
    245       ! deepest/shallowest W level Above/Bellow ~10m 
    246       zrefdep = 10. - ( 0.1*MINVAL(e3w_0) )                          ! ref. depth with tolerance (10% of minimum layer thickness) 
    247       nlb10 = MINLOC( gdepw_0, mask = gdepw_0 > zrefdep, dim = 1 )   ! shallowest W level Bellow ~10m 
    248       nla10 = nlb10 - 1                                              ! deepest    W level Above  ~10m 
     264      ! deepest/shallowest W level Above/Below ~10m 
     265      zrefdep = 10._wp - 0.1_wp * MINVAL( e3w_0 )                    ! ref. depth with tolerance (10% of minimum layer thickness) 
     266      nlb10 = MINLOC( gdepw_0, mask = gdepw_0 > zrefdep, dim = 1 )   ! shallowest W level Below ~10m 
     267      nla10 = nlb10 - 1                                              ! deepest    W level Above ~10m 
    249268!!gm end bug 
    250269 
    … …  
    256275      ENDIF 
    257276      DO jk = 1, jpk                      ! control positivity 
    258          IF( e3w_0  (jk) <= 0.e0 .OR. e3t_0  (jk) <= 0.e0 )   CALL ctl_stop( ' e3w or e3t =< 0 '    ) 
    259          IF( gdepw_0(jk) <  0.e0 .OR. gdept_0(jk) <  0.e0 )   CALL ctl_stop( ' gdepw or gdept < 0 ' ) 
     277         IF( e3w_0  (jk) <= 0._wp .OR. e3t_0  (jk) <= 0._wp )   CALL ctl_stop( 'dom:zgr_z: e3w or e3t =< 0 '    ) 
     278         IF( gdepw_0(jk) <  0._wp .OR. gdept_0(jk) <  0._wp )   CALL ctl_stop( 'dom:zgr_z: gdepw or gdept < 0 ' ) 
    260279      END DO 
    261280      ! 
    … …  
    289308      !!      ntopo= 1 :   mbathy is read in 'bathy_level.nc' NetCDF file 
    290309      !!                   bathy  is read in 'bathy_meter.nc' NetCDF file 
    291       !!      C A U T I O N : mbathy will be modified during the initializa- 
    292       !!      tion phase to become the number of non-zero w-levels of a water 
    293       !!      column, with a minimum value of 1. 
    294310      !! 
    295311      !! ** Action  : - mbathy: level bathymetry (in level index) 
    296312      !!              - bathy : meter bathymetry (in meters) 
    297313      !!---------------------------------------------------------------------- 
    298       USE iom 
    299       !! 
    300314      INTEGER  ::   ji, jj, jl, jk            ! dummy loop indices 
    301315      INTEGER  ::   inum                      ! temporary logical unit 
    302316      INTEGER  ::   ii_bump, ij_bump, ih      ! bump center position 
    303       INTEGER  ::   ii0, ii1, ij0, ij1        ! indices 
     317      INTEGER  ::   ii0, ii1, ij0, ij1, ik    ! local indices 
    304318      REAL(wp) ::   r_bump , h_bump , h_oce   ! bump characteristics  
    305       REAL(wp) ::   zi     , zj     , zh      ! temporary scalars 
     319      REAL(wp) ::   zi, zj, zh, zhmin         ! local scalars 
    306320      INTEGER , DIMENSION(jpidta,jpjdta) ::   idta   ! global domain integer data 
    307321      REAL(wp), DIMENSION(jpidta,jpjdta) ::   zdta   ! global domain scalar data 
    … …  
    328342            ii_bump = jpidta / 2                           ! i-index of the bump center 
    329343            ij_bump = jpjdta / 2                           ! j-index of the bump center 
    330             r_bump  = 50000.e0                             ! bump radius (meters)        
    331             h_bump  = 2700.e0                              ! bump height (meters) 
     344            r_bump  = 50000._wp                            ! bump radius (meters)        
     345            h_bump  =  2700._wp                            ! bump height (meters) 
    332346            h_oce   = gdepw_0(jpk)                         ! background ocean depth (meters) 
    333347            IF(lwp) WRITE(numout,*) '            bump characteristics: ' 
    … …  
    350364               idta(:,:) = jpkm1 
    351365               DO jk = 1, jpkm1 
    352                   DO jj = 1, jpjdta 
    353                      DO ji = 1, jpidta 
    354                         IF( gdept_0(jk) < zdta(ji,jj) .AND. zdta(ji,jj) <= gdept_0(jk+1) )   idta(ji,jj) = jk 
    355                      END DO 
    356                   END DO 
     366                  WHERE( gdept_0(jk) < zdta(:,:) .AND. zdta(:,:) <= gdept_0(jk+1) )   idta(:,:) = jk 
    357367               END DO 
    358368            ENDIF 
    … …  
    361371         !                                            ! Caution : idta on the global domain: use of jperio, not nperio 
    362372         IF( jperio == 1 .OR. jperio == 4 .OR. jperio == 6 ) THEN 
    363             idta( :    , 1    ) = -1                ;      zdta( :    , 1    ) = -1.e0 
    364             idta( :    ,jpjdta) =  0                ;      zdta( :    ,jpjdta) =  0.e0 
     373            idta( :    , 1    ) = -1                ;      zdta( :    , 1    ) = -1._wp 
     374            idta( :    ,jpjdta) =  0                ;      zdta( :    ,jpjdta) =  0._wp 
    365375         ELSEIF( jperio == 2 ) THEN 
    366376            idta( :    , 1    ) = idta( : ,  3  )   ;      zdta( :    , 1    ) = zdta( : ,  3  ) 
    367             idta( :    ,jpjdta) = 0                 ;      zdta( :    ,jpjdta) =  0.e0 
    368             idta( 1    , :    ) = 0                 ;      zdta( 1    , :    ) =  0.e0 
    369             idta(jpidta, :    ) = 0                 ;      zdta(jpidta, :    ) =  0.e0 
     377            idta( :    ,jpjdta) = 0                 ;      zdta( :    ,jpjdta) =  0._wp 
     378            idta( 1    , :    ) = 0                 ;      zdta( 1    , :    ) =  0._wp 
     379            idta(jpidta, :    ) = 0                 ;      zdta(jpidta, :    ) =  0._wp 
    370380         ELSE 
    371             ih = 0                                  ;      zh = 0.e0 
     381            ih = 0                                  ;      zh = 0._wp 
    372382            IF( ln_sco )   ih = jpkm1               ;      IF( ln_sco )   zh = h_oce 
    373383            idta( :    , 1    ) = ih                ;      zdta( :    , 1    ) =  zh 
    … …  
    379389         !                                            ! local domain level and meter bathymetries (mbathy,bathy) 
    380390         mbathy(:,:) = 0                                   ! set to zero extra halo points 
    381          bathy (:,:) = 0.e0                                ! (require for mpp case) 
     391         bathy (:,:) = 0._wp                               ! (require for mpp case) 
    382392         DO jj = 1, nlcj                                   ! interior values 
    383393            DO ji = 1, nlci 
    … …  
    392402         ! 
    393403         IF( ln_zco )   THEN                          ! zco : read level bathymetry  
    394             CALL iom_open( 'bathy_level.nc', inum )   
    395             CALL iom_get ( inum, jpdom_data, 'Bathy_level', bathy ) 
    396             CALL iom_close (inum) 
     404            CALL iom_open ( 'bathy_level.nc', inum )   
     405            CALL iom_get  ( inum, jpdom_data, 'Bathy_level', bathy ) 
     406            CALL iom_close( inum ) 
    397407            mbathy(:,:) = INT( bathy(:,:) ) 
    398408            !                                                ! ===================== 
    399409            IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN    ! ORCA R2 configuration 
    400410               !                                             ! ===================== 
    401                IF( n_cla == 0 ) THEN 
    402                   ! 
     411               IF( nn_cla == 0 ) THEN 
    403412                  ii0 = 140   ;   ii1 = 140                  ! Gibraltar Strait open  
    404413                  ij0 = 102   ;   ij1 = 102                  ! (Thomson, Ocean Modelling, 1995) 
    … …  
    409418                  END DO 
    410419                  IF(lwp) WRITE(numout,*) 
    411                   IF(lwp) WRITE(numout,*) '             orca_r2: Gibraltar strait open at i=',ii0,' j=',ij0 
     420                  IF(lwp) WRITE(numout,*) '      orca_r2: Gibraltar strait open at i=',ii0,' j=',ij0 
    412421                  ! 
    413422                  ii0 = 160   ;   ii1 = 160                  ! Bab el mandeb Strait open 
    … …  
    419428                  END DO 
    420429                  IF(lwp) WRITE(numout,*) 
    421                   IF(lwp) WRITE(numout,*) '             orca_r2: Bab el Mandeb strait open at i=',ii0,' j=',ij0 
    422                   ! 
     430                  IF(lwp) WRITE(numout,*) '      orca_r2: Bab el Mandeb strait open at i=',ii0,' j=',ij0 
    423431               ENDIF 
    424432               ! 
    … …  
    427435         ENDIF 
    428436         IF( ln_zps .OR. ln_sco )   THEN              ! zps or sco : read meter bathymetry 
    429             CALL iom_open( 'bathy_meter.nc', inum )  
    430             CALL iom_get ( inum, jpdom_data, 'Bathymetry', bathy ) 
    431             CALL iom_close (inum) 
     437            CALL iom_open ( 'bathy_meter.nc', inum )  
     438            CALL iom_get  ( inum, jpdom_data, 'Bathymetry', bathy ) 
     439            CALL iom_close( inum ) 
    432440            !                                                ! ===================== 
    433            IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN    ! ORCA R2 configuration 
    434               !                                             ! ===================== 
    435               IF( n_cla == 0 ) THEN 
    436                  ! 
    437                  ii0 = 140   ;   ii1 = 140                  ! Gibraltar Strait open  
    438                  ij0 = 102   ;   ij1 = 102                  ! (Thomson, Ocean Modelling, 1995) 
     441            IF( cp_cfg == "orca" .AND. jp_cfg == 1 ) THEN    ! ORCA R1 configuration 
     442               ii0 = 142   ;   ii1 = 142                     ! ===================== 
     443               ij0 =  51   ;   ij1 =  53                      
     444               DO ji = mi0(ii0), mi1(ii1)                    ! Close Halmera Strait 
     445                  DO jj = mj0(ij0), mj1(ij1) 
     446                     bathy(ji,jj) = 0._wp 
     447                  END DO 
     448               END DO 
     449               IF(lwp) WRITE(numout,*) 
     450               IF(lwp) WRITE(numout,*) '      orca_r1: Halmera strait closed at i=',ii0,' j=',ij0,'->',ij1 
     451            ENDIF 
     452            !                                                ! ===================== 
     453            IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN    ! ORCA R2 configuration 
     454               !                                             ! ===================== 
     455              IF( nn_cla == 0 ) THEN 
     456                 ii0 = 140   ;   ii1 = 140                   ! Gibraltar Strait open  
     457                 ij0 = 102   ;   ij1 = 102                   ! (Thomson, Ocean Modelling, 1995) 
    439458                 DO ji = mi0(ii0), mi1(ii1) 
    440459                    DO jj = mj0(ij0), mj1(ij1) 
    441                        bathy(ji,jj) = 284.e0 
     460                       bathy(ji,jj) = 284._wp 
    442461                    END DO 
    443462                 END DO 
    444463                 IF(lwp) WRITE(numout,*) 
    445                  IF(lwp) WRITE(numout,*) '             orca_r2: Gibraltar strait open at i=',ii0,' j=',ij0 
     464                 IF(lwp) WRITE(numout,*) '      orca_r2: Gibraltar strait open at i=',ii0,' j=',ij0 
    446465                 ! 
    447                  ii0 = 160   ;   ii1 = 160                  ! Bab el mandeb Strait open 
    448                  ij0 = 88    ;   ij1 = 88                   ! (Thomson, Ocean Modelling, 1995) 
     466                 ii0 = 160   ;   ii1 = 160                   ! Bab el mandeb Strait open 
     467                 ij0 = 88    ;   ij1 = 88                    ! (Thomson, Ocean Modelling, 1995) 
    449468                 DO ji = mi0(ii0), mi1(ii1) 
    450469                    DO jj = mj0(ij0), mj1(ij1) 
    451                        bathy(ji,jj) = 137.e0 
     470                       bathy(ji,jj) = 137._wp 
    452471                    END DO 
    453472                 END DO 
    454473                 IF(lwp) WRITE(numout,*) 
    455474                 IF(lwp) WRITE(numout,*) '             orca_r2: Bab el Mandeb strait open at i=',ii0,' j=',ij0 
    456                  ! 
    457475              ENDIF 
    458476              ! 
    … …  
    473491               DO ji = ncsi1(jl), ncsi2(jl) 
    474492                  mbathy(ji,jj) = 0                   ! suppress closed seas and lakes from bathymetry 
    475                   bathy (ji,jj) = 0.e0                 
     493                  bathy (ji,jj) = 0._wp                
    476494               END DO 
    477495            END DO 
    478496         END DO 
    479497      ENDIF 
    480 #if defined key_orca_lev10 
    481       !                                               ! ================= ! 
    482       mbathy(:,:) = 10 * mbathy(:,:)                  !   key_orca_lev10  ! 
    483       !                                               ! ================= ! 
    484       IF( ln_zps .OR. ln_sco )   CALL ctl_stop( ' CAUTION: 300 levels only with level bathymetry' ) 
    485 #endif 
    486       !                          
     498      ! 
     499      !                                               ! =========================== ! 
     500      !                                               !     set a minimum depth     ! 
     501      !                                               ! =========================== ! 
     502      IF( rn_hmin < 0._wp ) THEN    ;   ik = INT( rn_hmin ) + 1                                    ! from a nb of level 
     503      ELSE                          ;   ik = MINLOC( gdepw_0, mask = gdepw_0 > rn_hmin, dim = 1 )  ! from a depth 
     504      ENDIF 
     505      zhmin = gdepw_0(ik+1)                                                         ! minimum depth = ik+1 w-levels  
     506      WHERE( bathy(:,:) <= 0._wp )   ;   bathy(:,:) = 0._wp                         ! min=0     over the lands 
     507      ELSE WHERE                     ;   bathy(:,:) = MIN(  zhmin , bathy(:,:)  )   ! min=zhmin over the oceans 
     508      END WHERE 
     509      IF(lwp) write(numout,*) 'Minimum ocean depth: ', zhmin, ' minimum number of ocean levels : ', ik 
     510      ! 
    487511   END SUBROUTINE zgr_bat 
    488512 
    … …  
    601625      IF( lk_mpp ) THEN 
    602626         zbathy(:,:) = FLOAT( mbathy(:,:) ) 
    603          CALL lbc_lnk( zbathy, 'T', 1. ) 
     627         CALL lbc_lnk( zbathy, 'T', 1._wp ) 
    604628         mbathy(:,:) = INT( zbathy(:,:) ) 
    605629      ENDIF 
    … …  
    641665         !   ... mono- or macro-tasking: T-point, >0, 2D array, no slab 
    642666         zbathy(:,:) = FLOAT( mbathy(:,:) ) 
    643          CALL lbc_lnk( zbathy, 'T', 1. ) 
     667         CALL lbc_lnk( zbathy, 'T', 1._wp ) 
    644668         mbathy(:,:) = INT( zbathy(:,:) ) 
    645669      ENDIF 
    … …  
    672696 
    673697 
     698   SUBROUTINE zgr_bot_level 
     699      !!---------------------------------------------------------------------- 
     700      !!                    ***  ROUTINE zgr_bot_level  *** 
     701      !! 
     702      !! ** Purpose :   defines the vertical index of ocean bottom (mbk. arrays) 
     703      !! 
     704      !! ** Method  :   computes from mbathy with a minimum value of 1 over land 
     705      !! 
     706      !! ** Action  :   mbkt, mbku, mbkv :   vertical indices of the deeptest  
     707      !!                                     ocean level at t-, u- & v-points 
     708      !!                                     (min value = 1 over land) 
     709      !!---------------------------------------------------------------------- 
     710      INTEGER ::   ji, jj   ! dummy loop indices 
     711      REAL(wp), DIMENSION(jpi,jpj) ::   zmbk   ! 2D workspace  
     712      !!---------------------------------------------------------------------- 
     713      ! 
     714      IF(lwp) WRITE(numout,*) 
     715      IF(lwp) WRITE(numout,*) '    zgr_bot_level : ocean bottom k-index of T-, U-, V- and W-levels ' 
     716      IF(lwp) WRITE(numout,*) '    ~~~~~~~~~~~~~' 
     717      ! 
     718      mbkt(:,:) = MAX( mbathy(:,:) , 1 )    ! bottom k-index of T-level (=1 over land) 
     719      !                                     ! bottom k-index of W-level = mbkt+1 
     720      DO jj = 1, jpjm1                      ! bottom k-index of u- (v-) level 
     721         DO ji = 1, jpim1 
     722            mbku(ji,jj) = MIN(  mbkt(ji+1,jj  ) , mbkt(ji,jj)  ) 
     723            mbkv(ji,jj) = MIN(  mbkt(ji  ,jj+1) , mbkt(ji,jj)  ) 
     724         END DO 
     725      END DO 
     726      ! converte into REAL to use lbc_lnk ; impose a min value of 1 as a zero can be set in lbclnk  
     727      zmbk(:,:) = REAL( mbku(:,:), wp )   ;   CALL lbc_lnk(zmbk,'U',1.)   ;   mbku  (:,:) = MAX( INT( zmbk(:,:) ), 1 ) 
     728      zmbk(:,:) = REAL( mbkv(:,:), wp )   ;   CALL lbc_lnk(zmbk,'V',1.)   ;   mbkv  (:,:) = MAX( INT( zmbk(:,:) ), 1 ) 
     729      ! 
     730   END SUBROUTINE zgr_bot_level 
     731 
     732 
    674733   SUBROUTINE zgr_zco 
    675734      !!---------------------------------------------------------------------- 
    … …  
    678737      !! ** Purpose :   define the z-coordinate system 
    679738      !! 
    680       !! ** Method  :   set 3D coord. arrays to reference 1D array if lk_zco=T 
     739      !! ** Method  :   set 3D coord. arrays to reference 1D array  
    681740      !!---------------------------------------------------------------------- 
    682741      INTEGER  ::   jk 
    683742      !!---------------------------------------------------------------------- 
    684743      ! 
    685       IF( .NOT.lk_zco ) THEN 
    686          DO jk = 1, jpk 
    687             fsdept(:,:,jk) = gdept_0(jk) 
    688             fsdepw(:,:,jk) = gdepw_0(jk) 
    689             fsde3w(:,:,jk) = gdepw_0(jk) 
    690             fse3t (:,:,jk) = e3t_0(jk) 
    691             fse3u (:,:,jk) = e3t_0(jk) 
    692             fse3v (:,:,jk) = e3t_0(jk) 
    693             fse3f (:,:,jk) = e3t_0(jk) 
    694             fse3w (:,:,jk) = e3w_0(jk) 
    695             fse3uw(:,:,jk) = e3w_0(jk) 
    696             fse3vw(:,:,jk) = e3w_0(jk) 
    697          END DO 
    698       ENDIF 
     744      DO jk = 1, jpk 
     745         fsdept(:,:,jk) = gdept_0(jk) 
     746         fsdepw(:,:,jk) = gdepw_0(jk) 
     747         fsde3w(:,:,jk) = gdepw_0(jk) 
     748         fse3t (:,:,jk) = e3t_0(jk) 
     749         fse3u (:,:,jk) = e3t_0(jk) 
     750         fse3v (:,:,jk) = e3t_0(jk) 
     751         fse3f (:,:,jk) = e3t_0(jk) 
     752         fse3w (:,:,jk) = e3w_0(jk) 
     753         fse3uw(:,:,jk) = e3w_0(jk) 
     754         fse3vw(:,:,jk) = e3w_0(jk) 
     755      END DO 
    699756      ! 
    700757   END SUBROUTINE zgr_zco 
    701758 
    702 #if defined key_zco 
    703    !!---------------------------------------------------------------------- 
    704    !!   'key_zco' :                                              "pure" zco (gdep & e3 are 1D arrays) 
    705    !!---------------------------------------------------------------------- 
    706    SUBROUTINE zgr_zps      ! Empty routine 
    707    END SUBROUTINE zgr_zps 
    708    SUBROUTINE zgr_sco      ! Empty routine 
    709    END SUBROUTINE zgr_sco 
    710  
    711 #else 
    712    !!---------------------------------------------------------------------- 
    713    !!   Default option :                      zco, zps and/or sco available (gedp & e3 are 3D arrays) 
    714    !!---------------------------------------------------------------------- 
    715759 
    716760   SUBROUTINE zgr_zps 
    … …  
    764808      REAL(wp) ::   ze3tp , ze3wp    ! Last ocean level thickness at T- and W-points 
    765809      REAL(wp) ::   zdepwp, zdepth   ! Ajusted ocean depth to avoid too small e3t 
    766       REAL(wp) ::   zmax, zmin       ! Maximum and minimum depth 
     810      REAL(wp) ::   zmax             ! Maximum depth 
    767811      REAL(wp) ::   zdiff            ! temporary scalar 
     812      REAL(wp) ::   zrefdep          ! temporary scalar 
    768813      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zprt   ! 3D workspace 
    769814      !!--------------------------------------------------------------------- 
    … …  
    774819      IF(lwp) WRITE(numout,*) '              mbathy is recomputed : bathy_level file is NOT used' 
    775820 
    776       ll_print=.FALSE.                     ! Local variable for debugging 
    777 !!    ll_print=.TRUE. 
     821      ll_print = .FALSE.                   ! Local variable for debugging 
     822!!    ll_print = .TRUE. 
    778823       
    779824      IF(lwp .AND. ll_print) THEN          ! control print of the ocean depth 
    … …  
    786831      ! bathymetry in level (from bathy_meter) 
    787832      ! =================== 
    788       zmax = gdepw_0(jpk) + e3t_0(jpk)     ! maximum depth (i.e. the last ocean level thickness <= 2*e3t_0(jpkm1) ) 
    789       zmin = gdepw_0(4)                    ! minimum depth = 3 levels 
    790  
    791       mbathy(:,:) = jpkm1                  ! initialize mbathy to the maximum ocean level available 
    792  
    793       !                                    ! storage of land and island's number (zera and negative values) in mbathy 
    794       WHERE( bathy(:,:) <= 0. )   mbathy(:,:) = INT( bathy(:,:) ) 
    795  
    796       ! bounded value of bathy 
    797 !!gm  bathy(:,:) = MIN(  zmax,  MAX( bathy(:,:), zmin )  )     !!gm : simpler   as zmin is > 0 
    798       DO jj = 1, jpj 
    799          DO ji= 1, jpi 
    800             IF( bathy(ji,jj) <= 0. ) THEN   ;   bathy(ji,jj) = 0.e0 
    801             ELSE                            ;   bathy(ji,jj) = MIN(  zmax,  MAX( bathy(ji,jj), zmin )  ) 
    802             ENDIF 
    803          END DO 
    804       END DO 
     833      zmax = gdepw_0(jpk) + e3t_0(jpk)          ! maximum depth (i.e. the last ocean level thickness <= 2*e3t_0(jpkm1) ) 
     834      bathy(:,:) = MIN( zmax ,  bathy(:,:) )    ! bounded value of bathy (min already set at the end of zgr_bat) 
     835      WHERE( bathy(:,:) == 0._wp )   ;   mbathy(:,:) = 0       ! land  : set mbathy to 0 
     836      ELSE WHERE                     ;   mbathy(:,:) = jpkm1   ! ocean : initialize mbathy to the max ocean level 
     837      END WHERE 
    805838 
    806839      ! Compute mbathy for ocean points (i.e. the number of ocean levels) 
    … …  
    810843      DO jk = jpkm1, 1, -1 
    811844         zdepth = gdepw_0(jk) + MIN( e3zps_min, e3t_0(jk)*e3zps_rat ) 
    812          DO jj = 1, jpj 
    813             DO ji = 1, jpi 
    814                IF( 0. < bathy(ji,jj) .AND. bathy(ji,jj) <= zdepth )   mbathy(ji,jj) = jk-1 
    815             END DO 
    816          END DO 
     845         WHERE( 0._wp < bathy(:,:) .AND. bathy(:,:) <= zdepth )   mbathy(:,:) = jk-1 
    817846      END DO 
    818847 
    … …  
    824853         e3w  (:,:,jk) = e3w_0  (jk) 
    825854      END DO 
    826       hdept(:,:) = gdept(:,:,2 ) 
    827       hdepw(:,:) = gdepw(:,:,3 )      
    828855      !  
    829856      DO jj = 1, jpj 
    … …  
    835862                  zdepwp = bathy(ji,jj) 
    836863                  ze3tp  = bathy(ji,jj) - gdepw_0(ik) 
    837                   ze3wp = 0.5 * e3w_0(ik) * ( 1. + ( ze3tp/e3t_0(ik) ) ) 
     864                  ze3wp = 0.5_wp * e3w_0(ik) * ( 1._wp + ( ze3tp/e3t_0(ik) ) ) 
    838865                  e3t(ji,jj,ik  ) = ze3tp 
    839866                  e3t(ji,jj,ik+1) = ze3tp 
    … …  
    855882                  e3t  (ji,jj,ik) = e3t_0  (ik) * ( gdepw  (ji,jj,ik+1) - gdepw_0(ik) )   &  
    856883                     &                          / ( gdepw_0(      ik+1) - gdepw_0(ik) )  
    857                   e3w  (ji,jj,ik) = 0.5 * ( gdepw(ji,jj,ik+1) + gdepw_0(ik+1) - 2.*gdepw_0(ik) )   & 
    858                      &                  * ( e3w_0(ik) / ( gdepw_0(ik+1) - gdepw_0(ik) ) ) 
     884                  e3w  (ji,jj,ik) = 0.5_wp * ( gdepw(ji,jj,ik+1) + gdepw_0(ik+1) - 2._wp * gdepw_0(ik) )   & 
     885                     &                     * ( e3w_0(ik) / ( gdepw_0(ik+1) - gdepw_0(ik) ) ) 
    859886                  !       ... on ik+1 
    860887                  e3w  (ji,jj,ik+1) = e3t  (ji,jj,ik) 
    … …  
    871898            ik = mbathy(ji,jj) 
    872899            IF( ik > 0 ) THEN               ! ocean point only 
    873                hdept(ji,jj) = gdept(ji,jj,ik  ) 
    874                hdepw(ji,jj) = gdepw(ji,jj,ik+1) 
    875900               e3tp (ji,jj) = e3t(ji,jj,ik  ) 
    876901               e3wp (ji,jj) = e3w(ji,jj,ik  ) 
    877902               ! test 
    878903               zdiff= gdepw(ji,jj,ik+1) - gdept(ji,jj,ik  ) 
    879                IF( zdiff <= 0. .AND. lwp ) THEN  
     904               IF( zdiff <= 0._wp .AND. lwp ) THEN  
    880905                  it = it + 1 
    881906                  WRITE(numout,*) ' it      = ', it, ' ik      = ', ik, ' (i,j) = ', ji, jj 
    … …  
    890915      ! Scale factors and depth at U-, V-, UW and VW-points 
    891916      DO jk = 1, jpk                        ! initialisation to z-scale factors 
    892          e3u (:,:,jk)  = e3t_0(jk) 
    893          e3v (:,:,jk)  = e3t_0(jk) 
    894          e3uw(:,:,jk)  = e3w_0(jk) 
    895          e3vw(:,:,jk)  = e3w_0(jk) 
     917         e3u (:,:,jk) = e3t_0(jk) 
     918         e3v (:,:,jk) = e3t_0(jk) 
     919         e3uw(:,:,jk) = e3w_0(jk) 
     920         e3vw(:,:,jk) = e3w_0(jk) 
    896921      END DO 
    897922      DO jk = 1,jpk                         ! Computed as the minimum of neighbooring scale factors 
    898923         DO jj = 1, jpjm1 
    899924            DO ji = 1, fs_jpim1   ! vector opt. 
    900                e3u (ji,jj,jk) = MIN( e3t(ji,jj,jk), e3t(ji+1,jj,jk)) 
    901                e3v (ji,jj,jk) = MIN( e3t(ji,jj,jk), e3t(ji,jj+1,jk)) 
     925               e3u (ji,jj,jk) = MIN( e3t(ji,jj,jk), e3t(ji+1,jj,jk) ) 
     926               e3v (ji,jj,jk) = MIN( e3t(ji,jj,jk), e3t(ji,jj+1,jk) ) 
    902927               e3uw(ji,jj,jk) = MIN( e3w(ji,jj,jk), e3w(ji+1,jj,jk) ) 
    903928               e3vw(ji,jj,jk) = MIN( e3w(ji,jj,jk), e3w(ji,jj+1,jk) ) 
    … …  
    905930         END DO 
    906931      END DO 
    907       !                                     ! Boundary conditions 
    908       CALL lbc_lnk( e3u , 'U', 1. )   ;   CALL lbc_lnk( e3uw, 'U', 1. ) 
    909       CALL lbc_lnk( e3v , 'V', 1. )   ;   CALL lbc_lnk( e3vw, 'V', 1. ) 
     932      CALL lbc_lnk( e3u , 'U', 1._wp )   ;   CALL lbc_lnk( e3uw, 'U', 1._wp )   ! lateral boundary conditions 
     933      CALL lbc_lnk( e3v , 'V', 1._wp )   ;   CALL lbc_lnk( e3vw, 'V', 1._wp ) 
    910934      ! 
    911935      DO jk = 1, jpk                        ! set to z-scale factor if zero (i.e. along closed boundaries) 
    912          WHERE( e3u (:,:,jk) == 0.e0 )   e3u (:,:,jk) = e3t_0(jk) 
    913          WHERE( e3v (:,:,jk) == 0.e0 )   e3v (:,:,jk) = e3t_0(jk) 
    914          WHERE( e3uw(:,:,jk) == 0.e0 )   e3uw(:,:,jk) = e3w_0(jk) 
    915          WHERE( e3vw(:,:,jk) == 0.e0 )   e3vw(:,:,jk) = e3w_0(jk) 
     936         WHERE( e3u (:,:,jk) == 0._wp )   e3u (:,:,jk) = e3t_0(jk) 
     937         WHERE( e3v (:,:,jk) == 0._wp )   e3v (:,:,jk) = e3t_0(jk) 
     938         WHERE( e3uw(:,:,jk) == 0._wp )   e3uw(:,:,jk) = e3w_0(jk) 
     939         WHERE( e3vw(:,:,jk) == 0._wp )   e3vw(:,:,jk) = e3w_0(jk) 
    916940      END DO 
    917941       
    918942      ! Scale factor at F-point 
    919943      DO jk = 1, jpk                        ! initialisation to z-scale factors 
    920          e3f (:,:,jk) = e3t_0(jk) 
     944         e3f(:,:,jk) = e3t_0(jk) 
    921945      END DO 
    922946      DO jk = 1, jpk                        ! Computed as the minimum of neighbooring V-scale factors 
    … …  
    927951         END DO 
    928952      END DO 
    929       CALL lbc_lnk( e3f, 'F', 1. )          ! Boundary conditions 
     953      CALL lbc_lnk( e3f, 'F', 1._wp )       ! Lateral boundary conditions 
    930954      ! 
    931955      DO jk = 1, jpk                        ! set to z-scale factor if zero (i.e. along closed boundaries) 
    932          WHERE( e3f(:,:,jk) == 0.e0 )   e3f(:,:,jk) = e3t_0(jk) 
     956         WHERE( e3f(:,:,jk) == 0._wp )   e3f(:,:,jk) = e3t_0(jk) 
    933957      END DO 
    934958!!gm  bug ? :  must be a do loop with mj0,mj1 
    … …  
    941965 
    942966      ! Control of the sign 
    943       IF( MINVAL( e3t  (:,:,:) ) <= 0.e0 )   CALL ctl_stop( '    zgr_zps :   e r r o r   e3t   <= 0' ) 
    944       IF( MINVAL( e3w  (:,:,:) ) <= 0.e0 )   CALL ctl_stop( '    zgr_zps :   e r r o r   e3w   <= 0' ) 
    945       IF( MINVAL( gdept(:,:,:) ) <  0.e0 )   CALL ctl_stop( '    zgr_zps :   e r r o r   gdepw <  0' ) 
    946       IF( MINVAL( gdepw(:,:,:) ) <  0.e0 )   CALL ctl_stop( '    zgr_zps :   e r r o r   gdepw <  0' ) 
     967      IF( MINVAL( e3t  (:,:,:) ) <= 0._wp )   CALL ctl_stop( '    zgr_zps :   e r r o r   e3t   <= 0' ) 
     968      IF( MINVAL( e3w  (:,:,:) ) <= 0._wp )   CALL ctl_stop( '    zgr_zps :   e r r o r   e3w   <= 0' ) 
     969      IF( MINVAL( gdept(:,:,:) ) <  0._wp )   CALL ctl_stop( '    zgr_zps :   e r r o r   gdepw <  0' ) 
     970      IF( MINVAL( gdepw(:,:,:) ) <  0._wp )   CALL ctl_stop( '    zgr_zps :   e r r o r   gdepw <  0' ) 
    947971      
    948972      ! Compute gdep3w (vertical sum of e3w) 
    949       gdep3w(:,:,1) = 0.5 * e3w(:,:,1) 
     973      gdep3w(:,:,1) = 0.5_wp * e3w(:,:,1) 
    950974      DO jk = 2, jpk 
    951975         gdep3w(:,:,jk) = gdep3w(:,:,jk-1) + e3w(:,:,jk)  
    952976      END DO 
    953            
     977         
    954978      !                                               ! ================= ! 
    955979      IF(lwp .AND. ll_print) THEN                     !   Control print   ! 
    … …  
    9791003         WRITE(numout,*) 'domzgr gdep3w(mbathy)'   ;   CALL prihre(zprt(:,:,1),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout) 
    9801004      ENDIF   
    981       !                                      
     1005      ! 
    9821006   END SUBROUTINE zgr_zps 
    9831007 
    … …  
    9931017      !!                T-points at integer values (between 1 and jpk) 
    9941018      !!                W-points at integer values - 1/2 (between 0.5 and jpk-0.5) 
    995       !! 
    996       !! Reference  :   ??? 
    997       !!---------------------------------------------------------------------- 
    998       REAL(wp), INTENT(in   ) ::   pk   ! continuous "k" coordinate 
    999       REAL(wp)                ::   pf   ! sigma value 
    1000       !!---------------------------------------------------------------------- 
    1001       ! 
    1002       pf =   (   TANH( rn_theta * ( -(pk-0.5) / REAL(jpkm1) + rn_thetb )  )      & 
     1019      !!---------------------------------------------------------------------- 
     1020      REAL(wp), INTENT(in) ::   pk   ! continuous "k" coordinate 
     1021      REAL(wp)             ::   pf   ! sigma value 
     1022      !!---------------------------------------------------------------------- 
     1023      ! 
     1024      pf =   (   TANH( rn_theta * ( -(pk-0.5_wp) / REAL(jpkm1) + rn_thetb )  )   & 
    10031025         &     - TANH( rn_thetb * rn_theta                                )  )   & 
    1004          & * (   COSH( rn_theta                           )                   & 
    1005          &     + COSH( rn_theta * ( 2.e0 * rn_thetb - 1.e0 ) )  )                & 
    1006          & / ( 2.e0 * SINH( rn_theta ) ) 
     1026         & * (   COSH( rn_theta                           )                      & 
     1027         &     + COSH( rn_theta * ( 2._wp * rn_thetb - 1._wp ) )  )              & 
     1028         & / ( 2._wp * SINH( rn_theta ) ) 
    10071029      ! 
    10081030   END FUNCTION fssig 
    … …  
    10191041      !!                T-points at integer values (between 1 and jpk) 
    10201042      !!                W-points at integer values - 1/2 (between 0.5 and jpk-0.5) 
    1021       !! 
    1022       !! Reference  :   ??? 
    1023       !!---------------------------------------------------------------------- 
    1024       REAL(wp), INTENT(in   ) ::   pk1   ! continuous "k" coordinate 
    1025       REAL(wp), INTENT(in   ) ::   pbb   ! Stretching coefficient 
    1026       REAL(wp)                ::   pf1   ! sigma value 
     1043      !!---------------------------------------------------------------------- 
     1044      REAL(wp), INTENT(in) ::   pk1   ! continuous "k" coordinate 
     1045      REAL(wp), INTENT(in) ::   pbb   ! Stretching coefficient 
     1046      REAL(wp)             ::   pf1   ! sigma value 
    10271047      !!---------------------------------------------------------------------- 
    10281048      ! 
    10291049      IF ( rn_theta == 0 ) then      ! uniform sigma 
    1030          pf1 = -(pk1-0.5) / REAL( jpkm1 ) 
     1050         pf1 = - ( pk1 - 0.5_wp ) / REAL( jpkm1 ) 
    10311051      ELSE                        ! stretched sigma 
    1032          pf1 =   (1.0-pbb) * (sinh( rn_theta*(-(pk1-0.5)/REAL(jpkm1)) ) ) / sinh(rn_theta) + & 
    1033             &    pbb * ( (tanh( rn_theta*( (-(pk1-0.5)/REAL(jpkm1)) + 0.5) ) - tanh(0.5*rn_theta) ) / & 
    1034             &    (2*tanh(0.5*rn_theta) ) ) 
     1052         pf1 =   ( 1._wp - pbb ) * ( SINH( rn_theta*(-(pk1-0.5_wp)/REAL(jpkm1)) ) ) / SINH( rn_theta )              & 
     1053            &  + pbb * (  (TANH( rn_theta*( (-(pk1-0.5_wp)/REAL(jpkm1)) + 0.5_wp) ) - TANH( 0.5_wp * rn_theta )  ) & 
     1054            &        / ( 2._wp * TANH( 0.5_wp * rn_theta ) ) ) 
    10351055      ENDIF 
    10361056      ! 
    … …  
    11091129      ENDIF 
    11101130 
    1111       gsigw3  = 0.0d0   ;   gsigt3  = 0.0d0   ; gsi3w3 = 0.0d0 
    1112       esigt3  = 0.0d0   ;   esigw3  = 0.0d0  
    1113       esigtu3 = 0.0d0   ;   esigtv3 = 0.0d0   ; esigtf3 = 0.0d0 
    1114       esigwu3 = 0.0d0   ;   esigwv3 = 0.0d0 
     1131      gsigw3  = 0._wp   ;   gsigt3  = 0._wp   ;   gsi3w3  = 0._wp 
     1132      esigt3  = 0._wp   ;   esigw3  = 0._wp  
     1133      esigtu3 = 0._wp   ;   esigtv3 = 0._wp   ;   esigtf3 = 0._wp 
     1134      esigwu3 = 0._wp   ;   esigwv3 = 0._wp 
    11151135 
    11161136      hift(:,:) = rn_sbot_min                     ! set the minimum depth for the s-coordinate 
    … …  
    11241144      DO jj = 1, jpj 
    11251145         DO ji = 1, jpi 
    1126            IF (bathy(ji,jj) .gt. 0.0) THEN 
     1146           IF( bathy(ji,jj) > 0._wp ) THEN 
    11271147              bathy(ji,jj) = MAX( rn_sbot_min, bathy(ji,jj) ) 
    11281148           ENDIF 
    … …  
    11401160      !  
    11411161      ! Smooth the bathymetry (if required) 
    1142       scosrf(:,:) = 0.e0              ! ocean surface depth (here zero: no under ice-shelf sea) 
     1162      scosrf(:,:) = 0._wp             ! ocean surface depth (here zero: no under ice-shelf sea) 
    11431163      scobot(:,:) = bathy(:,:)        ! ocean bottom  depth 
    11441164      ! 
    11451165      jl = 0 
    1146       zrmax = 1.e0 
     1166      zrmax = 1._wp 
    11471167      !                                                     ! ================ ! 
    1148       DO WHILE ( jl <= 10000 .AND. zrmax > rn_rmax )          !  Iterative loop  ! 
     1168      DO WHILE( jl <= 10000 .AND. zrmax > rn_rmax )         !  Iterative loop  ! 
    11491169         !                                                  ! ================ ! 
    11501170         jl = jl + 1 
    1151          zrmax = 0.e0 
    1152          zmsk(:,:) = 0.e0 
     1171         zrmax = 0._wp 
     1172         zmsk(:,:) = 0._wp 
    11531173         DO jj = 1, nlcj 
    11541174            DO ji = 1, nlci 
    … …  
    11581178               zrj(ji,jj) = ABS( zenv(ji  ,ijp1) - zenv(ji,jj) ) / ( zenv(ji  ,ijp1) + zenv(ji,jj) ) 
    11591179               zrmax = MAX( zrmax, zri(ji,jj), zrj(ji,jj) ) 
    1160                IF( zri(ji,jj) > rn_rmax )   zmsk(ji  ,jj  ) = 1.0 
    1161                IF( zri(ji,jj) > rn_rmax )   zmsk(iip1,jj  ) = 1.0 
    1162                IF( zrj(ji,jj) > rn_rmax )   zmsk(ji  ,jj  ) = 1.0 
    1163                IF( zrj(ji,jj) > rn_rmax )   zmsk(ji  ,ijp1) = 1.0 
     1180               IF( zri(ji,jj) > rn_rmax )   zmsk(ji  ,jj  ) = 1._wp 
     1181               IF( zri(ji,jj) > rn_rmax )   zmsk(iip1,jj  ) = 1._wp 
     1182               IF( zrj(ji,jj) > rn_rmax )   zmsk(ji  ,jj  ) = 1._wp 
     1183               IF( zrj(ji,jj) > rn_rmax )   zmsk(ji  ,ijp1) = 1._wp 
    11641184            END DO 
    11651185         END DO 
    11661186         IF( lk_mpp )   CALL mpp_max( zrmax )   ! max over the global domain 
    11671187         ! lateral boundary condition on zmsk: keep 1 along closed boundary (use of MAX) 
    1168          ztmp(:,:) = zmsk(:,:)   ;   CALL lbc_lnk( zmsk, 'T', 1. ) 
     1188         ztmp(:,:) = zmsk(:,:)   ;   CALL lbc_lnk( zmsk, 'T', 1._wp ) 
    11691189         DO jj = 1, nlcj 
    11701190            DO ji = 1, nlci 
    … …  
    11811201               iim1 = MAX( ji-1,  1  )      ! first line (ji=nlci) 
    11821202               ijm1 = MAX( jj-1,  1  )      ! first raw  (jj=nlcj) 
    1183                IF( zmsk(ji,jj) == 1.0 ) THEN 
     1203               IF( zmsk(ji,jj) == 1._wp ) THEN 
    11841204                  ztmp(ji,jj) =   (                                                                                   & 
    11851205             &      zenv(iim1,ijp1)*zmsk(iim1,ijp1) + zenv(ji,ijp1)*zmsk(ji,ijp1) + zenv(iip1,ijp1)*zmsk(iip1,ijp1)   & 
    1186              &    + zenv(iim1,jj  )*zmsk(iim1,jj  ) + zenv(ji,jj  )*    2.e0      + zenv(iip1,jj  )*zmsk(iip1,jj  )   & 
     1206             &    + zenv(iim1,jj  )*zmsk(iim1,jj  ) + zenv(ji,jj  )*    2._wp     + zenv(iip1,jj  )*zmsk(iip1,jj  )   & 
    11871207             &    + zenv(iim1,ijm1)*zmsk(iim1,ijm1) + zenv(ji,ijm1)*zmsk(ji,ijm1) + zenv(iip1,ijm1)*zmsk(iip1,ijm1)   & 
    11881208             &                    ) / (                                                                               & 
    11891209             &                      zmsk(iim1,ijp1) +               zmsk(ji,ijp1) +                 zmsk(iip1,ijp1)   & 
    1190              &    +                 zmsk(iim1,jj  ) +                   2.e0      +                 zmsk(iip1,jj  )   & 
     1210             &    +                 zmsk(iim1,jj  ) +                   2._wp     +                 zmsk(iip1,jj  )   & 
    11911211             &    +                 zmsk(iim1,ijm1) +               zmsk(ji,ijm1) +                 zmsk(iip1,ijm1)   & 
    11921212             &                        ) 
    … …  
    11971217         DO jj = 1, nlcj 
    11981218            DO ji = 1, nlci 
    1199                IF( zmsk(ji,jj) == 1.0 )   zenv(ji,jj) = MAX( ztmp(ji,jj), bathy(ji,jj) ) 
     1219               IF( zmsk(ji,jj) == 1._wp )   zenv(ji,jj) = MAX( ztmp(ji,jj), bathy(ji,jj) ) 
    12001220            END DO 
    12011221         END DO 
    12021222         ! 
    12031223         ! Apply lateral boundary condition   CAUTION: kept the value when the lbc field is zero 
    1204          ztmp(:,:) = zenv(:,:)   ;   CALL lbc_lnk( zenv, 'T', 1. ) 
     1224         ztmp(:,:) = zenv(:,:)   ;   CALL lbc_lnk( zenv, 'T', 1._wp ) 
    12051225         DO jj = 1, nlcj 
    12061226            DO ji = 1, nlci 
    1207                IF( zenv(ji,jj) == 0.e0 )   zenv(ji,jj) = ztmp(ji,jj) 
     1227               IF( zenv(ji,jj) == 0._wp )   zenv(ji,jj) = ztmp(ji,jj) 
    12081228            END DO 
    12091229         END DO 
    … …  
    12141234      !                                        ! envelop bathymetry saved in hbatt 
    12151235      hbatt(:,:) = zenv(:,:)  
    1216       IF( MINVAL( gphit(:,:) ) * MAXVAL( gphit(:,:) ) <= 0.e0 ) THEN 
     1236      IF( MINVAL( gphit(:,:) ) * MAXVAL( gphit(:,:) ) <= 0._wp ) THEN 
    12171237         CALL ctl_warn( ' s-coordinates are tapered in vicinity of the Equator' ) 
    12181238         DO jj = 1, jpj 
    12191239            DO ji = 1, jpi 
    1220                ztaper = EXP( -(gphit(ji,jj)/8)**2 ) 
    1221                hbatt(ji,jj) = rn_sbot_max * ztaper + hbatt(ji,jj) * (1.-ztaper) 
     1240               ztaper = EXP( -(gphit(ji,jj)/8._wp)**2 ) 
     1241               hbatt(ji,jj) = rn_sbot_max * ztaper + hbatt(ji,jj) * ( 1._wp - ztaper ) 
    12221242            END DO 
    12231243         END DO 
    … …  
    12281248         WRITE(numout,*) ' domzgr: hbatt field; ocean depth in meters' 
    12291249         WRITE(numout,*) 
    1230          CALL prihre( hbatt(1,1), jpi, jpj, 1, jpi, 1, 1, jpj, 1, 0., numout ) 
     1250         CALL prihre( hbatt(1,1), jpi, jpj, 1, jpi, 1, 1, jpj, 1, 0._wp, numout ) 
    12311251         IF( nprint == 1 )   THEN         
    12321252            WRITE(numout,*) ' bathy  MAX ', MAXVAL( bathy(:,:) ), ' MIN ', MINVAL( bathy(:,:) ) 
    … …  
    12471267      DO jj = 1, jpjm1 
    12481268        DO ji = 1, jpim1   ! NO vector opt. 
    1249            hbatu(ji,jj) = 0.5 * ( hbatt(ji  ,jj) + hbatt(ji+1,jj  ) ) 
    1250            hbatv(ji,jj) = 0.5 * ( hbatt(ji  ,jj) + hbatt(ji  ,jj+1) ) 
    1251            hbatf(ji,jj) = 0.25* ( hbatt(ji  ,jj) + hbatt(ji  ,jj+1)   & 
    1252               &                 + hbatt(ji+1,jj) + hbatt(ji+1,jj+1) ) 
     1269           hbatu(ji,jj) = 0.50_wp * ( hbatt(ji  ,jj) + hbatt(ji+1,jj  ) ) 
     1270           hbatv(ji,jj) = 0.50_wp * ( hbatt(ji  ,jj) + hbatt(ji  ,jj+1) ) 
     1271           hbatf(ji,jj) = 0.25_wp * ( hbatt(ji  ,jj) + hbatt(ji  ,jj+1)   & 
     1272              &                     + hbatt(ji+1,jj) + hbatt(ji+1,jj+1) ) 
    12531273        END DO 
    12541274      END DO 
    … …  
    12561276      ! Apply lateral boundary condition 
    12571277!!gm  ! CAUTION: retain non zero value in the initial file this should be OK for orca cfg, not for EEL 
    1258       zhbat(:,:) = hbatu(:,:)   ;   CALL lbc_lnk( hbatu, 'U', 1. ) 
     1278      zhbat(:,:) = hbatu(:,:)   ;   CALL lbc_lnk( hbatu, 'U', 1._wp ) 
    12591279      DO jj = 1, jpj 
    12601280         DO ji = 1, jpi 
    1261             IF( hbatu(ji,jj) == 0.e0 ) THEN 
    1262                IF( zhbat(ji,jj) == 0.e0 )   hbatu(ji,jj) = rn_sbot_min 
    1263                IF( zhbat(ji,jj) /= 0.e0 )   hbatu(ji,jj) = zhbat(ji,jj) 
     1281            IF( hbatu(ji,jj) == 0._wp ) THEN 
     1282               IF( zhbat(ji,jj) == 0._wp )   hbatu(ji,jj) = rn_sbot_min 
     1283               IF( zhbat(ji,jj) /= 0._wp )   hbatu(ji,jj) = zhbat(ji,jj) 
    12641284            ENDIF 
    12651285         END DO 
    12661286      END DO 
    1267       zhbat(:,:) = hbatv(:,:)   ;   CALL lbc_lnk( hbatv, 'V', 1. ) 
     1287      zhbat(:,:) = hbatv(:,:)   ;   CALL lbc_lnk( hbatv, 'V', 1._wp ) 
    12681288      DO jj = 1, jpj 
    12691289         DO ji = 1, jpi 
    1270             IF( hbatv(ji,jj) == 0.e0 ) THEN 
    1271                IF( zhbat(ji,jj) == 0.e0 )   hbatv(ji,jj) = rn_sbot_min 
    1272                IF( zhbat(ji,jj) /= 0.e0 )   hbatv(ji,jj) = zhbat(ji,jj) 
     1290            IF( hbatv(ji,jj) == 0._wp ) THEN 
     1291               IF( zhbat(ji,jj) == 0._wp )   hbatv(ji,jj) = rn_sbot_min 
     1292               IF( zhbat(ji,jj) /= 0._wp )   hbatv(ji,jj) = zhbat(ji,jj) 
    12731293            ENDIF 
    12741294         END DO 
    12751295      END DO 
    1276       zhbat(:,:) = hbatf(:,:)   ;   CALL lbc_lnk( hbatf, 'F', 1. ) 
     1296      zhbat(:,:) = hbatf(:,:)   ;   CALL lbc_lnk( hbatf, 'F', 1._wp ) 
    12771297      DO jj = 1, jpj 
    12781298         DO ji = 1, jpi 
    1279             IF( hbatf(ji,jj) == 0.e0 ) THEN 
    1280                IF( zhbat(ji,jj) == 0.e0 )   hbatf(ji,jj) = rn_sbot_min 
    1281                IF( zhbat(ji,jj) /= 0.e0 )   hbatf(ji,jj) = zhbat(ji,jj) 
     1299            IF( hbatf(ji,jj) == 0._wp ) THEN 
     1300               IF( zhbat(ji,jj) == 0._wp )   hbatf(ji,jj) = rn_sbot_min 
     1301               IF( zhbat(ji,jj) /= 0._wp )   hbatf(ji,jj) = zhbat(ji,jj) 
    12821302            ENDIF 
    12831303         END DO 
    … …  
    13131333            DO jj = 1, jpj 
    13141334 
    1315              IF (hbatt(ji,jj).GT.rn_hc) THEN !deep water, stretched sigma 
     1335               IF( hbatt(ji,jj) > rn_hc ) THEN    !deep water, stretched sigma 
     1336                  DO jk = 1, jpk 
     1337                     gsigw3(ji,jj,jk) = -fssig1( REAL(jk,wp)-0.5_wp, rn_bb ) 
     1338                     gsigt3(ji,jj,jk) = -fssig1( REAL(jk,wp)       , rn_bb ) 
     1339                  END DO 
     1340               ELSE ! shallow water, uniform sigma 
     1341                  DO jk = 1, jpk 
     1342                     gsigw3(ji,jj,jk) =   REAL(jk-1,wp)            / REAL(jpk-1,wp) 
     1343                     gsigt3(ji,jj,jk) = ( REAL(jk-1,wp) + 0.5_wp ) / REAL(jpk-1,wp) 
     1344                  END DO 
     1345               ENDIF 
     1346               IF( nprint == 1 .AND. lwp )   WRITE(numout,*) 'gsigw3 1 jpk    ', gsigw3(ji,jj,1), gsigw3(ji,jj,jpk) 
     1347               ! 
     1348               DO jk = 1, jpkm1 
     1349                  esigt3(ji,jj,jk  ) = gsigw3(ji,jj,jk+1) - gsigw3(ji,jj,jk) 
     1350                  esigw3(ji,jj,jk+1) = gsigt3(ji,jj,jk+1) - gsigt3(ji,jj,jk) 
     1351               END DO 
     1352               esigw3(ji,jj,1  ) = 2._wp * ( gsigt3(ji,jj,1  ) - gsigw3(ji,jj,1  ) ) 
     1353               esigt3(ji,jj,jpk) = 2._wp * ( gsigt3(ji,jj,jpk) - gsigw3(ji,jj,jpk) ) 
     1354               ! 
     1355               ! Coefficients for vertical depth as the sum of e3w scale factors 
     1356               gsi3w3(ji,jj,1) = 0.5_wp * esigw3(ji,jj,1) 
     1357               DO jk = 2, jpk 
     1358                  gsi3w3(ji,jj,jk) = gsi3w3(ji,jj,jk-1) + esigw3(ji,jj,jk) 
     1359               END DO 
     1360               ! 
    13161361               DO jk = 1, jpk 
    1317                   gsigw3(ji,jj,jk) = -fssig1( REAL(jk,wp)-0.5_wp, rn_bb ) 
    1318                   gsigt3(ji,jj,jk) = -fssig1( REAL(jk,wp)       , rn_bb ) 
     1362                  zcoeft = ( REAL(jk,wp) - 0.5_wp ) / REAL(jpkm1,wp) 
     1363                  zcoefw = ( REAL(jk,wp) - 1.0_wp ) / REAL(jpkm1,wp) 
     1364                  gdept (ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*gsigt3(ji,jj,jk)+rn_hc*zcoeft ) 
     1365                  gdepw (ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*gsigw3(ji,jj,jk)+rn_hc*zcoefw ) 
     1366                  gdep3w(ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*gsi3w3(ji,jj,jk)+rn_hc*zcoeft ) 
    13191367               END DO 
    1320              ELSE ! shallow water, uniform sigma 
     1368               ! 
     1369            END DO   ! for all jj's 
     1370         END DO    ! for all ji's 
     1371 
     1372         DO ji = 1, jpi 
     1373            DO jj = 1, jpj 
    13211374               DO jk = 1, jpk 
    1322                  gsigw3(ji,jj,jk) =  REAL(jk-1,wp)     /REAL(jpk-1,wp) 
    1323                  gsigt3(ji,jj,jk) = (REAL(jk-1,wp)+0.5)/REAL(jpk-1,wp) 
     1375                  esigtu3(ji,jj,jk) = ( hbatt(ji,jj)*esigt3(ji,jj,jk)+hbatt(ji+1,jj)*esigt3(ji+1,jj,jk) )   & 
     1376                     &              / ( hbatt(ji,jj)+hbatt(ji+1,jj) ) 
     1377                  esigtv3(ji,jj,jk) = ( hbatt(ji,jj)*esigt3(ji,jj,jk)+hbatt(ji,jj+1)*esigt3(ji,jj+1,jk) )   & 
     1378                     &              / ( hbatt(ji,jj)+hbatt(ji,jj+1) ) 
     1379                  esigtf3(ji,jj,jk) = ( hbatt(ji,jj)*esigt3(ji,jj,jk)+hbatt(ji+1,jj)*esigt3(ji+1,jj,jk)     & 
     1380                     &                + hbatt(ji,jj+1)*esigt3(ji,jj+1,jk)+hbatt(ji+1,jj+1)*esigt3(ji+1,jj+1,jk) )   & 
     1381                     &              / ( hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1) ) 
     1382                  esigwu3(ji,jj,jk) = ( hbatt(ji,jj)*esigw3(ji,jj,jk)+hbatt(ji+1,jj)*esigw3(ji+1,jj,jk) )   & 
     1383                     &              / ( hbatt(ji,jj)+hbatt(ji+1,jj) ) 
     1384                  esigwv3(ji,jj,jk) = ( hbatt(ji,jj)*esigw3(ji,jj,jk)+hbatt(ji,jj+1)*esigw3(ji,jj+1,jk) )   & 
     1385                     &              / ( hbatt(ji,jj)+hbatt(ji,jj+1) ) 
     1386                  ! 
     1387                  e3t(ji,jj,jk) = ( (hbatt(ji,jj)-rn_hc)*esigt3 (ji,jj,jk) + rn_hc/FLOAT(jpkm1) ) 
     1388                  e3u(ji,jj,jk) = ( (hbatu(ji,jj)-rn_hc)*esigtu3(ji,jj,jk) + rn_hc/FLOAT(jpkm1) ) 
     1389                  e3v(ji,jj,jk) = ( (hbatv(ji,jj)-rn_hc)*esigtv3(ji,jj,jk) + rn_hc/FLOAT(jpkm1) ) 
     1390                  e3f(ji,jj,jk) = ( (hbatf(ji,jj)-rn_hc)*esigtf3(ji,jj,jk) + rn_hc/FLOAT(jpkm1) ) 
     1391                  ! 
     1392                  e3w (ji,jj,jk) = ( (hbatt(ji,jj)-rn_hc)*esigw3 (ji,jj,jk) + rn_hc/FLOAT(jpkm1) ) 
     1393                  e3uw(ji,jj,jk) = ( (hbatu(ji,jj)-rn_hc)*esigwu3(ji,jj,jk) + rn_hc/FLOAT(jpkm1) ) 
     1394                  e3vw(ji,jj,jk) = ( (hbatv(ji,jj)-rn_hc)*esigwv3(ji,jj,jk) + rn_hc/FLOAT(jpkm1) ) 
    13241395               END DO 
    1325              ENDIF 
    1326              IF( nprint == 1 .AND. lwp ) WRITE(numout,*) 'gsigw3 1 jpk    ', gsigw3(ji,jj,1), gsigw3(ji,jj,jpk) 
    1327  
    1328  
    1329              DO jk = 1, jpkm1 
    1330                 esigt3(ji,jj,jk) = gsigw3(ji,jj,jk+1) - gsigw3(ji,jj,jk) 
    1331                 esigw3(ji,jj,jk+1) = gsigt3(ji,jj,jk+1) - gsigt3(ji,jj,jk) 
    1332              END DO 
    1333              esigw3(ji,jj,1  ) = 2.0_wp * (gsigt3(ji,jj,1  ) - gsigw3(ji,jj,1  )) 
    1334              esigt3(ji,jj,jpk) = 2.0_wp * (gsigt3(ji,jj,jpk) - gsigw3(ji,jj,jpk)) 
    1335  
    1336              ! Coefficients for vertical depth as the sum of e3w scale factors 
    1337              gsi3w3(ji,jj,1) = 0.5 * esigw3(ji,jj,1) 
    1338              DO jk = 2, jpk 
    1339                 gsi3w3(ji,jj,jk) = gsi3w3(ji,jj,jk-1) + esigw3(ji,jj,jk) 
    1340              END DO 
    1341  
    1342              DO jk = 1, jpk 
    1343                 zcoeft = ( REAL(jk,wp) - 0.5 ) / REAL(jpkm1,wp) 
    1344                 zcoefw = ( REAL(jk,wp) - 1.0 ) / REAL(jpkm1,wp) 
    1345                 gdept (ji,jj,jk) = (scosrf(ji,jj)+(hbatt(ji,jj)-rn_hc)*gsigt3(ji,jj,jk)+rn_hc*zcoeft) 
    1346                 gdepw (ji,jj,jk) = (scosrf(ji,jj)+(hbatt(ji,jj)-rn_hc)*gsigw3(ji,jj,jk)+rn_hc*zcoefw) 
    1347                 gdep3w(ji,jj,jk) = (scosrf(ji,jj)+(hbatt(ji,jj)-rn_hc)*gsi3w3(ji,jj,jk)+rn_hc*zcoeft) 
    1348              END DO 
    1349  
    1350            ENDDO   ! for all jj's 
    1351          ENDDO    ! for all ji's 
    1352  
    1353  
    1354          DO ji=1,jpi 
    1355            DO jj=1,jpj 
    1356  
    1357              DO jk = 1, jpk 
    1358                 esigtu3(ji,jj,jk) = ( hbatt(ji,jj)*esigt3(ji,jj,jk)+hbatt(ji+1,jj)*esigt3(ji+1,jj,jk) ) / & 
    1359                                    ( hbatt(ji,jj)+hbatt(ji+1,jj) ) 
    1360                 esigtv3(ji,jj,jk) = ( hbatt(ji,jj)*esigt3(ji,jj,jk)+hbatt(ji,jj+1)*esigt3(ji,jj+1,jk) ) / & 
    1361                                    ( hbatt(ji,jj)+hbatt(ji,jj+1) ) 
    1362                 esigtf3(ji,jj,jk) = ( hbatt(ji,jj)*esigt3(ji,jj,jk)+hbatt(ji+1,jj)*esigt3(ji+1,jj,jk) +  & 
    1363                                      hbatt(ji,jj+1)*esigt3(ji,jj+1,jk)+hbatt(ji+1,jj+1)*esigt3(ji+1,jj+1,jk) ) / & 
    1364                                    ( hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1) ) 
    1365                 esigwu3(ji,jj,jk) = ( hbatt(ji,jj)*esigw3(ji,jj,jk)+hbatt(ji+1,jj)*esigw3(ji+1,jj,jk) ) / & 
    1366                                    ( hbatt(ji,jj)+hbatt(ji+1,jj) ) 
    1367                 esigwv3(ji,jj,jk) = ( hbatt(ji,jj)*esigw3(ji,jj,jk)+hbatt(ji,jj+1)*esigw3(ji,jj+1,jk) ) / & 
    1368                                    ( hbatt(ji,jj)+hbatt(ji,jj+1) ) 
    1369  
    1370                 e3t(ji,jj,jk)=((hbatt(ji,jj)-rn_hc)*esigt3(ji,jj,jk) + rn_hc/FLOAT(jpkm1)) 
    1371                 e3u(ji,jj,jk)=((hbatu(ji,jj)-rn_hc)*esigtu3(ji,jj,jk) + rn_hc/FLOAT(jpkm1)) 
    1372                 e3v(ji,jj,jk)=((hbatv(ji,jj)-rn_hc)*esigtv3(ji,jj,jk) + rn_hc/FLOAT(jpkm1)) 
    1373                 e3f(ji,jj,jk)=((hbatf(ji,jj)-rn_hc)*esigtf3(ji,jj,jk) + rn_hc/FLOAT(jpkm1)) 
    1374                 ! 
    1375                 e3w (ji,jj,jk)=((hbatt(ji,jj)-rn_hc)*esigw3(ji,jj,jk) + rn_hc/FLOAT(jpkm1)) 
    1376                 e3uw(ji,jj,jk)=((hbatu(ji,jj)-rn_hc)*esigwu3(ji,jj,jk) + rn_hc/FLOAT(jpkm1)) 
    1377                 e3vw(ji,jj,jk)=((hbatv(ji,jj)-rn_hc)*esigwv3(ji,jj,jk) + rn_hc/FLOAT(jpkm1)) 
    1378              END DO 
    1379  
    1380            ENDDO 
    1381          ENDDO 
    1382  
     1396            END DO 
     1397         END DO 
     1398         ! 
    13831399      ELSE   ! not ln_s_sigma 
    1384  
     1400         ! 
    13851401         DO jk = 1, jpk 
    13861402           gsigw(jk) = -fssig( REAL(jk,wp)-0.5_wp ) 
    1387            gsigt(jk) = -fssig( REAL(jk,wp)     ) 
    1388          END DO 
    1389          IF( nprint == 1 .AND. lwp ) WRITE(numout,*) 'gsigw 1 jpk    ', gsigw(1), gsigw(jpk) 
     1403           gsigt(jk) = -fssig( REAL(jk,wp)        ) 
     1404         END DO 
     1405         IF( nprint == 1 .AND. lwp )   WRITE(numout,*) 'gsigw 1 jpk    ', gsigw(1), gsigw(jpk) 
    13901406         ! 
    1391      ! Coefficients for vertical scale factors at w-, t- levels 
     1407         ! Coefficients for vertical scale factors at w-, t- levels 
    13921408!!gm bug :  define it from analytical function, not like juste bellow.... 
    13931409!!gm        or betteroffer the 2 possibilities.... 
    … …  
    13961412            esigw(jk+1) = gsigt(jk+1) - gsigt(jk) 
    13971413         END DO 
    1398          esigw( 1 ) = 2.0_wp * (gsigt(1) - gsigw(1))  
    1399          esigt(jpk) = 2.0_wp * (gsigt(jpk) - gsigw(jpk)) 
     1414         esigw( 1 ) = 2._wp * ( gsigt(1  ) - gsigw(1  ) )  
     1415         esigt(jpk) = 2._wp * ( gsigt(jpk) - gsigw(jpk) ) 
    14001416 
    14011417!!gm  original form 
    … …  
    14071423         ! 
    14081424         ! Coefficients for vertical depth as the sum of e3w scale factors 
    1409          gsi3w(1) = 0.5 * esigw(1) 
     1425         gsi3w(1) = 0.5_wp * esigw(1) 
    14101426         DO jk = 2, jpk 
    14111427            gsi3w(jk) = gsi3w(jk-1) + esigw(jk) 
    … …  
    14131429!!gm: depuw, depvw can be suppressed (modif in ldfslp) and depw=dep3w can be set (save 3 3D arrays) 
    14141430         DO jk = 1, jpk 
    1415             zcoeft = ( FLOAT(jk) - 0.5 ) / FLOAT(jpkm1) 
    1416             zcoefw = ( FLOAT(jk) - 1.0 ) / FLOAT(jpkm1) 
    1417             gdept (:,:,jk) = (scosrf(:,:)+(hbatt(:,:)-hift(:,:))*gsigt(jk)+hift(:,:)*zcoeft) 
    1418             gdepw (:,:,jk) = (scosrf(:,:)+(hbatt(:,:)-hift(:,:))*gsigw(jk)+hift(:,:)*zcoefw) 
    1419             gdep3w(:,:,jk) = (scosrf(:,:)+(hbatt(:,:)-hift(:,:))*gsi3w(jk)+hift(:,:)*zcoeft) 
     1431            zcoeft = ( REAL(jk,wp) - 0.5_wp ) / REAL(jpkm1,wp) 
     1432            zcoefw = ( REAL(jk,wp) - 1.0_wp ) / REAL(jpkm1,wp) 
     1433            gdept (:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*gsigt(jk) + hift(:,:)*zcoeft ) 
     1434            gdepw (:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*gsigw(jk) + hift(:,:)*zcoefw ) 
     1435            gdep3w(:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*gsi3w(jk) + hift(:,:)*zcoeft ) 
    14201436         END DO 
    14211437!!gm: e3uw, e3vw can be suppressed  (modif in dynzdf, dynzdf_iso, zdfbfr) (save 2 3D arrays) 
    … …  
    14231439            DO ji = 1, jpi 
    14241440               DO jk = 1, jpk 
    1425                  e3t(ji,jj,jk)=((hbatt(ji,jj)-hift(ji,jj))*esigt(jk) + hift(ji,jj)/FLOAT(jpkm1)) 
    1426                  e3u(ji,jj,jk)=((hbatu(ji,jj)-hifu(ji,jj))*esigt(jk) + hifu(ji,jj)/FLOAT(jpkm1)) 
    1427                  e3v(ji,jj,jk)=((hbatv(ji,jj)-hifv(ji,jj))*esigt(jk) + hifv(ji,jj)/FLOAT(jpkm1)) 
    1428                  e3f(ji,jj,jk)=((hbatf(ji,jj)-hiff(ji,jj))*esigt(jk) + hiff(ji,jj)/FLOAT(jpkm1)) 
     1441                 e3t(ji,jj,jk) = ( (hbatt(ji,jj)-hift(ji,jj))*esigt(jk) + hift(ji,jj)/REAL(jpkm1,wp) ) 
     1442                 e3u(ji,jj,jk) = ( (hbatu(ji,jj)-hifu(ji,jj))*esigt(jk) + hifu(ji,jj)/REAL(jpkm1,wp) ) 
     1443                 e3v(ji,jj,jk) = ( (hbatv(ji,jj)-hifv(ji,jj))*esigt(jk) + hifv(ji,jj)/REAL(jpkm1,wp) ) 
     1444                 e3f(ji,jj,jk) = ( (hbatf(ji,jj)-hiff(ji,jj))*esigt(jk) + hiff(ji,jj)/REAL(jpkm1,wp) ) 
    14291445                 ! 
    1430                  e3w (ji,jj,jk)=((hbatt(ji,jj)-hift(ji,jj))*esigw(jk) + hift(ji,jj)/FLOAT(jpkm1)) 
    1431                  e3uw(ji,jj,jk)=((hbatu(ji,jj)-hifu(ji,jj))*esigw(jk) + hifu(ji,jj)/FLOAT(jpkm1)) 
    1432                  e3vw(ji,jj,jk)=((hbatv(ji,jj)-hifv(ji,jj))*esigw(jk) + hifv(ji,jj)/FLOAT(jpkm1)) 
     1446                 e3w (ji,jj,jk) = ( (hbatt(ji,jj)-hift(ji,jj))*esigw(jk) + hift(ji,jj)/REAL(jpkm1,wp) ) 
     1447                 e3uw(ji,jj,jk) = ( (hbatu(ji,jj)-hifu(ji,jj))*esigw(jk) + hifu(ji,jj)/REAL(jpkm1,wp) ) 
     1448                 e3vw(ji,jj,jk) = ( (hbatv(ji,jj)-hifv(ji,jj))*esigw(jk) + hifv(ji,jj)/REAL(jpkm1,wp) ) 
    14331449               END DO 
    14341450            END DO 
    14351451         END DO 
    1436  
     1452         ! 
    14371453      ENDIF ! ln_s_sigma 
    14381454 
    … …  
    14571473            DO jk = 1, jpkm1 
    14581474               IF( scobot(ji,jj) >= fsdept(ji,jj,jk) )   mbathy(ji,jj) = MAX( 2, jk ) 
    1459                IF( scobot(ji,jj) == 0.e0             )   mbathy(ji,jj) = 0 
     1475               IF( scobot(ji,jj) == 0._wp            )   mbathy(ji,jj) = 0 
    14601476            END DO 
    14611477         END DO 
    … …  
    14631479      IF( nprint == 1 .AND. lwp ) WRITE(numout,*) ' MIN val mbathy h90 ', MINVAL( mbathy(:,:) ),   & 
    14641480         &                                                       ' MAX ', MAXVAL( mbathy(:,:) ) 
    1465  
    14661481 
    14671482      !                                               ! ============= 
    … …  
    15231538         DO jj = 1, jpj 
    15241539            DO ji = 1, jpi 
    1525                IF( fse3w(ji,jj,jk) <= 0. .OR. fse3t(ji,jj,jk) <= 0. ) THEN 
     1540               IF( fse3w(ji,jj,jk) <= 0._wp .OR. fse3t(ji,jj,jk) <= 0._wp ) THEN 
    15261541                  WRITE(ctmp1,*) 'zgr_sco :   e3w   or e3t   =< 0  at point (i,j,k)= ', ji, jj, jk 
    15271542                  CALL ctl_stop( ctmp1 ) 
    15281543               ENDIF 
    1529                IF( fsdepw(ji,jj,jk) < 0. .OR. fsdept(ji,jj,jk) < 0. ) THEN 
     1544               IF( fsdepw(ji,jj,jk) < 0._wp .OR. fsdept(ji,jj,jk) < 0._wp ) THEN 
    15301545                  WRITE(ctmp1,*) 'zgr_sco :   gdepw or gdept =< 0  at point (i,j,k)= ', ji, jj, jk 
    15311546                  CALL ctl_stop( ctmp1 ) 
    … …  
    15381553   END SUBROUTINE zgr_sco 
    15391554 
    1540 #endif 
    15411555 
    15421556   !!====================================================================== 

  • trunk/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr_substitute.h90

    • Property svn:eol-style deleted
    r1565 r2528  
    88   !!            3.1  !  2009-02  (G. Madec, M. Leclair)  pure z* coordinate 
    99   !!---------------------------------------------------------------------- 
    10 #if defined key_zco 
    11 ! reference for pure z-coordinate (1D - no i,j and time dependency) 
    12 #   define  fsdept_0(i,j,k)  gdept_0(k) 
    13 #   define  fsdepw_0(i,j,k)  gdepw_0(k) 
    14 #   define  fsde3w_0(i,j,k)  gdepw_0(k) 
    15 #   define  fse3t_0(i,j,k)   e3t_0(k) 
    16 #   define  fse3u_0(i,j,k)   e3t_0(k) 
    17 #   define  fse3v_0(i,j,k)   e3t_0(k) 
    18 #   define  fse3f_0(i,j,k)   e3t_0(k) 
    19 #   define  fse3w_0(i,j,k)   e3w_0(k) 
    20 #   define  fse3uw_0(i,j,k)  e3w_0(k) 
    21 #   define  fse3vw_0(i,j,k)  e3w_0(k) 
    22 #else 
    2310! reference for s- or zps-coordinate (3D no time dependency) 
    2411#   define  fsdept_0(i,j,k)  gdept(i,j,k) 
    … …  
    3219#   define  fse3uw_0(i,j,k)  e3uw(i,j,k) 
    3320#   define  fse3vw_0(i,j,k)  e3vw(i,j,k) 
    34 #endif 
    3521#if defined key_vvl 
    3622! s* or z*-coordinate (3D + time dependency) + use of additional now arrays (..._1) 
    … …  
    4632#   define  fse3vw(i,j,k)  e3vw_1(i,j,k) 
    4733 
    48 #   define  fsdept_b(i,j,k)  (fsdept_0(i,j,k)*(1.+sshb(i,j)*mut(i,j,k))) 
    49 #   define  fsdepw_b(i,j,k)  (fsdepw_0(i,j,k)*(1.+sshb(i,j)*mut(i,j,k))) 
    50 #   define  fsde3w_b(i,j,k)  (fsde3w_0(i,j,k)*(1.+sshb(i,j)*mut(i,j,k))-sshb(i,j)) 
    51 #   define  fse3t_b(i,j,k)   (fse3t_0(i,j,k)*(1.+sshb(i,j)*mut(i,j,k))) 
    52 #   define  fse3u_b(i,j,k)   (fse3u_0(i,j,k)*(1.+sshu_b(i,j)*muu(i,j,k))) 
    53 #   define  fse3v_b(i,j,k)   (fse3v_0(i,j,k)*(1.+sshv_b(i,j)*muv(i,j,k))) 
    54 #   define  fse3f_b(i,j,k)   (fse3f_0(i,j,k)*(1.+sshf_b(i,j)*muf(i,j,k))) 
    55 #   define  fse3w_b(i,j,k)   (fse3w_0(i,j,k)*(1.+sshb(i,j)*mut(i,j,k))) 
     34#   define  fse3t_b(i,j,k)   e3t_b(i,j,k) 
     35#   define  fse3u_b(i,j,k)   e3u_b(i,j,k) 
     36#   define  fse3v_b(i,j,k)   e3v_b(i,j,k) 
    5637#   define  fse3uw_b(i,j,k)  (fse3uw_0(i,j,k)*(1.+sshu_b(i,j)*muu(i,j,k))) 
    5738#   define  fse3vw_b(i,j,k)  (fse3vw_0(i,j,k)*(1.+sshv_b(i,j)*muv(i,j,k))) 
    … …  
    7051#   define  fse3t_m(i,j,k)   (fse3t_0(i,j,k)*(1.+ssh_m(i,j)*mut(i,j,k))) 
    7152 
    72 #   define  fsdept_a(i,j,k)  (fsdept_0(i,j,k)*(1.+ssha(i,j)*mut(i,j,k))) 
    73 #   define  fsdepw_a(i,j,k)  (fsdepw_0(i,j,k)*(1.+ssha(i,j)*mut(i,j,k))) 
    74 #   define  fsde3w_a(i,j,k)  (fsde3w_0(i,j,k)*(1.+ssha(i,j)*mut(i,j,k))-ssha(i,j)) 
    7553#   define  fse3t_a(i,j,k)   (fse3t_0(i,j,k)*(1.+ssha(i,j)*mut(i,j,k))) 
    7654#   define  fse3u_a(i,j,k)   (fse3u_0(i,j,k)*(1.+sshu_a(i,j)*muu(i,j,k))) 
    7755#   define  fse3v_a(i,j,k)   (fse3v_0(i,j,k)*(1.+sshv_a(i,j)*muv(i,j,k))) 
    78 #   define  fse3f_a(i,j,k)   (fse3f_0(i,j,k)*(1.+sshf_a(i,j)*muf(i,j,k))) 
    79 #   define  fse3w_a(i,j,k)   (fse3w_0(i,j,k)*(1.+ssha(i,j)*mut(i,j,k))) 
    80 #   define  fse3uw_a(i,j,k)  (fse3uw_0(i,j,k)*(1.+sshu_a(i,j)*muu(i,j,k))) 
    81 #   define  fse3vw_a(i,j,k)  (fse3vw_0(i,j,k)*(1.+sshv_a(i,j)*muv(i,j,k))) 
    8256 
    8357#else 
    … …  
    9468#   define  fse3vw(i,j,k)  fse3vw_0(i,j,k) 
    9569 
    96 #   define  fsdept_b(i,j,k)  fsdept_0(i,j,k) 
    97 #   define  fsdepw_b(i,j,k)  fsdepw_0(i,j,k) 
    98 #   define  fsde3w_b(i,j,k)  fsde3w_0(i,j,k) 
    9970#   define  fse3t_b(i,j,k)   fse3t_0(i,j,k) 
    10071#   define  fse3u_b(i,j,k)   fse3u_0(i,j,k) 
    10172#   define  fse3v_b(i,j,k)   fse3v_0(i,j,k) 
    102 #   define  fse3f_b(i,j,k)   fse3f_0(i,j,k) 
    103 #   define  fse3w_b(i,j,k)   fse3w_0(i,j,k) 
    10473#   define  fse3uw_b(i,j,k)  fse3uw_0(i,j,k) 
    10574#   define  fse3vw_b(i,j,k)  fse3vw_0(i,j,k) 
    … …  
    11887#   define  fse3t_m(i,j,k)   fse3t_0(i,j,k) 
    11988 
    120 #   define  fsdept_a(i,j,k)  fsdept_0(i,j,k) 
    121 #   define  fsdepw_a(i,j,k)  fsdepw_0(i,j,k) 
    122 #   define  fsde3w_a(i,j,k)  fsde3w_0(i,j,k) 
    12389#   define  fse3t_a(i,j,k)   fse3t_0(i,j,k) 
    12490#   define  fse3u_a(i,j,k)   fse3u_0(i,j,k) 
    12591#   define  fse3v_a(i,j,k)   fse3v_0(i,j,k) 
    126 #   define  fse3f_a(i,j,k)   fse3f_0(i,j,k) 
    127 #   define  fse3w_a(i,j,k)   fse3w_0(i,j,k) 
    128 #   define  fse3uw_a(i,j,k)  fse3uw_0(i,j,k) 
    129 #   define  fse3vw_a(i,j,k)  fse3vw_0(i,j,k) 
    13092#endif 
    13193   !!---------------------------------------------------------------------- 
    132    !! NEMO/OPA 3.2 , LOCEAN-IPSL (2009)  
     94   !! NEMO/OPA 3.3 , NEMO Consortium (2010) 
    13395   !! $Id$ 
    134    !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 
     96   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 
    13597   !!---------------------------------------------------------------------- 

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

    • Property svn:eol-style deleted
    r2277 r2528  
    8181   REAL(wp), PUBLIC ::   stefan  =   5.67e-8_wp   !: Stefan-Boltzmann constant  
    8282   !!---------------------------------------------------------------------- 
    83    !! NEMO/OPA 3.2 , LOCEAN-IPSL (2009)  
     83   !! NEMO/OPA 3.3 , NEMO Consortium (2010) 
    8484   !! $Id$  
    85    !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 
     85   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 
    8686   !!---------------------------------------------------------------------- 
    8787    
Note: See TracChangeset for help on using the changeset viewer.