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 6827 – NEMO

Changeset 6827


Ignore:
Timestamp:
2016-08-01T15:37:15+02:00 (8 years ago)
Author:
flavoni
Message:

#1692 - branch SIMPLIF_2_usrdef: commit routines Fortran to create domain_cfg.nc file, to have domain (input) files for new simplification branch. To be moved in TOOLS directory

Location:
branches/2016/dev_r6409_SIMPLIF_2_usrdef_tools/NEMOGCM/NEMO/OPA_SRC
Files:
21 deleted
13 edited

Legend:

Unmodified
Added
Removed
  • branches/2016/dev_r6409_SIMPLIF_2_usrdef_tools/NEMOGCM/NEMO/OPA_SRC/DOM/closea.F90

    r5836 r6827  
    2121   USE phycst          ! physical constants 
    2222   USE in_out_manager  ! I/O manager 
    23    USE sbc_oce         ! ocean surface boundary conditions 
    2423   USE lib_fortran,    ONLY: glob_sum, DDPDD 
    2524   USE lbclnk          ! lateral boundary condition - MPP exchanges 
     
    3130 
    3231   PUBLIC dom_clo      ! routine called by domain module 
    33    PUBLIC sbc_clo      ! routine called by step module 
    34    PUBLIC clo_rnf      ! routine called by sbcrnf module 
    35    PUBLIC clo_ups      ! routine called in traadv_cen2(_jki) module 
    3632   PUBLIC clo_bat      ! routine called in domzgr module 
    3733 
     
    185181 
    186182 
    187    SUBROUTINE sbc_clo( kt ) 
    188       !!--------------------------------------------------------------------- 
    189       !!                  ***  ROUTINE sbc_clo  *** 
    190       !!                     
    191       !! ** Purpose :   Special handling of closed seas 
    192       !! 
    193       !! ** Method  :   Water flux is forced to zero over closed sea 
    194       !!      Excess is shared between remaining ocean, or 
    195       !!      put as run-off in open ocean. 
    196       !! 
    197       !! ** Action  :   emp updated surface freshwater fluxes and associated heat content at kt 
    198       !!---------------------------------------------------------------------- 
    199       INTEGER, INTENT(in) ::   kt   ! ocean model time step 
    200       ! 
    201       INTEGER             ::   ji, jj, jc, jn   ! dummy loop indices 
    202       REAL(wp), PARAMETER ::   rsmall = 1.e-20_wp    ! Closed sea correction epsilon 
    203       REAL(wp)            ::   zze2, ztmp, zcorr     !  
    204       REAL(wp)            ::   zcoef, zcoef1         !  
    205       COMPLEX(wp)         ::   ctmp  
    206       REAL(wp), DIMENSION(jpncs) ::   zfwf   ! 1D workspace 
    207       !!---------------------------------------------------------------------- 
    208       ! 
    209       IF( nn_timing == 1 )  CALL timing_start('sbc_clo') 
    210       !                                                   !------------------! 
    211       IF( kt == nit000 ) THEN                             !  Initialisation  ! 
    212          !                                                !------------------! 
    213          IF(lwp) WRITE(numout,*) 
    214          IF(lwp) WRITE(numout,*)'sbc_clo : closed seas ' 
    215          IF(lwp) WRITE(numout,*)'~~~~~~~' 
    216  
    217          surf(:) = 0.e0_wp 
    218          ! 
    219          surf(jpncs+1) = glob_sum( e1e2t(:,:) )   ! surface of the global ocean 
    220          ! 
    221          !                                        ! surface of closed seas  
    222          IF( lk_mpp_rep ) THEN                         ! MPP reproductible calculation 
    223             DO jc = 1, jpncs 
    224                ctmp = CMPLX( 0.e0, 0.e0, wp ) 
    225                DO jj = ncsj1(jc), ncsj2(jc) 
    226                   DO ji = ncsi1(jc), ncsi2(jc) 
    227                      ztmp = e1e2t(ji,jj) * tmask_i(ji,jj) 
    228                      CALL DDPDD( CMPLX( ztmp, 0.e0, wp ), ctmp ) 
    229                   END DO  
    230                END DO  
    231                IF( lk_mpp )   CALL mpp_sum( ctmp ) 
    232                surf(jc) = REAL(ctmp,wp) 
    233             END DO 
    234          ELSE                                          ! Standard calculation            
    235             DO jc = 1, jpncs 
    236                DO jj = ncsj1(jc), ncsj2(jc) 
    237                   DO ji = ncsi1(jc), ncsi2(jc) 
    238                      surf(jc) = surf(jc) + e1e2t(ji,jj) * tmask_i(ji,jj)      ! surface of closed seas 
    239                   END DO  
    240                END DO  
    241             END DO  
    242             IF( lk_mpp )   CALL mpp_sum ( surf, jpncs )       ! mpp: sum over all the global domain 
    243          ENDIF 
    244  
    245          IF(lwp) WRITE(numout,*)'     Closed sea surfaces' 
    246          DO jc = 1, jpncs 
    247             IF(lwp)WRITE(numout,FMT='(1I3,4I4,5X,F16.2)') jc, ncsi1(jc), ncsi2(jc), ncsj1(jc), ncsj2(jc), surf(jc) 
    248          END DO 
    249  
    250          ! jpncs+1 : surface of sea, closed seas excluded 
    251          DO jc = 1, jpncs 
    252             surf(jpncs+1) = surf(jpncs+1) - surf(jc) 
    253          END DO            
    254          ! 
    255       ENDIF 
    256       !                                                   !--------------------! 
    257       !                                                   !  update emp        ! 
    258       zfwf = 0.e0_wp                                      !--------------------! 
    259       IF( lk_mpp_rep ) THEN                         ! MPP reproductible calculation 
    260          DO jc = 1, jpncs 
    261             ctmp = CMPLX( 0.e0, 0.e0, wp ) 
    262             DO jj = ncsj1(jc), ncsj2(jc) 
    263                DO ji = ncsi1(jc), ncsi2(jc) 
    264                   ztmp = e1e2t(ji,jj) * ( emp(ji,jj)-rnf(ji,jj) ) * tmask_i(ji,jj) 
    265                   CALL DDPDD( CMPLX( ztmp, 0.e0, wp ), ctmp ) 
    266                END DO   
    267             END DO  
    268             IF( lk_mpp )   CALL mpp_sum( ctmp ) 
    269             zfwf(jc) = REAL(ctmp,wp) 
    270          END DO 
    271       ELSE                                          ! Standard calculation            
    272          DO jc = 1, jpncs 
    273             DO jj = ncsj1(jc), ncsj2(jc) 
    274                DO ji = ncsi1(jc), ncsi2(jc) 
    275                   zfwf(jc) = zfwf(jc) + e1e2t(ji,jj) * ( emp(ji,jj)-rnf(ji,jj) ) * tmask_i(ji,jj)  
    276                END DO   
    277             END DO  
    278          END DO 
    279          IF( lk_mpp )   CALL mpp_sum ( zfwf(:) , jpncs )       ! mpp: sum over all the global domain 
    280       ENDIF 
    281  
    282       IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN      ! Black Sea case for ORCA_R2 configuration 
    283          zze2    = ( zfwf(3) + zfwf(4) ) * 0.5_wp 
    284          zfwf(3) = zze2 
    285          zfwf(4) = zze2 
    286       ENDIF 
    287  
    288       zcorr = 0._wp 
    289  
    290       DO jc = 1, jpncs 
    291          ! 
    292          ! The following if avoids the redistribution of the round off 
    293          IF ( ABS(zfwf(jc) / surf(jpncs+1) ) > rsmall) THEN 
    294             ! 
    295             IF( ncstt(jc) == 0 ) THEN           ! water/evap excess is shared by all open ocean 
    296                zcoef    = zfwf(jc) / surf(jpncs+1) 
    297                zcoef1   = rcp * zcoef 
    298                emp(:,:) = emp(:,:) + zcoef 
    299                qns(:,:) = qns(:,:) - zcoef1 * sst_m(:,:) 
    300                ! accumulate closed seas correction 
    301                zcorr    = zcorr    + zcoef 
    302                ! 
    303             ELSEIF( ncstt(jc) == 1 ) THEN       ! Excess water in open sea, at outflow location, excess evap shared 
    304                IF ( zfwf(jc) <= 0.e0_wp ) THEN  
    305                    DO jn = 1, ncsnr(jc) 
    306                      ji = mi0(ncsir(jc,jn)) 
    307                      jj = mj0(ncsjr(jc,jn)) ! Location of outflow in open ocean 
    308                      IF (      ji > 1 .AND. ji < jpi   & 
    309                          .AND. jj > 1 .AND. jj < jpj ) THEN  
    310                          zcoef      = zfwf(jc) / ( REAL(ncsnr(jc)) * e1e2t(ji,jj) ) 
    311                          zcoef1     = rcp * zcoef 
    312                          emp(ji,jj) = emp(ji,jj) + zcoef 
    313                          qns(ji,jj) = qns(ji,jj) - zcoef1 * sst_m(ji,jj) 
    314                      ENDIF  
    315                    END DO  
    316                ELSE  
    317                    zcoef    = zfwf(jc) / surf(jpncs+1) 
    318                    zcoef1   = rcp * zcoef 
    319                    emp(:,:) = emp(:,:) + zcoef 
    320                    qns(:,:) = qns(:,:) - zcoef1 * sst_m(:,:) 
    321                    ! accumulate closed seas correction 
    322                    zcorr    = zcorr    + zcoef 
    323                ENDIF 
    324             ELSEIF( ncstt(jc) == 2 ) THEN       ! Excess e-p-r (either sign) goes to open ocean, at outflow location 
    325                DO jn = 1, ncsnr(jc) 
    326                   ji = mi0(ncsir(jc,jn)) 
    327                   jj = mj0(ncsjr(jc,jn)) ! Location of outflow in open ocean 
    328                   IF(      ji > 1 .AND. ji < jpi    & 
    329                      .AND. jj > 1 .AND. jj < jpj ) THEN  
    330                      zcoef      = zfwf(jc) / ( REAL(ncsnr(jc)) *  e1e2t(ji,jj) ) 
    331                      zcoef1     = rcp * zcoef 
    332                      emp(ji,jj) = emp(ji,jj) + zcoef 
    333                      qns(ji,jj) = qns(ji,jj) - zcoef1 * sst_m(ji,jj) 
    334                   ENDIF  
    335                END DO  
    336             ENDIF  
    337             ! 
    338             DO jj = ncsj1(jc), ncsj2(jc) 
    339                DO ji = ncsi1(jc), ncsi2(jc) 
    340                   zcoef      = zfwf(jc) / surf(jc) 
    341                   zcoef1     = rcp * zcoef 
    342                   emp(ji,jj) = emp(ji,jj) - zcoef 
    343                   qns(ji,jj) = qns(ji,jj) + zcoef1 * sst_m(ji,jj) 
    344                END DO   
    345             END DO  
    346             ! 
    347          END IF 
    348       END DO  
    349  
    350       IF ( ABS(zcorr) > rsmall ) THEN      ! remove the global correction from the closed seas 
    351          DO jc = 1, jpncs                  ! only if it is large enough 
    352             DO jj = ncsj1(jc), ncsj2(jc) 
    353                DO ji = ncsi1(jc), ncsi2(jc) 
    354                   emp(ji,jj) = emp(ji,jj) - zcorr 
    355                   qns(ji,jj) = qns(ji,jj) + rcp * zcorr * sst_m(ji,jj) 
    356                END DO   
    357              END DO  
    358           END DO 
    359       ENDIF 
    360       ! 
    361       emp (:,:) = emp (:,:) * tmask(:,:,1) 
    362       ! 
    363       CALL lbc_lnk( emp , 'T', 1._wp ) 
    364       ! 
    365       IF( nn_timing == 1 )  CALL timing_stop('sbc_clo') 
    366       ! 
    367    END SUBROUTINE sbc_clo 
    368  
    369  
    370    SUBROUTINE clo_rnf( p_rnfmsk ) 
    371       !!--------------------------------------------------------------------- 
    372       !!                  ***  ROUTINE sbc_rnf  *** 
    373       !!                     
    374       !! ** Purpose :   allow the treatment of closed sea outflow grid-points 
    375       !!                to be the same as river mouth grid-points 
    376       !! 
    377       !! ** Method  :   set to 1 the runoff mask (mskrnf, see sbcrnf module) 
    378       !!                at the closed sea outflow grid-point. 
    379       !! 
    380       !! ** Action  :   update (p_)mskrnf (set 1 at closed sea outflow) 
    381       !!---------------------------------------------------------------------- 
    382       REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   p_rnfmsk   ! river runoff mask (rnfmsk array) 
    383       ! 
    384       INTEGER  ::   jc, jn, ji, jj      ! dummy loop indices 
    385       !!---------------------------------------------------------------------- 
    386       ! 
    387       DO jc = 1, jpncs 
    388          IF( ncstt(jc) >= 1 ) THEN            ! runoff mask set to 1 at closed sea outflows 
    389              DO jn = 1, 4 
    390                 DO jj =    mj0( ncsjr(jc,jn) ), mj1( ncsjr(jc,jn) ) 
    391                    DO ji = mi0( ncsir(jc,jn) ), mi1( ncsir(jc,jn) ) 
    392                       p_rnfmsk(ji,jj) = MAX( p_rnfmsk(ji,jj), 1.0_wp ) 
    393                    END DO 
    394                 END DO 
    395             END DO  
    396          ENDIF  
    397       END DO  
    398       ! 
    399    END SUBROUTINE clo_rnf 
    400  
    401     
    402    SUBROUTINE clo_ups( p_upsmsk ) 
    403       !!--------------------------------------------------------------------- 
    404       !!                  ***  ROUTINE sbc_rnf  *** 
    405       !!                     
    406       !! ** Purpose :   allow the treatment of closed sea outflow grid-points 
    407       !!                to be the same as river mouth grid-points 
    408       !! 
    409       !! ** Method  :   set to 0.5 the upstream mask (upsmsk, see traadv_cen2  
    410       !!                module) over the closed seas. 
    411       !! 
    412       !! ** Action  :   update (p_)upsmsk (set 0.5 over closed seas) 
    413       !!---------------------------------------------------------------------- 
    414       REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   p_upsmsk   ! upstream mask (upsmsk array) 
    415       ! 
    416       INTEGER  ::   jc, ji, jj      ! dummy loop indices 
    417       !!---------------------------------------------------------------------- 
    418       ! 
    419       DO jc = 1, jpncs 
    420          DO jj = ncsj1(jc), ncsj2(jc) 
    421             DO ji = ncsi1(jc), ncsi2(jc) 
    422                p_upsmsk(ji,jj) = 0.5_wp         ! mixed upstream/centered scheme over closed seas 
    423             END DO  
    424          END DO  
    425        END DO  
    426        ! 
    427    END SUBROUTINE clo_ups 
    428     
    429        
    430183   SUBROUTINE clo_bat( pbat, kbat ) 
    431184      !!--------------------------------------------------------------------- 
  • branches/2016/dev_r6409_SIMPLIF_2_usrdef_tools/NEMOGCM/NEMO/OPA_SRC/DOM/daymod.F90

    r6140 r6827  
    3333   USE ioipsl  , ONLY :   ymds2ju   ! for calendar 
    3434   USE prtctl         ! Print control 
    35    USE trc_oce , ONLY : lk_offline ! offline flag 
    3635   USE timing         ! Timing 
    37    USE restart        ! restart 
    3836 
    3937   IMPLICIT NONE 
     
    8886      ndt05   = NINT(0.5 * rdt  ) 
    8987 
    90       IF( .NOT. lk_offline ) CALL day_rst( nit000, 'READ' ) 
    9188 
    9289      ! set the calandar from ndastp (read in restart file and namelist) 
     
    281278      ENDIF 
    282279 
    283       IF( .NOT. lk_offline ) CALL rst_opn( kt )               ! Open the restart file if needed and control lrst_oce 
    284       IF( lrst_oce         ) CALL day_rst( kt, 'WRITE' )      ! write day restart information 
    285280      ! 
    286281      IF( nn_timing == 1 )  CALL timing_stop('day') 
     
    288283   END SUBROUTINE day 
    289284 
    290  
    291    SUBROUTINE day_rst( kt, cdrw ) 
    292       !!--------------------------------------------------------------------- 
    293       !!                   ***  ROUTINE ts_rst  *** 
    294       !! 
    295       !!  ** Purpose : Read or write calendar in restart file: 
    296       !! 
    297       !!  WRITE(READ) mode: 
    298       !!       kt        : number of time step since the begining of the experiment at the 
    299       !!                   end of the current(previous) run 
    300       !!       adatrj(0) : number of elapsed days since the begining of the experiment at the 
    301       !!                   end of the current(previous) run (REAL -> keep fractions of day) 
    302       !!       ndastp    : date at the end of the current(previous) run (coded as yyyymmdd integer) 
    303       !! 
    304       !!   According to namelist parameter nrstdt, 
    305       !!       nrstdt = 0  no control on the date (nit000 is  arbitrary). 
    306       !!       nrstdt = 1  we verify that nit000 is equal to the last 
    307       !!                   time step of previous run + 1. 
    308       !!       In both those options, the  exact duration of the experiment 
    309       !!       since the beginning (cumulated duration of all previous restart runs) 
    310       !!       is not stored in the restart and is assumed to be (nit000-1)*rdt. 
    311       !!       This is valid is the time step has remained constant. 
    312       !! 
    313       !!       nrstdt = 2  the duration of the experiment in days (adatrj) 
    314       !!                    has been stored in the restart file. 
    315       !!---------------------------------------------------------------------- 
    316       INTEGER         , INTENT(in) ::   kt         ! ocean time-step 
    317       CHARACTER(len=*), INTENT(in) ::   cdrw       ! "READ"/"WRITE" flag 
    318       ! 
    319       REAL(wp) ::   zkt, zndastp, zdayfrac, ksecs, ktime 
    320       INTEGER  ::   ihour, iminute 
    321       !!---------------------------------------------------------------------- 
    322  
    323       IF( TRIM(cdrw) == 'READ' ) THEN 
    324  
    325          IF( iom_varid( numror, 'kt', ldstop = .FALSE. ) > 0 ) THEN 
    326             ! Get Calendar informations 
    327             CALL iom_get( numror, 'kt', zkt )   ! last time-step of previous run 
    328             IF(lwp) THEN 
    329                WRITE(numout,*) ' *** Info read in restart : ' 
    330                WRITE(numout,*) '   previous time-step                               : ', NINT( zkt ) 
    331                WRITE(numout,*) ' *** restart option' 
    332                SELECT CASE ( nrstdt ) 
    333                CASE ( 0 )   ;   WRITE(numout,*) ' nrstdt = 0 : no control of nit000' 
    334                CASE ( 1 )   ;   WRITE(numout,*) ' nrstdt = 1 : no control the date at nit000 (use ndate0 read in the namelist)' 
    335                CASE ( 2 )   ;   WRITE(numout,*) ' nrstdt = 2 : calendar parameters read in restart' 
    336                END SELECT 
    337                WRITE(numout,*) 
    338             ENDIF 
    339             ! Control of date 
    340             IF( nit000 - NINT( zkt ) /= 1 .AND. nrstdt /= 0 )                                         & 
    341                  &   CALL ctl_stop( ' ===>>>> : problem with nit000 for the restart',                 & 
    342                  &                  ' verify the restart file or rerun with nrstdt = 0 (namelist)' ) 
    343             ! define ndastp and adatrj 
    344             IF ( nrstdt == 2 ) THEN 
    345                ! read the parameters corresponding to nit000 - 1 (last time step of previous run) 
    346                CALL iom_get( numror, 'ndastp', zndastp ) 
    347                ndastp = NINT( zndastp ) 
    348                CALL iom_get( numror, 'adatrj', adatrj  ) 
    349           CALL iom_get( numror, 'ntime', ktime ) 
    350           nn_time0=INT(ktime) 
    351                ! calculate start time in hours and minutes 
    352           zdayfrac=adatrj-INT(adatrj) 
    353           ksecs = NINT(zdayfrac*86400)        ! Nearest second to catch rounding errors in adatrj          
    354           ihour = INT(ksecs/3600) 
    355           iminute = ksecs/60-ihour*60 
    356             
    357                ! Add to nn_time0 
    358                nhour   =   nn_time0 / 100 
    359                nminute = ( nn_time0 - nhour * 100 ) 
    360           nminute=nminute+iminute 
    361            
    362           IF( nminute >= 60 ) THEN 
    363              nminute=nminute-60 
    364         nhour=nhour+1 
    365           ENDIF 
    366           nhour=nhour+ihour 
    367           IF( nhour >= 24 ) THEN 
    368         nhour=nhour-24 
    369              adatrj=adatrj+1 
    370           ENDIF           
    371           nn_time0 = nhour * 100 + nminute 
    372           adatrj = INT(adatrj)                    ! adatrj set to integer as nn_time0 updated           
    373             ELSE 
    374                ! parameters corresponding to nit000 - 1 (as we start the step loop with a call to day) 
    375                ndastp = ndate0        ! ndate0 read in the namelist in dom_nam 
    376                nhour   =   nn_time0 / 100 
    377                nminute = ( nn_time0 - nhour * 100 ) 
    378                IF( nhour*3600+nminute*60-ndt05 .lt. 0 )  ndastp=ndastp-1      ! Start hour is specified in the namelist (default 0) 
    379                adatrj = ( REAL( nit000-1, wp ) * rdt ) / rday 
    380                ! note this is wrong if time step has changed during run 
    381             ENDIF 
    382          ELSE 
    383             ! parameters corresponding to nit000 - 1 (as we start the step loop with a call to day) 
    384             ndastp = ndate0           ! ndate0 read in the namelist in dom_nam 
    385             nhour   =   nn_time0 / 100 
    386        nminute = ( nn_time0 - nhour * 100 ) 
    387             IF( nhour*3600+nminute*60-ndt05 .lt. 0 )  ndastp=ndastp-1      ! Start hour is specified in the namelist (default 0) 
    388             adatrj = ( REAL( nit000-1, wp ) * rdt ) / rday 
    389          ENDIF 
    390          IF( ABS(adatrj  - REAL(NINT(adatrj),wp)) < 0.1 / rday )   adatrj = REAL(NINT(adatrj),wp)   ! avoid truncation error 
    391          ! 
    392          IF(lwp) THEN 
    393             WRITE(numout,*) ' *** Info used values : ' 
    394             WRITE(numout,*) '   date ndastp                                      : ', ndastp 
    395             WRITE(numout,*) '   number of elapsed days since the begining of run : ', adatrj 
    396        WRITE(numout,*) '   nn_time0                                         : ',nn_time0 
    397             WRITE(numout,*) 
    398          ENDIF 
    399          ! 
    400       ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN 
    401          ! 
    402          IF( kt == nitrst ) THEN 
    403             IF(lwp) WRITE(numout,*) 
    404             IF(lwp) WRITE(numout,*) 'rst_write : write oce restart file  kt =', kt 
    405             IF(lwp) WRITE(numout,*) '~~~~~~~' 
    406          ENDIF 
    407          ! calendar control 
    408          CALL iom_rstput( kt, nitrst, numrow, 'kt'     , REAL( kt    , wp) )   ! time-step 
    409          CALL iom_rstput( kt, nitrst, numrow, 'ndastp' , REAL( ndastp, wp) )   ! date 
    410          CALL iom_rstput( kt, nitrst, numrow, 'adatrj' , adatrj            )   ! number of elapsed days since 
    411          !                                                                     ! the begining of the run [s] 
    412     CALL iom_rstput( kt, nitrst, numrow, 'ntime'  , REAL( nn_time0, wp) ) ! time 
    413       ENDIF 
    414       ! 
    415    END SUBROUTINE day_rst 
    416285 
    417286   !!====================================================================== 
  • branches/2016/dev_r6409_SIMPLIF_2_usrdef_tools/NEMOGCM/NEMO/OPA_SRC/DOM/dom_oce.F90

    r6140 r6827  
    166166   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   e1e2f , r1_e1e2f                !: associated metrics at f-point 
    167167   ! 
    168    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ff                              !: coriolis factor                   [1/s] 
     168   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ff_f, ff_t                      !: coriolis factor                   [1/s] 
    169169 
    170170   !!---------------------------------------------------------------------- 
     
    332332         &      e1e2v(jpi,jpj) , r1_e1e2v(jpi,jpj) , e1_e2v(jpi,jpj)                   ,     & 
    333333         &      e1e2f(jpi,jpj) , r1_e1e2f(jpi,jpj)                                     ,     & 
    334          &        ff (jpi,jpj)                                                         , STAT=ierr(3) ) 
     334         &       ff_f(jpi,jpj) ,     ff_t(jpi,jpj)                                     , STAT=ierr(3) ) 
    335335         ! 
    336336      ALLOCATE( gdept_0(jpi,jpj,jpk) , gdepw_0(jpi,jpj,jpk) , gde3w_0(jpi,jpj,jpk) ,     & 
  • branches/2016/dev_r6409_SIMPLIF_2_usrdef_tools/NEMOGCM/NEMO/OPA_SRC/DOM/domain.F90

    r6140 r6827  
    2424   USE oce             ! ocean variables 
    2525   USE dom_oce         ! domain: ocean 
    26    USE sbc_oce         ! surface boundary condition: ocean 
    2726   USE phycst          ! physical constants 
    2827   USE closea          ! closed seas 
     
    3332   USE domwri          ! domain: write the meshmask file 
    3433   USE domvvl          ! variable volume 
    35    USE c1d             ! 1D vertical configuration 
    36    USE dyncor_c1d      ! Coriolis term (c1d case)         (cor_c1d routine) 
    3734   ! 
    3835   USE in_out_manager  ! I/O manager 
     36   USE iom             !  
    3937   USE wrk_nemo        ! Memory Allocation 
    4038   USE lib_mpp         ! distributed memory computing library 
     
    137135      ENDIF 
    138136      ! 
    139       IF( lk_c1d         )   CALL cor_c1d       ! 1D configuration: Coriolis set at T-point 
    140       ! 
    141                              CALL dom_stp       ! time step 
    142       IF( nmsh /= 0 .AND. .NOT. ln_iscpl )                         CALL dom_wri      ! Create a domain file 
    143       IF( nmsh /= 0 .AND.       ln_iscpl .AND. .NOT. ln_rstart )   CALL dom_wri      ! Create a domain file 
    144       IF( .NOT.ln_rstart )   CALL dom_ctl       ! Domain control 
     137      CALL cfg_write         ! create the configuration file 
    145138      ! 
    146139      IF( nn_timing == 1 )   CALL timing_stop('dom_init') 
     
    465458   END SUBROUTINE dom_stiff 
    466459 
     460 
     461   SUBROUTINE cfg_write 
     462      !!---------------------------------------------------------------------- 
     463      !!                  ***  ROUTINE cfg_write  *** 
     464      !!                    
     465      !! ** Purpose :   Create the "domain_cfg" file, a NetCDF file which  
     466      !!              contains all the ocean domain informations required to  
     467      !!              define an ocean configuration. 
     468      !! 
     469      !! ** Method  :   Write in a file all the arrays required to set up an 
     470      !!              ocean configuration. 
     471      !! 
     472      !! ** output file :   domain_cfg.nc : domain size, characteristics, 
     473      !horizontal mesh, 
     474      !!                              Coriolis parameter, depth and vertical 
     475      !scale factors 
     476      !!---------------------------------------------------------------------- 
     477      INTEGER           ::   ji, jj, jk   ! dummy loop indices 
     478      INTEGER           ::   izco, izps, isco, icav 
     479      INTEGER           ::   inum     ! temprary units for 'domain_cfg.nc' file 
     480      CHARACTER(len=21) ::   clnam    ! filename (mesh and mask informations) 
     481      REAL(wp), DIMENSION(jpi,jpj) ::   z2d   ! workspace 
     482      !!---------------------------------------------------------------------- 
     483      ! 
     484      IF(lwp) WRITE(numout,*) 
     485      IF(lwp) WRITE(numout,*) 'cfg_write : create the "domain_cfg.nc" file containing all required configuration information' 
     486      IF(lwp) WRITE(numout,*) '~~~~~~~~~' 
     487      ! 
     488      !                       ! ============================= ! 
     489      !                       !  create 'domain_cfg.nc' file  ! 
     490      !                       ! ============================= ! 
     491      !          
     492      clnam = 'domain_cfg'  ! filename (configuration information) 
     493      CALL iom_open( TRIM(clnam), inum, ldwrt = .TRUE., kiolib = jprstlib ) 
     494       
     495      !                             !==  global domain size  ==! 
     496      ! 
     497      CALL iom_rstput( 0, 0, inum, 'jpiglo', REAL( jpiglo, wp), ktype = jp_i4 ) 
     498      CALL iom_rstput( 0, 0, inum, 'jpjglo', REAL( jpjglo, wp), ktype = jp_i4 ) 
     499      CALL iom_rstput( 0, 0, inum, 'jpkglo', REAL( jpk   , wp), ktype = jp_i4 ) 
     500      ! 
     501      !                             !==  domain characteristics  ==! 
     502      ! 
     503      !                                   ! lateral boundary of the global 
     504      !                                   domain 
     505      CALL iom_rstput( 0, 0, inum, 'jperio', REAL( jperio, wp), ktype = jp_i4 ) 
     506      ! 
     507      !                                   ! type of vertical coordinate 
     508      IF( ln_zco    ) THEN   ;   izco = 1   ;   ELSE   ;   izco = 0   ;   ENDIF 
     509      IF( ln_zps    ) THEN   ;   izps = 1   ;   ELSE   ;   izps = 0   ;   ENDIF 
     510      IF( ln_sco    ) THEN   ;   isco = 1   ;   ELSE   ;   isco = 0   ;   ENDIF 
     511      CALL iom_rstput( 0, 0, inum, 'ln_zco'   , REAL( izco, wp), ktype = jp_i4 ) 
     512      CALL iom_rstput( 0, 0, inum, 'ln_zps'   , REAL( izps, wp), ktype = jp_i4 ) 
     513      CALL iom_rstput( 0, 0, inum, 'ln_sco'   , REAL( isco, wp), ktype = jp_i4 ) 
     514      ! 
     515      !                                   ! ocean cavities under iceshelves 
     516      IF( ln_isfcav ) THEN   ;   icav = 1   ;   ELSE   ;   icav = 0   ;   ENDIF 
     517      CALL iom_rstput( 0, 0, inum, 'ln_isfcav', REAL( icav, wp), ktype = jp_i4 ) 
     518      ! 
     519      !                             !==  horizontal mesh  ! 
     520      ! 
     521      CALL iom_rstput( 0, 0, inum, 'glamt', glamt, ktype = jp_r8 )   ! latitude 
     522      CALL iom_rstput( 0, 0, inum, 'glamu', glamu, ktype = jp_r8 ) 
     523      CALL iom_rstput( 0, 0, inum, 'glamv', glamv, ktype = jp_r8 ) 
     524      CALL iom_rstput( 0, 0, inum, 'glamf', glamf, ktype = jp_r8 ) 
     525      !                                 
     526      CALL iom_rstput( 0, 0, inum, 'gphit', gphit, ktype = jp_r8 )   ! longitude 
     527      CALL iom_rstput( 0, 0, inum, 'gphiu', gphiu, ktype = jp_r8 ) 
     528      CALL iom_rstput( 0, 0, inum, 'gphiv', gphiv, ktype = jp_r8 ) 
     529      CALL iom_rstput( 0, 0, inum, 'gphif', gphif, ktype = jp_r8 ) 
     530      !                                 
     531      CALL iom_rstput( 0, 0, inum, 'e1t'  , e1t  , ktype = jp_r8 )   ! i-scale factors (e1.) 
     532      CALL iom_rstput( 0, 0, inum, 'e1u'  , e1u  , ktype = jp_r8 ) 
     533      CALL iom_rstput( 0, 0, inum, 'e1v'  , e1v  , ktype = jp_r8 ) 
     534      CALL iom_rstput( 0, 0, inum, 'e1f'  , e1f  , ktype = jp_r8 ) 
     535      ! 
     536      CALL iom_rstput( 0, 0, inum, 'e2t'  , e2t  , ktype = jp_r8 )   ! j-scale factors (e2.) 
     537      CALL iom_rstput( 0, 0, inum, 'e2u'  , e2u  , ktype = jp_r8 ) 
     538      CALL iom_rstput( 0, 0, inum, 'e2v'  , e2v  , ktype = jp_r8 ) 
     539      CALL iom_rstput( 0, 0, inum, 'e2f'  , e2f  , ktype = jp_r8 ) 
     540      ! 
     541      CALL iom_rstput( 0, 0, inum, 'ff_f' , ff_f , ktype = jp_r8 )   ! coriolis factor 
     542      CALL iom_rstput( 0, 0, inum, 'ff_t' , ff_t , ktype = jp_r8 ) 
     543      ! 
     544      !                             !==  vertical mesh - 3D mask  ==! 
     545      !                                                      
     546      CALL iom_rstput( 0, 0, inum, 'gdept_1d', gdept_1d, ktype = jp_r8 )   !  reference 1D-coordinate 
     547      CALL iom_rstput( 0, 0, inum, 'gdepw_1d', gdepw_1d, ktype = jp_r8 ) 
     548      CALL iom_rstput( 0, 0, inum, 'e3t_1d'  , e3t_1d  , ktype = jp_r8 ) 
     549      CALL iom_rstput( 0, 0, inum, 'e3w_1d'  , e3w_1d  , ktype = jp_r8 ) 
     550      !                                                      
     551      CALL iom_rstput( 0, 0, inum, 'gdept_0' , gdept_0 , ktype = jp_r8 )   !  depth (t- & w-points) 
     552      CALL iom_rstput( 0, 0, inum, 'gdepw_0' , gdepw_0 , ktype = jp_r8 ) 
     553      ! 
     554      CALL iom_rstput( 0, 0, inum, 'e3t_0'   , e3t_0   , ktype = jp_r8 )   !  vertical scale factors (e 
     555      CALL iom_rstput( 0, 0, inum, 'e3u_0'   , e3u_0   , ktype = jp_r8 ) 
     556      CALL iom_rstput( 0, 0, inum, 'e3v_0'   , e3v_0   , ktype = jp_r8 ) 
     557      CALL iom_rstput( 0, 0, inum, 'e3f_0'   , e3f_0   , ktype = jp_r8 ) 
     558      CALL iom_rstput( 0, 0, inum, 'e3w_0'   , e3w_0   , ktype = jp_r8 ) 
     559      CALL iom_rstput( 0, 0, inum, 'e3uw_0'  , e3uw_0  , ktype = jp_r8 ) 
     560      CALL iom_rstput( 0, 0, inum, 'e3vw_0'  , e3vw_0  , ktype = jp_r8 ) 
     561      !                                          
     562      !                             !==  ocean top and bottom level  ==! 
     563      ! 
     564      CALL iom_rstput( 0, 0, inum, 'bottom level' , REAL( mbkt, wp )*ssmask , ktype = jp_i4 )   ! nb of ocean T-points 
     565      CALL iom_rstput( 0, 0, inum, 'top    level' , REAL( mikt, wp )*ssmask , ktype = jp_i4 )   ! nb of ocean T-points (ISF) 
     566      ! 
     567      IF( ln_sco ) THEN             ! s-coordinate: store grid stiffness ratio (Not required anyway) 
     568         CALL dom_stiff 
     569         !SF  CALL dom_stiff( z2d ) !commented because at compilation error:  The 
     570                                   !number of actual arguments cannot be greater than the number of dummy 
     571                                   !arguments.   [DOM_STIFF] 
     572                                   !CALL dom_stiff( z2d ) 
     573                                   ! --------------^ compilation aborted for 
     574                                   !/workgpfs/rech/gzi/rgzi011/commit-simplif2-TOOLS/NEMOGCM/CONFIG/TEST/BLD/ppsrc/nemo/domain.f90 
     575         CALL iom_rstput( 0, 0, inum, 'stiffness', z2d )        !    ! Max. grid stiffness ratio 
     576      ENDIF 
     577      ! 
     578      !                                ! ============================ 
     579      !                                !        close the files  
     580      !                                ! ============================ 
     581      CALL iom_close( inum ) 
     582      ! 
     583   END SUBROUTINE cfg_write 
     584 
     585 
     586 
    467587   !!====================================================================== 
    468588END MODULE domain 
  • branches/2016/dev_r6409_SIMPLIF_2_usrdef_tools/NEMOGCM/NEMO/OPA_SRC/DOM/domcfg.F90

    r6140 r6827  
    1616   USE lib_mpp         ! distributed memory computing library 
    1717   USE timing          ! Timing 
    18    USE c1d             ! 1D configuration 
    19    USE domc1d          ! 1D configuration: column location 
    2018 
    2119   IMPLICIT NONE 
     
    8280      !!---------------------------------------------------------------------- 
    8381      !                              ! recalculate jpizoom/jpjzoom given lat/lon 
    84       IF( lk_c1d .AND. ln_c1d_locpt )  CALL dom_c1d( rn_lat1d, rn_lon1d ) 
    8582      ! 
    8683      !                        ! ============== ! 
  • branches/2016/dev_r6409_SIMPLIF_2_usrdef_tools/NEMOGCM/NEMO/OPA_SRC/DOM/domhgr.F90

    r6140 r6827  
    377377      CASE ( 0, 1, 4 )               ! mesh on the sphere 
    378378         ! 
    379          ff(:,:) = 2. * omega * SIN( rad * gphif(:,:) )  
     379         ff_f(:,:) = 2. * omega * SIN( rad * gphif(:,:) )  
     380         ff_t(:,:) = 2. * omega * SIN( rad * gphit(:,:) )     !    -        -       - at t-point 
    380381         ! 
    381382      CASE ( 2 )                     ! f-plane at ppgphi0  
    382383         ! 
    383          ff(:,:) = 2. * omega * SIN( rad * ppgphi0 ) 
    384          ! 
    385          IF(lwp) WRITE(numout,*) '          f-plane: Coriolis parameter = constant = ', ff(1,1) 
     384         ff_f(:,:) = 2. * omega * SIN( rad * ppgphi0 ) 
     385         ff_t(:,:) = 2. * omega * SIN( rad * ppgphi0 ) 
     386         ! 
     387         IF(lwp) WRITE(numout,*) '          f-plane: Coriolis parameter = constant = ', ff_f(1,1) 
    386388         ! 
    387389      CASE ( 3 )                     ! beta-plane 
     
    399401         zf0     = 2. * omega * SIN( rad * zphi0 )                              ! compute f0 1st point south 
    400402         ! 
    401          ff(:,:) = ( zf0  + zbeta * gphif(:,:) * 1.e+3 )                        ! f = f0 +beta* y ( y=0 at south) 
     403         ff_f(:,:) = ( zf0  + zbeta * gphif(:,:) * 1.e+3 )                        ! f = f0 +beta* y ( y=0 at south) 
     404         ff_t(:,:) = ( zf0  + zbeta * gphit(:,:) * 1.e+3 )                        ! f = f0 +beta* y ( y=0 at south) 
    402405         ! 
    403406         IF(lwp) THEN 
    404407            WRITE(numout,*)  
    405             WRITE(numout,*) '          Beta-plane: Beta parameter = constant = ', ff(nldi,nldj) 
    406             WRITE(numout,*) '          Coriolis parameter varies from ', ff(nldi,nldj),' to ', ff(nldi,nlej) 
     408            WRITE(numout,*) '          Beta-plane: Beta parameter = constant = ', ff_f(nldi,nldj) 
     409            WRITE(numout,*) '          Coriolis parameter varies from ', ff_f(nldi,nldj),' to ', ff_f(nldi,nlej) 
    407410         ENDIF 
    408411         IF( lk_mpp ) THEN  
    409             zminff=ff(nldi,nldj) 
    410             zmaxff=ff(nldi,nlej) 
     412            zminff=ff_f(nldi,nldj) 
     413            zmaxff=ff_f(nldi,nlej) 
    411414            CALL mpp_min( zminff )   ! min over the global domain 
    412415            CALL mpp_max( zmaxff )   ! max over the global domain 
     
    420423         zf0   = 2. * omega * SIN( rad * zphi0 )                            ! compute f0 1st point south 
    421424         ! 
    422          ff(:,:) = ( zf0 + zbeta * ABS( gphif(:,:) - zphi0 ) * rad * ra )   ! f = f0 +beta* y ( y=0 at south) 
     425         ff_f(:,:) = ( zf0 + zbeta * ABS( gphif(:,:) - zphi0 ) * rad * ra )   ! f = f0 +beta* y ( y=0 at south) 
     426         ff_t(:,:) = ( zf0 + zbeta * ABS( gphit(:,:) - zphi0 ) * rad * ra )   ! f = f0 +beta* y ( y=0 at south) 
    423427         ! 
    424428         IF(lwp) THEN 
    425429            WRITE(numout,*)  
    426430            WRITE(numout,*) '          Beta-plane and rotated domain : ' 
    427             WRITE(numout,*) '          Coriolis parameter varies in this processor from ', ff(nldi,nldj),' to ', ff(nldi,nlej) 
     431            WRITE(numout,*) '          Coriolis parameter varies in this processor from ', ff_f(nldi,nldj),' to ', ff_f(nldi,nlej) 
    428432         ENDIF 
    429433         ! 
    430434         IF( lk_mpp ) THEN  
    431             zminff=ff(nldi,nldj) 
    432             zmaxff=ff(nldi,nlej) 
     435            zminff=ff_f(nldi,nldj) 
     436            zmaxff=ff_f(nldi,nlej) 
    433437            CALL mpp_min( zminff )   ! min over the global domain 
    434438            CALL mpp_max( zmaxff )   ! max over the global domain 
  • branches/2016/dev_r6409_SIMPLIF_2_usrdef_tools/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90

    r6351 r6827  
    2222   USE phycst          ! physical constant 
    2323   USE dom_oce         ! ocean space and time domain 
    24    USE sbc_oce         ! ocean surface boundary condition 
    25    USE wet_dry         ! wetting and drying 
    26    USE restart         ! ocean restart 
    2724   ! 
    2825   USE in_out_manager  ! I/O manager 
     
    3734 
    3835   PUBLIC  dom_vvl_init       ! called by domain.F90 
    39    PUBLIC  dom_vvl_sf_nxt     ! called by step.F90 
    40    PUBLIC  dom_vvl_sf_swp     ! called by step.F90 
    41    PUBLIC  dom_vvl_interpol   ! called by dynnxt.F90 
    4236 
    4337   !                                                      !!* Namelist nam_vvl 
     
    133127      ! 
    134128      !                    ! Read or initialize e3t_(b/n), tilde_e3t_(b/n) and hdiv_lf 
    135       CALL dom_vvl_rst( nit000, 'READ' ) 
    136129      e3t_a(:,:,jpk) = e3t_0(:,:,jpk)  ! last level always inside the sea floor set one for all 
    137130      ! 
     
    246239 
    247240 
    248    SUBROUTINE dom_vvl_sf_nxt( kt, kcall )  
    249       !!---------------------------------------------------------------------- 
    250       !!                ***  ROUTINE dom_vvl_sf_nxt  *** 
    251       !!                    
    252       !! ** Purpose :  - compute the after scale factors used in tra_zdf, dynnxt, 
    253       !!                 tranxt and dynspg routines 
    254       !! 
    255       !! ** Method  :  - z_star case:  Repartition of ssh INCREMENT proportionnaly to the level thickness. 
    256       !!               - z_tilde_case: after scale factor increment =  
    257       !!                                    high frequency part of horizontal divergence 
    258       !!                                  + retsoring towards the background grid 
    259       !!                                  + thickness difusion 
    260       !!                               Then repartition of ssh INCREMENT proportionnaly 
    261       !!                               to the "baroclinic" level thickness. 
    262       !! 
    263       !! ** Action  :  - hdiv_lf    : restoring towards full baroclinic divergence in z_tilde case 
    264       !!               - tilde_e3t_a: after increment of vertical scale factor  
    265       !!                              in z_tilde case 
    266       !!               - e3(t/u/v)_a 
    267       !! 
    268       !! Reference  : Leclair, M., and Madec, G. 2011, Ocean Modelling. 
    269       !!---------------------------------------------------------------------- 
    270       INTEGER, INTENT( in )           ::   kt      ! time step 
    271       INTEGER, INTENT( in ), OPTIONAL ::   kcall   ! optional argument indicating call sequence 
    272       ! 
    273       INTEGER                ::   ji, jj, jk            ! dummy loop indices 
    274       INTEGER , DIMENSION(3) ::   ijk_max, ijk_min      ! temporary integers 
    275       REAL(wp)               ::   z2dt, z_tmin, z_tmax  ! local scalars 
    276       LOGICAL                ::   ll_do_bclinic         ! local logical 
    277       REAL(wp), POINTER, DIMENSION(:,:,:) ::   ze3t 
    278       REAL(wp), POINTER, DIMENSION(:,:  ) ::   zht, z_scale, zwu, zwv, zhdiv 
    279       !!---------------------------------------------------------------------- 
    280       ! 
    281       IF( ln_linssh )   RETURN      ! No calculation in linear free surface 
    282       ! 
    283       IF( nn_timing == 1 )   CALL timing_start('dom_vvl_sf_nxt') 
    284       ! 
    285       CALL wrk_alloc( jpi,jpj,zht,   z_scale, zwu, zwv, zhdiv ) 
    286       CALL wrk_alloc( jpi,jpj,jpk,   ze3t ) 
    287  
    288       IF( kt == nit000 ) THEN 
    289          IF(lwp) WRITE(numout,*) 
    290          IF(lwp) WRITE(numout,*) 'dom_vvl_sf_nxt : compute after scale factors' 
    291          IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~' 
    292       ENDIF 
    293  
    294       ll_do_bclinic = .TRUE. 
    295       IF( PRESENT(kcall) ) THEN 
    296          IF( kcall == 2 .AND. ln_vvl_ztilde )   ll_do_bclinic = .FALSE. 
    297       ENDIF 
    298  
    299       ! ******************************* ! 
    300       ! After acale factors at t-points ! 
    301       ! ******************************* ! 
    302       !                                                ! --------------------------------------------- ! 
    303       !                                                ! z_star coordinate and barotropic z-tilde part ! 
    304       !                                                ! --------------------------------------------- ! 
    305       ! 
    306       z_scale(:,:) = ( ssha(:,:) - sshb(:,:) ) * ssmask(:,:) / ( ht_0(:,:) + sshn(:,:) + 1. - ssmask(:,:) ) 
    307       DO jk = 1, jpkm1 
    308          ! formally this is the same as e3t_a = e3t_0*(1+ssha/ht_0) 
    309          e3t_a(:,:,jk) = e3t_b(:,:,jk) + e3t_n(:,:,jk) * z_scale(:,:) * tmask(:,:,jk) 
    310       END DO 
    311       ! 
    312       IF( ln_vvl_ztilde .OR. ln_vvl_layer .AND. ll_do_bclinic ) THEN   ! z_tilde or layer coordinate ! 
    313          !                                                            ! ------baroclinic part------ ! 
    314          ! I - initialization 
    315          ! ================== 
    316  
    317          ! 1 - barotropic divergence 
    318          ! ------------------------- 
    319          zhdiv(:,:) = 0._wp 
    320          zht(:,:)   = 0._wp 
    321          DO jk = 1, jpkm1 
    322             zhdiv(:,:) = zhdiv(:,:) + e3t_n(:,:,jk) * hdivn(:,:,jk) 
    323             zht  (:,:) = zht  (:,:) + e3t_n(:,:,jk) * tmask(:,:,jk) 
    324          END DO 
    325          zhdiv(:,:) = zhdiv(:,:) / ( zht(:,:) + 1. - tmask_i(:,:) ) 
    326  
    327          ! 2 - Low frequency baroclinic horizontal divergence  (z-tilde case only) 
    328          ! -------------------------------------------------- 
    329          IF( ln_vvl_ztilde ) THEN 
    330             IF( kt > nit000 ) THEN 
    331                DO jk = 1, jpkm1 
    332                   hdiv_lf(:,:,jk) = hdiv_lf(:,:,jk) - rdt * frq_rst_hdv(:,:)   & 
    333                      &          * ( hdiv_lf(:,:,jk) - e3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) ) 
    334                END DO 
    335             ENDIF 
    336          ENDIF 
    337  
    338          ! II - after z_tilde increments of vertical scale factors 
    339          ! ======================================================= 
    340          tilde_e3t_a(:,:,:) = 0._wp  ! tilde_e3t_a used to store tendency terms 
    341  
    342          ! 1 - High frequency divergence term 
    343          ! ---------------------------------- 
    344          IF( ln_vvl_ztilde ) THEN     ! z_tilde case 
    345             DO jk = 1, jpkm1 
    346                tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - ( e3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) - hdiv_lf(:,:,jk) ) 
    347             END DO 
    348          ELSE                         ! layer case 
    349             DO jk = 1, jpkm1 
    350                tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) -   e3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) * tmask(:,:,jk) 
    351             END DO 
    352          ENDIF 
    353  
    354          ! 2 - Restoring term (z-tilde case only) 
    355          ! ------------------ 
    356          IF( ln_vvl_ztilde ) THEN 
    357             DO jk = 1, jpk 
    358                tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - frq_rst_e3t(:,:) * tilde_e3t_b(:,:,jk) 
    359             END DO 
    360          ENDIF 
    361  
    362          ! 3 - Thickness diffusion term 
    363          ! ---------------------------- 
    364          zwu(:,:) = 0._wp 
    365          zwv(:,:) = 0._wp 
    366          DO jk = 1, jpkm1        ! a - first derivative: diffusive fluxes 
    367             DO jj = 1, jpjm1 
    368                DO ji = 1, fs_jpim1   ! vector opt. 
    369                   un_td(ji,jj,jk) = rn_ahe3 * umask(ji,jj,jk) * e2_e1u(ji,jj)           & 
    370                      &            * ( tilde_e3t_b(ji,jj,jk) - tilde_e3t_b(ji+1,jj  ,jk) ) 
    371                   vn_td(ji,jj,jk) = rn_ahe3 * vmask(ji,jj,jk) * e1_e2v(ji,jj)           &  
    372                      &            * ( tilde_e3t_b(ji,jj,jk) - tilde_e3t_b(ji  ,jj+1,jk) ) 
    373                   zwu(ji,jj) = zwu(ji,jj) + un_td(ji,jj,jk) 
    374                   zwv(ji,jj) = zwv(ji,jj) + vn_td(ji,jj,jk) 
    375                END DO 
    376             END DO 
    377          END DO 
    378          DO jj = 1, jpj          ! b - correction for last oceanic u-v points 
    379             DO ji = 1, jpi 
    380                un_td(ji,jj,mbku(ji,jj)) = un_td(ji,jj,mbku(ji,jj)) - zwu(ji,jj) 
    381                vn_td(ji,jj,mbkv(ji,jj)) = vn_td(ji,jj,mbkv(ji,jj)) - zwv(ji,jj) 
    382             END DO 
    383          END DO 
    384          DO jk = 1, jpkm1        ! c - second derivative: divergence of diffusive fluxes 
    385             DO jj = 2, jpjm1 
    386                DO ji = fs_2, fs_jpim1   ! vector opt. 
    387                   tilde_e3t_a(ji,jj,jk) = tilde_e3t_a(ji,jj,jk) + (   un_td(ji-1,jj  ,jk) - un_td(ji,jj,jk)    & 
    388                      &                                          +     vn_td(ji  ,jj-1,jk) - vn_td(ji,jj,jk)    & 
    389                      &                                            ) * r1_e1e2t(ji,jj) 
    390                END DO 
    391             END DO 
    392          END DO 
    393          !                       ! d - thickness diffusion transport: boundary conditions 
    394          !                             (stored for tracer advction and continuity equation) 
    395          CALL lbc_lnk( un_td , 'U' , -1._wp) 
    396          CALL lbc_lnk( vn_td , 'V' , -1._wp) 
    397  
    398          ! 4 - Time stepping of baroclinic scale factors 
    399          ! --------------------------------------------- 
    400          ! Leapfrog time stepping 
    401          ! ~~~~~~~~~~~~~~~~~~~~~~ 
    402          IF( neuler == 0 .AND. kt == nit000 ) THEN 
    403             z2dt =  rdt 
    404          ELSE 
    405             z2dt = 2.0_wp * rdt 
    406          ENDIF 
    407          CALL lbc_lnk( tilde_e3t_a(:,:,:), 'T', 1._wp ) 
    408          tilde_e3t_a(:,:,:) = tilde_e3t_b(:,:,:) + z2dt * tmask(:,:,:) * tilde_e3t_a(:,:,:) 
    409  
    410          ! Maximum deformation control 
    411          ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~ 
    412          ze3t(:,:,jpk) = 0._wp 
    413          DO jk = 1, jpkm1 
    414             ze3t(:,:,jk) = tilde_e3t_a(:,:,jk) / e3t_0(:,:,jk) * tmask(:,:,jk) * tmask_i(:,:) 
    415          END DO 
    416          z_tmax = MAXVAL( ze3t(:,:,:) ) 
    417          IF( lk_mpp )   CALL mpp_max( z_tmax )                 ! max over the global domain 
    418          z_tmin = MINVAL( ze3t(:,:,:) ) 
    419          IF( lk_mpp )   CALL mpp_min( z_tmin )                 ! min over the global domain 
    420          ! - ML - test: for the moment, stop simulation for too large e3_t variations 
    421          IF( ( z_tmax >  rn_zdef_max ) .OR. ( z_tmin < - rn_zdef_max ) ) THEN 
    422             IF( lk_mpp ) THEN 
    423                CALL mpp_maxloc( ze3t, tmask, z_tmax, ijk_max(1), ijk_max(2), ijk_max(3) ) 
    424                CALL mpp_minloc( ze3t, tmask, z_tmin, ijk_min(1), ijk_min(2), ijk_min(3) ) 
    425             ELSE 
    426                ijk_max = MAXLOC( ze3t(:,:,:) ) 
    427                ijk_max(1) = ijk_max(1) + nimpp - 1 
    428                ijk_max(2) = ijk_max(2) + njmpp - 1 
    429                ijk_min = MINLOC( ze3t(:,:,:) ) 
    430                ijk_min(1) = ijk_min(1) + nimpp - 1 
    431                ijk_min(2) = ijk_min(2) + njmpp - 1 
    432             ENDIF 
    433             IF (lwp) THEN 
    434                WRITE(numout, *) 'MAX( tilde_e3t_a(:,:,:) / e3t_0(:,:,:) ) =', z_tmax 
    435                WRITE(numout, *) 'at i, j, k=', ijk_max 
    436                WRITE(numout, *) 'MIN( tilde_e3t_a(:,:,:) / e3t_0(:,:,:) ) =', z_tmin 
    437                WRITE(numout, *) 'at i, j, k=', ijk_min             
    438                CALL ctl_warn('MAX( ABS( tilde_e3t_a(:,:,:) ) / e3t_0(:,:,:) ) too high') 
    439             ENDIF 
    440          ENDIF 
    441          ! - ML - end test 
    442          ! - ML - Imposing these limits will cause a baroclinicity error which is corrected for below 
    443          tilde_e3t_a(:,:,:) = MIN( tilde_e3t_a(:,:,:),   rn_zdef_max * e3t_0(:,:,:) ) 
    444          tilde_e3t_a(:,:,:) = MAX( tilde_e3t_a(:,:,:), - rn_zdef_max * e3t_0(:,:,:) ) 
    445  
    446          ! 
    447          ! "tilda" change in the after scale factor 
    448          ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 
    449          DO jk = 1, jpkm1 
    450             dtilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - tilde_e3t_b(:,:,jk) 
    451          END DO 
    452          ! III - Barotropic repartition of the sea surface height over the baroclinic profile 
    453          ! ================================================================================== 
    454          ! add ( ssh increment + "baroclinicity error" ) proportionly to e3t(n) 
    455          ! - ML - baroclinicity error should be better treated in the future 
    456          !        i.e. locally and not spread over the water column. 
    457          !        (keep in mind that the idea is to reduce Eulerian velocity as much as possible) 
    458          zht(:,:) = 0. 
    459          DO jk = 1, jpkm1 
    460             zht(:,:)  = zht(:,:) + tilde_e3t_a(:,:,jk) * tmask(:,:,jk) 
    461          END DO 
    462          z_scale(:,:) =  - zht(:,:) / ( ht_0(:,:) + sshn(:,:) + 1. - ssmask(:,:) ) 
    463          DO jk = 1, jpkm1 
    464             dtilde_e3t_a(:,:,jk) = dtilde_e3t_a(:,:,jk) + e3t_n(:,:,jk) * z_scale(:,:) * tmask(:,:,jk) 
    465          END DO 
    466  
    467       ENDIF 
    468  
    469       IF( ln_vvl_ztilde .OR. ln_vvl_layer )  THEN   ! z_tilde or layer coordinate ! 
    470       !                                           ! ---baroclinic part--------- ! 
    471          DO jk = 1, jpkm1 
    472             e3t_a(:,:,jk) = e3t_a(:,:,jk) + dtilde_e3t_a(:,:,jk) * tmask(:,:,jk) 
    473          END DO 
    474       ENDIF 
    475  
    476       IF( ln_vvl_dbg .AND. .NOT. ll_do_bclinic ) THEN   ! - ML - test: control prints for debuging 
    477          ! 
    478          IF( lwp ) WRITE(numout, *) 'kt =', kt 
    479          IF ( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN 
    480             z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( zht(:,:) ) ) 
    481             IF( lk_mpp ) CALL mpp_max( z_tmax )                             ! max over the global domain 
    482             IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(SUM(tilde_e3t_a))) =', z_tmax 
    483          END IF 
    484          ! 
    485          zht(:,:) = 0.0_wp 
    486          DO jk = 1, jpkm1 
    487             zht(:,:) = zht(:,:) + e3t_n(:,:,jk) * tmask(:,:,jk) 
    488          END DO 
    489          z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + sshn(:,:) - zht(:,:) ) ) 
    490          IF( lk_mpp ) CALL mpp_max( z_tmax )                                ! max over the global domain 
    491          IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+sshn-SUM(e3t_n))) =', z_tmax 
    492          ! 
    493          zht(:,:) = 0.0_wp 
    494          DO jk = 1, jpkm1 
    495             zht(:,:) = zht(:,:) + e3t_a(:,:,jk) * tmask(:,:,jk) 
    496          END DO 
    497          z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + ssha(:,:) - zht(:,:) ) ) 
    498          IF( lk_mpp ) CALL mpp_max( z_tmax )                                ! max over the global domain 
    499          IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+ssha-SUM(e3t_a))) =', z_tmax 
    500          ! 
    501          zht(:,:) = 0.0_wp 
    502          DO jk = 1, jpkm1 
    503             zht(:,:) = zht(:,:) + e3t_b(:,:,jk) * tmask(:,:,jk) 
    504          END DO 
    505          z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + sshb(:,:) - zht(:,:) ) ) 
    506          IF( lk_mpp ) CALL mpp_max( z_tmax )                                ! max over the global domain 
    507          IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+sshb-SUM(e3t_b))) =', z_tmax 
    508          ! 
    509          z_tmax = MAXVAL( tmask(:,:,1) *  ABS( sshb(:,:) ) ) 
    510          IF( lk_mpp ) CALL mpp_max( z_tmax )                                ! max over the global domain 
    511          IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(sshb))) =', z_tmax 
    512          ! 
    513          z_tmax = MAXVAL( tmask(:,:,1) *  ABS( sshn(:,:) ) ) 
    514          IF( lk_mpp ) CALL mpp_max( z_tmax )                                ! max over the global domain 
    515          IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(sshn))) =', z_tmax 
    516          ! 
    517          z_tmax = MAXVAL( tmask(:,:,1) *  ABS( ssha(:,:) ) ) 
    518          IF( lk_mpp ) CALL mpp_max( z_tmax )                                ! max over the global domain 
    519          IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(ssha))) =', z_tmax 
    520       END IF 
    521  
    522       ! *********************************** ! 
    523       ! After scale factors at u- v- points ! 
    524       ! *********************************** ! 
    525  
    526       CALL dom_vvl_interpol( e3t_a(:,:,:), e3u_a(:,:,:), 'U' ) 
    527       CALL dom_vvl_interpol( e3t_a(:,:,:), e3v_a(:,:,:), 'V' ) 
    528  
    529       ! *********************************** ! 
    530       ! After depths at u- v points         ! 
    531       ! *********************************** ! 
    532  
    533       hu_a(:,:) = e3u_a(:,:,1) * umask(:,:,1) 
    534       hv_a(:,:) = e3v_a(:,:,1) * vmask(:,:,1) 
    535       DO jk = 2, jpkm1 
    536          hu_a(:,:) = hu_a(:,:) + e3u_a(:,:,jk) * umask(:,:,jk) 
    537          hv_a(:,:) = hv_a(:,:) + e3v_a(:,:,jk) * vmask(:,:,jk) 
    538       END DO 
    539       !                                        ! Inverse of the local depth 
    540 !!gm BUG ?  don't understand the use of umask_i here ..... 
    541       r1_hu_a(:,:) = ssumask(:,:) / ( hu_a(:,:) + 1._wp - ssumask(:,:) ) 
    542       r1_hv_a(:,:) = ssvmask(:,:) / ( hv_a(:,:) + 1._wp - ssvmask(:,:) ) 
    543       ! 
    544       CALL wrk_dealloc( jpi,jpj,       zht, z_scale, zwu, zwv, zhdiv ) 
    545       CALL wrk_dealloc( jpi,jpj,jpk,   ze3t ) 
    546       ! 
    547       IF( nn_timing == 1 )  CALL timing_stop('dom_vvl_sf_nxt') 
    548       ! 
    549    END SUBROUTINE dom_vvl_sf_nxt 
    550  
    551  
    552    SUBROUTINE dom_vvl_sf_swp( kt ) 
    553       !!---------------------------------------------------------------------- 
    554       !!                ***  ROUTINE dom_vvl_sf_swp  *** 
    555       !!                    
    556       !! ** Purpose :  compute time filter and swap of scale factors  
    557       !!               compute all depths and related variables for next time step 
    558       !!               write outputs and restart file 
    559       !! 
    560       !! ** Method  :  - swap of e3t with trick for volume/tracer conservation 
    561       !!               - reconstruct scale factor at other grid points (interpolate) 
    562       !!               - recompute depths and water height fields 
    563       !! 
    564       !! ** Action  :  - e3t_(b/n), tilde_e3t_(b/n) and e3(u/v)_n ready for next time step 
    565       !!               - Recompute: 
    566       !!                    e3(u/v)_b        
    567       !!                    e3w_n            
    568       !!                    e3(u/v)w_b       
    569       !!                    e3(u/v)w_n       
    570       !!                    gdept_n, gdepw_n  and gde3w_n 
    571       !!                    h(u/v) and h(u/v)r 
    572       !! 
    573       !! Reference  : Leclair, M., and G. Madec, 2009, Ocean Modelling. 
    574       !!              Leclair, M., and G. Madec, 2011, Ocean Modelling. 
    575       !!---------------------------------------------------------------------- 
    576       INTEGER, INTENT( in ) ::   kt   ! time step 
    577       ! 
    578       INTEGER  ::   ji, jj, jk   ! dummy loop indices 
    579       REAL(wp) ::   zcoef        ! local scalar 
    580       !!---------------------------------------------------------------------- 
    581       ! 
    582       IF( ln_linssh )   RETURN      ! No calculation in linear free surface 
    583       ! 
    584       IF( nn_timing == 1 )  CALL timing_start('dom_vvl_sf_swp') 
    585       ! 
    586       IF( kt == nit000 )   THEN 
    587          IF(lwp) WRITE(numout,*) 
    588          IF(lwp) WRITE(numout,*) 'dom_vvl_sf_swp : - time filter and swap of scale factors' 
    589          IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~   - interpolate scale factors and compute depths for next time step' 
    590       ENDIF 
    591       ! 
    592       ! Time filter and swap of scale factors 
    593       ! ===================================== 
    594       ! - ML - e3(t/u/v)_b are allready computed in dynnxt. 
    595       IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN 
    596          IF( neuler == 0 .AND. kt == nit000 ) THEN 
    597             tilde_e3t_b(:,:,:) = tilde_e3t_n(:,:,:) 
    598          ELSE 
    599             tilde_e3t_b(:,:,:) = tilde_e3t_n(:,:,:) &  
    600             &         + atfp * ( tilde_e3t_b(:,:,:) - 2.0_wp * tilde_e3t_n(:,:,:) + tilde_e3t_a(:,:,:) ) 
    601          ENDIF 
    602          tilde_e3t_n(:,:,:) = tilde_e3t_a(:,:,:) 
    603       ENDIF 
    604       gdept_b(:,:,:) = gdept_n(:,:,:) 
    605       gdepw_b(:,:,:) = gdepw_n(:,:,:) 
    606  
    607       e3t_n(:,:,:) = e3t_a(:,:,:) 
    608       e3u_n(:,:,:) = e3u_a(:,:,:) 
    609       e3v_n(:,:,:) = e3v_a(:,:,:) 
    610  
    611       ! Compute all missing vertical scale factor and depths 
    612       ! ==================================================== 
    613       ! Horizontal scale factor interpolations 
    614       ! -------------------------------------- 
    615       ! - ML - e3u_b and e3v_b are allready computed in dynnxt 
    616       ! - JC - hu_b, hv_b, hur_b, hvr_b also 
    617        
    618       CALL dom_vvl_interpol( e3u_n(:,:,:), e3f_n(:,:,:), 'F'  ) 
    619        
    620       ! Vertical scale factor interpolations 
    621       CALL dom_vvl_interpol( e3t_n(:,:,:),  e3w_n(:,:,:), 'W'  ) 
    622       CALL dom_vvl_interpol( e3u_n(:,:,:), e3uw_n(:,:,:), 'UW' ) 
    623       CALL dom_vvl_interpol( e3v_n(:,:,:), e3vw_n(:,:,:), 'VW' ) 
    624       CALL dom_vvl_interpol( e3t_b(:,:,:),  e3w_b(:,:,:), 'W'  ) 
    625       CALL dom_vvl_interpol( e3u_b(:,:,:), e3uw_b(:,:,:), 'UW' ) 
    626       CALL dom_vvl_interpol( e3v_b(:,:,:), e3vw_b(:,:,:), 'VW' ) 
    627  
    628       ! t- and w- points depth (set the isf depth as it is in the initial step) 
    629       gdept_n(:,:,1) = 0.5_wp * e3w_n(:,:,1) 
    630       gdepw_n(:,:,1) = 0.0_wp 
    631       gde3w_n(:,:,1) = gdept_n(:,:,1) - sshn(:,:) 
    632       DO jk = 2, jpk 
    633          DO jj = 1,jpj 
    634             DO ji = 1,jpi 
    635               !    zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk))   ! 0 everywhere tmask = wmask, ie everywhere expect at jk = mikt 
    636                                                                  ! 1 for jk = mikt 
    637                zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk)) 
    638                gdepw_n(ji,jj,jk) = gdepw_n(ji,jj,jk-1) + e3t_n(ji,jj,jk-1) 
    639                gdept_n(ji,jj,jk) =    zcoef  * ( gdepw_n(ji,jj,jk  ) + 0.5 * e3w_n(ji,jj,jk) )  & 
    640                    &             + (1-zcoef) * ( gdept_n(ji,jj,jk-1) +       e3w_n(ji,jj,jk) )  
    641                gde3w_n(ji,jj,jk) = gdept_n(ji,jj,jk) - sshn(ji,jj) 
    642             END DO 
    643          END DO 
    644       END DO 
    645  
    646       ! Local depth and Inverse of the local depth of the water 
    647       ! ------------------------------------------------------- 
    648       hu_n(:,:) = hu_a(:,:)   ;   r1_hu_n(:,:) = r1_hu_a(:,:) 
    649       hv_n(:,:) = hv_a(:,:)   ;   r1_hv_n(:,:) = r1_hv_a(:,:) 
    650       ! 
    651       ht_n(:,:) = e3t_n(:,:,1) * tmask(:,:,1) 
    652       DO jk = 2, jpkm1 
    653          ht_n(:,:) = ht_n(:,:) + e3t_n(:,:,jk) * tmask(:,:,jk) 
    654       END DO 
    655  
    656       ! write restart file 
    657       ! ================== 
    658       IF( lrst_oce )   CALL dom_vvl_rst( kt, 'WRITE' ) 
    659       ! 
    660       IF( nn_timing == 1 )   CALL timing_stop('dom_vvl_sf_swp') 
    661       ! 
    662    END SUBROUTINE dom_vvl_sf_swp 
    663  
    664  
    665241   SUBROUTINE dom_vvl_interpol( pe3_in, pe3_out, pout ) 
    666242      !!--------------------------------------------------------------------- 
     
    684260      IF( nn_timing == 1 )   CALL timing_start('dom_vvl_interpol') 
    685261      ! 
    686       IF(ln_wd) THEN 
    687         zlnwd = 1.0_wp 
    688       ELSE 
    689         zlnwd = 0.0_wp 
    690       END IF 
     262      zlnwd = 0.0_wp 
    691263      ! 
    692264      SELECT CASE ( pout )    !==  type of interpolation  ==! 
     
    772344      ! 
    773345   END SUBROUTINE dom_vvl_interpol 
    774  
    775  
    776    SUBROUTINE dom_vvl_rst( kt, cdrw ) 
    777       !!--------------------------------------------------------------------- 
    778       !!                   ***  ROUTINE dom_vvl_rst  *** 
    779       !!                      
    780       !! ** Purpose :   Read or write VVL file in restart file 
    781       !! 
    782       !! ** Method  :   use of IOM library 
    783       !!                if the restart does not contain vertical scale factors, 
    784       !!                they are set to the _0 values 
    785       !!                if the restart does not contain vertical scale factors increments (z_tilde), 
    786       !!                they are set to 0. 
    787       !!---------------------------------------------------------------------- 
    788       INTEGER         , INTENT(in) ::   kt     ! ocean time-step 
    789       CHARACTER(len=*), INTENT(in) ::   cdrw   ! "READ"/"WRITE" flag 
    790       ! 
    791       INTEGER ::   ji, jj, jk 
    792       INTEGER ::   id1, id2, id3, id4, id5     ! local integers 
    793       !!---------------------------------------------------------------------- 
    794       ! 
    795       IF( nn_timing == 1 )  CALL timing_start('dom_vvl_rst') 
    796       IF( TRIM(cdrw) == 'READ' ) THEN        ! Read/initialise  
    797          !                                   ! =============== 
    798          IF( ln_rstart ) THEN                   !* Read the restart file 
    799             CALL rst_read_open                  !  open the restart file if necessary 
    800             CALL iom_get( numror, jpdom_autoglo, 'sshn'   , sshn    ) 
    801             ! 
    802             id1 = iom_varid( numror, 'e3t_b', ldstop = .FALSE. ) 
    803             id2 = iom_varid( numror, 'e3t_n', ldstop = .FALSE. ) 
    804             id3 = iom_varid( numror, 'tilde_e3t_b', ldstop = .FALSE. ) 
    805             id4 = iom_varid( numror, 'tilde_e3t_n', ldstop = .FALSE. ) 
    806             id5 = iom_varid( numror, 'hdiv_lf', ldstop = .FALSE. ) 
    807             !                             ! --------- ! 
    808             !                             ! all cases ! 
    809             !                             ! --------- ! 
    810             IF( MIN( id1, id2 ) > 0 ) THEN       ! all required arrays exist 
    811                CALL iom_get( numror, jpdom_autoglo, 'e3t_b', e3t_b(:,:,:) ) 
    812                CALL iom_get( numror, jpdom_autoglo, 'e3t_n', e3t_n(:,:,:) ) 
    813                ! needed to restart if land processor not computed  
    814                IF(lwp) write(numout,*) 'dom_vvl_rst : e3t_b and e3t_n found in restart files' 
    815                WHERE ( tmask(:,:,:) == 0.0_wp )  
    816                   e3t_n(:,:,:) = e3t_0(:,:,:) 
    817                   e3t_b(:,:,:) = e3t_0(:,:,:) 
    818                END WHERE 
    819                IF( neuler == 0 ) THEN 
    820                   e3t_b(:,:,:) = e3t_n(:,:,:) 
    821                ENDIF 
    822             ELSE IF( id1 > 0 ) THEN 
    823                IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t_n not found in restart files' 
    824                IF(lwp) write(numout,*) 'e3t_n set equal to e3t_b.' 
    825                IF(lwp) write(numout,*) 'neuler is forced to 0' 
    826                CALL iom_get( numror, jpdom_autoglo, 'e3t_b', e3t_b(:,:,:) ) 
    827                e3t_n(:,:,:) = e3t_b(:,:,:) 
    828                neuler = 0 
    829             ELSE IF( id2 > 0 ) THEN 
    830                IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t_b not found in restart files' 
    831                IF(lwp) write(numout,*) 'e3t_b set equal to e3t_n.' 
    832                IF(lwp) write(numout,*) 'neuler is forced to 0' 
    833                CALL iom_get( numror, jpdom_autoglo, 'e3t_n', e3t_n(:,:,:) ) 
    834                e3t_b(:,:,:) = e3t_n(:,:,:) 
    835                neuler = 0 
    836             ELSE 
    837                IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t_n not found in restart file' 
    838                IF(lwp) write(numout,*) 'Compute scale factor from sshn' 
    839                IF(lwp) write(numout,*) 'neuler is forced to 0' 
    840                DO jk = 1, jpk 
    841                   e3t_n(:,:,jk) =  e3t_0(:,:,jk) * ( ht_0(:,:) + sshn(:,:) ) & 
    842                       &                          / ( ht_0(:,:) + 1._wp - ssmask(:,:) ) * tmask(:,:,jk)   & 
    843                       &          + e3t_0(:,:,jk)                               * (1._wp -tmask(:,:,jk)) 
    844                END DO 
    845                e3t_b(:,:,:) = e3t_n(:,:,:) 
    846                neuler = 0 
    847             ENDIF 
    848             !                             ! ----------- ! 
    849             IF( ln_vvl_zstar ) THEN       ! z_star case ! 
    850                !                          ! ----------- ! 
    851                IF( MIN( id3, id4 ) > 0 ) THEN 
    852                   CALL ctl_stop( 'dom_vvl_rst: z_star cannot restart from a z_tilde or layer run' ) 
    853                ENDIF 
    854                !                          ! ----------------------- ! 
    855             ELSE                          ! z_tilde and layer cases ! 
    856                !                          ! ----------------------- ! 
    857                IF( MIN( id3, id4 ) > 0 ) THEN  ! all required arrays exist 
    858                   CALL iom_get( numror, jpdom_autoglo, 'tilde_e3t_b', tilde_e3t_b(:,:,:) ) 
    859                   CALL iom_get( numror, jpdom_autoglo, 'tilde_e3t_n', tilde_e3t_n(:,:,:) ) 
    860                ELSE                            ! one at least array is missing 
    861                   tilde_e3t_b(:,:,:) = 0.0_wp 
    862                   tilde_e3t_n(:,:,:) = 0.0_wp 
    863                ENDIF 
    864                !                          ! ------------ ! 
    865                IF( ln_vvl_ztilde ) THEN   ! z_tilde case ! 
    866                   !                       ! ------------ ! 
    867                   IF( id5 > 0 ) THEN  ! required array exists 
    868                      CALL iom_get( numror, jpdom_autoglo, 'hdiv_lf', hdiv_lf(:,:,:) ) 
    869                   ELSE                ! array is missing 
    870                      hdiv_lf(:,:,:) = 0.0_wp 
    871                   ENDIF 
    872                ENDIF 
    873             ENDIF 
    874             ! 
    875          ELSE                                   !* Initialize at "rest" 
    876             e3t_b(:,:,:) = e3t_0(:,:,:) 
    877             e3t_n(:,:,:) = e3t_0(:,:,:) 
    878             sshn(:,:) = 0.0_wp 
    879  
    880             IF( ln_wd ) THEN 
    881               DO jj = 1, jpj 
    882                 DO ji = 1, jpi 
    883                   IF( e3t_0(ji,jj,1) <= 0.5_wp * rn_wdmin1 ) THEN 
    884                      e3t_b(ji,jj,:) = 0.5_wp * rn_wdmin1  
    885                      e3t_n(ji,jj,:) = 0.5_wp * rn_wdmin1  
    886                      e3t_a(ji,jj,:) = 0.5_wp * rn_wdmin1  
    887                      sshb(ji,jj) = rn_wdmin1 - bathy(ji,jj) 
    888                      sshn(ji,jj) = rn_wdmin1 - bathy(ji,jj) 
    889                      ssha(ji,jj) = rn_wdmin1 - bathy(ji,jj) 
    890                   ENDIF 
    891                 ENDDO 
    892               ENDDO 
    893             END IF 
    894  
    895             IF( ln_vvl_ztilde .OR. ln_vvl_layer) THEN 
    896                tilde_e3t_b(:,:,:) = 0.0_wp 
    897                tilde_e3t_n(:,:,:) = 0.0_wp 
    898                IF( ln_vvl_ztilde ) hdiv_lf(:,:,:) = 0.0_wp 
    899             END IF 
    900          ENDIF 
    901          ! 
    902       ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN   ! Create restart file 
    903          !                                   ! =================== 
    904          IF(lwp) WRITE(numout,*) '---- dom_vvl_rst ----' 
    905          !                                           ! --------- ! 
    906          !                                           ! all cases ! 
    907          !                                           ! --------- ! 
    908          CALL iom_rstput( kt, nitrst, numrow, 'e3t_b', e3t_b(:,:,:) ) 
    909          CALL iom_rstput( kt, nitrst, numrow, 'e3t_n', e3t_n(:,:,:) ) 
    910          !                                           ! ----------------------- ! 
    911          IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN  ! z_tilde and layer cases ! 
    912             !                                        ! ----------------------- ! 
    913             CALL iom_rstput( kt, nitrst, numrow, 'tilde_e3t_b', tilde_e3t_b(:,:,:) ) 
    914             CALL iom_rstput( kt, nitrst, numrow, 'tilde_e3t_n', tilde_e3t_n(:,:,:) ) 
    915          END IF 
    916          !                                           ! -------------!     
    917          IF( ln_vvl_ztilde ) THEN                    ! z_tilde case ! 
    918             !                                        ! ------------ ! 
    919             CALL iom_rstput( kt, nitrst, numrow, 'hdiv_lf', hdiv_lf(:,:,:) ) 
    920          ENDIF 
    921          ! 
    922       ENDIF 
    923       ! 
    924       IF( nn_timing == 1 )  CALL timing_stop('dom_vvl_rst') 
    925       ! 
    926    END SUBROUTINE dom_vvl_rst 
    927346 
    928347 
     
    990409      ! 
    991410      IF( ioptio /= 1 )   CALL ctl_stop( 'Choose ONE vertical coordinate in namelist nam_vvl' ) 
    992       IF( .NOT. ln_vvl_zstar .AND. ln_isf ) CALL ctl_stop( 'Only vvl_zstar has been tested with ice shelf cavity' ) 
    993411      ! 
    994412      IF(lwp) THEN                   ! Print the choice 
  • branches/2016/dev_r6409_SIMPLIF_2_usrdef_tools/NEMOGCM/NEMO/OPA_SRC/DOM/domwri.F90

    r6689 r6827  
    249249      CALL iom_rstput( 0, 0, inum3, 'e2f', e2f, ktype = jp_r8 ) 
    250250       
    251       CALL iom_rstput( 0, 0, inum3, 'ff', ff, ktype = jp_r8 )           !    ! coriolis factor 
     251      CALL iom_rstput( 0, 0, inum3, 'ff_f', ff_f, ktype = jp_r8 )           !    ! coriolis factor 
     252      CALL iom_rstput( 0, 0, inum3, 'ff_t', ff_t, ktype = jp_r8 )           !    ! coriolis factor 
    252253       
    253254      ! note that mbkt is set to 1 over land ==> use surface tmask 
  • branches/2016/dev_r6409_SIMPLIF_2_usrdef_tools/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90

    r6492 r6827  
    3737   USE oce               ! ocean variables 
    3838   USE dom_oce           ! ocean domain 
    39    USE wet_dry           ! wetting and drying 
    4039   USE closea            ! closed seas 
    41    USE c1d               ! 1D vertical configuration 
    4240   ! 
    4341   USE in_out_manager    ! I/O manager 
     
    147145                          CALL zgr_z            ! Reference z-coordinate system (always called) 
    148146                          CALL zgr_bat          ! Bathymetry fields (levels and meters) 
    149       IF( lk_c1d      )   CALL lbc_lnk( bathy , 'T', 1._wp )   ! 1D config.: same bathy value over the 3x3 domain 
    150147      IF( ln_zco      )   CALL zgr_zco          ! z-coordinate 
    151148      IF( ln_zps      )   CALL zgr_zps          ! Partial step z-coordinate 
     
    155152      ! ----------------------------------- 
    156153      IF( lzoom       )   CALL zgr_bat_zoom     ! correct mbathy in case of zoom subdomain 
    157       IF( .NOT.lk_c1d )   CALL zgr_bat_ctl      ! check bathymetry (mbathy) and suppress isolated ocean points 
     154                          CALL zgr_bat_ctl      ! check bathymetry (mbathy) and suppress isolated ocean points 
    158155                          CALL zgr_bot_level    ! deepest ocean level for t-, u- and v-points 
    159156                          CALL zgr_top_level    ! shallowest ocean level for T-, U-, V- points 
    160       ! 
    161       IF( lk_c1d ) THEN                         ! 1D config.: same mbathy value over the 3x3 domain 
    162          ibat = mbathy(2,2) 
    163          mbathy(:,:) = ibat 
    164       END IF 
    165157      ! 
    166158      IF( nprint == 1 .AND. lwp )   THEN 
     
    19501942      bathy(:,:) = MIN( rn_sbot_max, bathy(:,:) ) 
    19511943 
    1952       IF( .NOT.ln_wd ) THEN                       
    19531944         DO jj = 1, jpj 
    19541945            DO ji = 1, jpi 
     
    19561947            END DO 
    19571948         END DO 
    1958       END IF 
    19591949      !                                        ! ============================= 
    19601950      !                                        ! Define the envelop bathymetry   (hbatt) 
     
    19631953      zenv(:,:) = bathy(:,:) 
    19641954      ! 
    1965       IF( .NOT.ln_wd ) THEN     
    19661955      ! set first land point adjacent to a wet cell to sbot_min as this needs to be included in smoothing 
    19671956         DO jj = 1, jpj 
     
    19861975            END DO 
    19871976         END DO 
    1988       END IF 
    19891977 
    19901978      ! apply lateral boundary condition   CAUTION: keep the value when the lbc field is zero 
     
    20802068      IF(lwp) THEN 
    20812069         WRITE(numout,*) 
    2082          IF( .NOT.ln_wd ) THEN 
    20832070           WRITE(numout,*) ' zgr_sco: minimum depth of the envelop topography set to : ', rn_sbot_min 
    2084          ELSE 
    2085            WRITE(numout,*) ' zgr_sco: minimum positive depth of the envelop topography set to : ', rn_sbot_min 
    2086            WRITE(numout,*) ' zgr_sco: minimum negative depth of the envelop topography set to : ', -rn_wdld 
    2087          ENDIF 
    20882071      ENDIF 
    20892072      hbatu(:,:) = rn_sbot_min 
     
    20992082      END DO 
    21002083 
    2101       IF( ln_wd ) THEN               !avoid the zero depth on T- (U-,V-,F-) points 
    2102         DO jj = 1, jpj 
    2103           DO ji = 1, jpi 
    2104             IF(ABS(hbatt(ji,jj)) < rn_wdmin1) & 
    2105               & hbatt(ji,jj) = SIGN(1._wp, hbatt(ji,jj)) * rn_wdmin1 
    2106             IF(ABS(hbatu(ji,jj)) < rn_wdmin1) & 
    2107               & hbatu(ji,jj) = SIGN(1._wp, hbatu(ji,jj)) * rn_wdmin1 
    2108             IF(ABS(hbatv(ji,jj)) < rn_wdmin1) & 
    2109               & hbatv(ji,jj) = SIGN(1._wp, hbatv(ji,jj)) * rn_wdmin1 
    2110             IF(ABS(hbatf(ji,jj)) < rn_wdmin1) & 
    2111               & hbatf(ji,jj) = SIGN(1._wp, hbatf(ji,jj)) * rn_wdmin1 
    2112           END DO 
    2113         END DO 
    2114       END IF 
    21152084      !  
    21162085      ! Apply lateral boundary condition 
     
    21462115 
    21472116!!bug:  key_helsinki a verifer 
    2148       IF( .NOT.ln_wd ) THEN 
    21492117        hift(:,:) = MIN( hift(:,:), hbatt(:,:) ) 
    21502118        hifu(:,:) = MIN( hifu(:,:), hbatu(:,:) ) 
    21512119        hifv(:,:) = MIN( hifv(:,:), hbatv(:,:) ) 
    21522120        hiff(:,:) = MIN( hiff(:,:), hbatf(:,:) ) 
    2153       END IF 
    21542121 
    21552122      IF( nprint == 1 .AND. lwp )   THEN 
     
    21932160      CALL lbc_lnk( e3vw_0, 'V', 1._wp ) 
    21942161      ! 
    2195       IF( .NOT.ln_wd ) THEN 
    21962162        WHERE( e3t_0 (:,:,:) == 0._wp )   e3t_0 (:,:,:) = 1._wp 
    21972163        WHERE( e3u_0 (:,:,:) == 0._wp )   e3u_0 (:,:,:) = 1._wp 
     
    22012167        WHERE( e3uw_0(:,:,:) == 0._wp )   e3uw_0(:,:,:) = 1._wp 
    22022168        WHERE( e3vw_0(:,:,:) == 0._wp )   e3vw_0(:,:,:) = 1._wp 
    2203       END IF 
    22042169 
    22052170#if defined key_agrif 
     
    22342199               IF( scobot(ji,jj) >= gdept_n(ji,jj,jk) )   mbathy(ji,jj) = MAX( 2, jk ) 
    22352200            END DO 
    2236             IF( ln_wd ) THEN 
    2237               IF( scobot(ji,jj) <= -(rn_wdld - rn_wdmin2)) THEN 
    2238                 mbathy(ji,jj) = 0 
    2239               ELSEIF(scobot(ji,jj) <= rn_wdmin1) THEN 
    2240                 mbathy(ji,jj) = 2 
    2241               ENDIF 
    2242             ELSE 
    2243               IF( scobot(ji,jj) == 0._wp )   mbathy(ji,jj) = 0 
    2244             ENDIF 
    22452201         END DO 
    22462202      END DO 
     
    24252381                   & MAX(hbatt(ji,jj), hbatt(ji+1,jj), hbatt(ji,jj+1), hbatt(ji+1,jj+1)) 
    24262382            DO jk = 1, jpk 
    2427                IF( ln_wd .AND. (ztmpu1 < 0._wp.OR.ABS(ztmpu) < rn_wdmin1) ) THEN 
    2428                  z_esigtu3(ji,jj,jk) = 0.5_wp * ( z_esigt3(ji,jj,jk) + z_esigt3(ji+1,jj,jk) ) 
    2429                  z_esigwu3(ji,jj,jk) = 0.5_wp * ( z_esigw3(ji,jj,jk) + z_esigw3(ji+1,jj,jk) ) 
    2430                ELSE 
    24312383                 z_esigtu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk) ) & 
    24322384                        &              / ztmpu 
    24332385                 z_esigwu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigw3(ji+1,jj,jk) ) & 
    24342386                        &              / ztmpu 
    2435                END IF 
    2436  
    2437                IF( ln_wd .AND. (ztmpv1 < 0._wp.OR.ABS(ztmpv) < rn_wdmin1) ) THEN 
    2438                  z_esigtv3(ji,jj,jk) = 0.5_wp * ( z_esigt3(ji,jj,jk) + z_esigt3(ji,jj+1,jk) ) 
    2439                  z_esigwv3(ji,jj,jk) = 0.5_wp * ( z_esigw3(ji,jj,jk) + z_esigw3(ji,jj+1,jk) ) 
    2440                ELSE 
     2387 
    24412388                 z_esigtv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk) ) & 
    24422389                        &              / ztmpv 
    24432390                 z_esigwv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigw3(ji,jj+1,jk) ) & 
    24442391                        &              / ztmpv 
    2445                END IF 
    2446  
    2447                IF( ln_wd .AND. (ztmpf1 < 0._wp.OR.ABS(ztmpf) < rn_wdmin1) ) THEN 
    2448                  z_esigtf3(ji,jj,jk) = 0.25_wp * ( z_esigt3(ji,jj  ,jk) + z_esigt3(ji+1,jj  ,jk)               & 
    2449                         &                        + z_esigt3(ji,jj+1,jk) + z_esigt3(ji+1,jj+1,jk) ) 
    2450                ELSE 
     2392 
    24512393                 z_esigtf3(ji,jj,jk) = ( hbatt(ji  ,jj  )*z_esigt3(ji  ,jj  ,jk)                               & 
    24522394                        &            +   hbatt(ji+1,jj  )*z_esigt3(ji+1,jj  ,jk)                               & 
    24532395                        &            +   hbatt(ji  ,jj+1)*z_esigt3(ji  ,jj+1,jk)                               & 
    24542396                        &            +   hbatt(ji+1,jj+1)*z_esigt3(ji+1,jj+1,jk) ) / ztmpf 
    2455                END IF 
    24562397 
    24572398               ! 
     
    25822523                  & MAX(hbatt(ji,jj), hbatt(ji+1,jj), hbatt(ji,jj+1), hbatt(ji+1,jj+1)) 
    25832524           DO jk = 1, jpk 
    2584               IF( ln_wd .AND. (ztmpu1 < 0._wp.OR.ABS(ztmpu) < rn_wdmin1) ) THEN 
    2585                 z_esigtu3(ji,jj,jk) = 0.5_wp * ( z_esigt3(ji,jj,jk) + z_esigt3(ji+1,jj,jk) ) 
    2586                 z_esigwu3(ji,jj,jk) = 0.5_wp * ( z_esigw3(ji,jj,jk) + z_esigw3(ji+1,jj,jk) ) 
    2587               ELSE 
    25882525                z_esigtu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk) ) & 
    25892526                       &              / ztmpu 
    25902527                z_esigwu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigw3(ji+1,jj,jk) ) & 
    25912528                       &              / ztmpu 
    2592               END IF 
    2593  
    2594               IF( ln_wd .AND. (ztmpv1 < 0._wp.OR.ABS(ztmpv) < rn_wdmin1) ) THEN 
    2595                 z_esigtv3(ji,jj,jk) = 0.5_wp * ( z_esigt3(ji,jj,jk) + z_esigt3(ji,jj+1,jk) ) 
    2596                 z_esigwv3(ji,jj,jk) = 0.5_wp * ( z_esigw3(ji,jj,jk) + z_esigw3(ji,jj+1,jk) ) 
    2597               ELSE 
    25982529                z_esigtv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk) ) & 
    25992530                       &              / ztmpv 
    26002531                z_esigwv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigw3(ji,jj+1,jk) ) & 
    26012532                       &              / ztmpv 
    2602               END IF 
    2603  
    2604               IF( ln_wd .AND. (ztmpf1 < 0._wp.OR.ABS(ztmpf) < rn_wdmin1) ) THEN 
    2605                 z_esigtf3(ji,jj,jk) = 0.25_wp * ( z_esigt3(ji,jj,jk)   + z_esigt3(ji+1,jj,jk)                 & 
    2606                        &                        + z_esigt3(ji,jj+1,jk) + z_esigt3(ji+1,jj+1,jk) ) 
    2607               ELSE 
     2533 
    26082534                z_esigtf3(ji,jj,jk) = ( hbatt(ji  ,jj  )*z_esigt3(ji  ,jj  ,jk)                               & 
    26092535                       &              + hbatt(ji+1,jj  )*z_esigt3(ji+1,jj  ,jk)                               & 
    26102536                       &              + hbatt(ji  ,jj+1)*z_esigt3(ji  ,jj+1,jk)                               & 
    26112537                       &              + hbatt(ji+1,jj+1)*z_esigt3(ji+1,jj+1,jk) ) / ztmpf 
    2612               END IF 
    26132538 
    26142539             ! Code prior to wetting and drying option (for reference) 
  • branches/2016/dev_r6409_SIMPLIF_2_usrdef_tools/NEMOGCM/NEMO/OPA_SRC/IOM/iom.F90

    r6519 r6827  
    2020   !!-------------------------------------------------------------------- 
    2121   USE dom_oce         ! ocean space and time domain 
    22    USE c1d             ! 1D vertical configuration 
    23    USE flo_oce         ! floats module declarations 
    2422   USE lbclnk          ! lateal boundary condition / mpp exchanges 
    2523   USE iom_def         ! iom variables definitions 
    2624   USE iom_nf90        ! NetCDF format with native NetCDF library 
    2725   USE in_out_manager  ! I/O manager 
    28    USE lib_mpp           ! MPP library 
    29 #if defined key_iomput 
    30    USE sbc_oce, ONLY :   nn_fsbc         ! ocean space and time domain 
    31    USE trc_oce, ONLY :   nn_dttrc        !  !: frequency of step on passive tracers 
    32    USE icb_oce, ONLY :   nclasses, class_num       !  !: iceberg classes 
    33 #if defined key_lim3 
    34    USE ice    , ONLY :   jpl 
    35 #elif defined key_lim2 
    36    USE par_ice_2 
    37 #endif 
     26   USE lib_mpp         ! MPP library 
    3827   USE domngb          ! ocean space and time domain 
    3928   USE phycst          ! physical constants 
    40    USE dianam          ! build name of file 
     29 !SF  USE dianam          ! build name of file 
    4130   USE xios 
    42 # endif 
    4331   USE ioipsl, ONLY :  ju2ymds    ! for calendar 
    44    USE crs             ! Grid coarsening 
    4532 
    4633   IMPLICIT NONE 
     
    139126      ENDIF 
    140127 
    141       IF( TRIM(cdname) == TRIM(cxios_context)//"_crs" ) THEN   
    142          CALL dom_grid_crs   ! Save the parent grid information  & Switch to coarse grid domain 
    143          ! 
    144          CALL set_grid( "T", glamt_crs, gphit_crs )  
    145          CALL set_grid( "U", glamu_crs, gphiu_crs )  
    146          CALL set_grid( "V", glamv_crs, gphiv_crs )  
    147          CALL set_grid( "W", glamt_crs, gphit_crs )  
    148          CALL set_grid_znl( gphit_crs ) 
    149           ! 
    150          CALL dom_grid_glo   ! Return to parent grid domain 
    151          ! 
    152          IF( ln_cfmeta ) THEN   ! Add additional grid metadata 
    153             CALL iom_set_domain_attr("grid_T", area = e1e2t_crs(nldi:nlei, nldj:nlej)) 
    154             CALL iom_set_domain_attr("grid_U", area = e1u_crs(nldi:nlei, nldj:nlej) * e2u_crs(nldi:nlei, nldj:nlej)) 
    155             CALL iom_set_domain_attr("grid_V", area = e1v_crs(nldi:nlei, nldj:nlej) * e2v_crs(nldi:nlei, nldj:nlej)) 
    156             CALL iom_set_domain_attr("grid_W", area = e1e2t_crs(nldi:nlei, nldj:nlej)) 
    157             CALL set_grid_bounds( "T", glamf_crs, gphif_crs, glamt_crs, gphit_crs ) 
    158             CALL set_grid_bounds( "U", glamv_crs, gphiv_crs, glamu_crs, gphiu_crs ) 
    159             CALL set_grid_bounds( "V", glamu_crs, gphiu_crs, glamv_crs, gphiv_crs ) 
    160             CALL set_grid_bounds( "W", glamf_crs, gphif_crs, glamt_crs, gphit_crs ) 
    161          ENDIF 
    162       ENDIF 
    163  
    164128      ! vertical grid definition 
    165129      CALL iom_set_axis_attr( "deptht", gdept_1d ) 
     
    180144      CALL iom_set_axis_attr( "depthw", bounds=z_bnds ) 
    181145 
    182 # if defined key_floats 
    183       CALL iom_set_axis_attr( "nfloat", (/ (REAL(ji,wp), ji=1,nfloat) /) ) 
    184 # endif 
    185 #if defined key_lim3 || defined key_lim2 
    186       CALL iom_set_axis_attr( "ncatice", (/ (REAL(ji,wp), ji=1,jpl) /) ) 
    187 #endif 
    188       CALL iom_set_axis_attr( "icbcla", class_num ) 
    189146      CALL iom_set_axis_attr( "iax_20C", (/ REAL(20,wp) /) ) 
    190147      CALL iom_set_axis_attr( "iax_28C", (/ REAL(28,wp) /) ) 
     
    882839            ENDIF 
    883840             
    884             ! C1D case : always call lbc_lnk to replicate the central value over the whole 3X3 domain 
    885             IF( lk_c1d .AND. PRESENT(pv_r2d) )   CALL lbc_lnk( pv_r2d,'Z',1. ) 
    886             IF( lk_c1d .AND. PRESENT(pv_r3d) )   CALL lbc_lnk( pv_r3d,'Z',1. ) 
    887      
    888841            !--- Apply scale_factor and offset 
    889842            zscf = iom_file(kiomid)%scf(idvar)      ! scale factor 
     
    14591412      ! frequency of the call of iom_put (attribut: freq_op) 
    14601413      WRITE(cl1,'(i1)')        1   ;   CALL iom_set_field_attr('field_definition', freq_op = cl1//'ts', freq_offset='0ts') 
    1461       WRITE(cl1,'(i1)')  nn_fsbc   ;   CALL iom_set_field_attr('SBC'             , freq_op = cl1//'ts', freq_offset='0ts') 
    1462       WRITE(cl1,'(i1)')  nn_fsbc   ;   CALL iom_set_field_attr('SBC_scalar'      , freq_op = cl1//'ts', freq_offset='0ts') 
    1463       WRITE(cl1,'(i1)') nn_dttrc   ;   CALL iom_set_field_attr('ptrc_T'          , freq_op = cl1//'ts', freq_offset='0ts') 
    1464       WRITE(cl1,'(i1)') nn_dttrc   ;   CALL iom_set_field_attr('diad_T'          , freq_op = cl1//'ts', freq_offset='0ts') 
    14651414        
    14661415      ! output file names (attribut: name) 
  • branches/2016/dev_r6409_SIMPLIF_2_usrdef_tools/NEMOGCM/NEMO/OPA_SRC/TRA/eosbn2.F90

    r6505 r6827  
    3939   USE dom_oce        ! ocean space and time domain 
    4040   USE phycst         ! physical constants 
    41    USE stopar         ! Stochastic T/S fluctuations 
    42    USE stopts         ! Stochastic T/S fluctuations 
    4341   ! 
    4442   USE in_out_manager ! I/O manager 
     
    335333      CASE( np_teos10, np_eos80 )                !==  polynomial TEOS-10 / EOS-80 ==! 
    336334         ! 
    337          ! Stochastic equation of state 
    338          IF ( ln_sto_eos ) THEN 
    339             ALLOCATE(zn0_sto(1:2*nn_sto_eos)) 
    340             ALLOCATE(zn_sto(1:2*nn_sto_eos)) 
    341             ALLOCATE(zsign(1:2*nn_sto_eos)) 
    342             DO jsmp = 1, 2*nn_sto_eos, 2 
    343               zsign(jsmp)   = 1._wp 
    344               zsign(jsmp+1) = -1._wp 
    345             END DO 
    346             ! 
    347             DO jk = 1, jpkm1 
    348                DO jj = 1, jpj 
    349                   DO ji = 1, jpi 
    350                      ! 
    351                      ! compute density (2*nn_sto_eos) times: 
    352                      ! (1) for t+dt, s+ds (with the random TS fluctutation computed in sto_pts) 
    353                      ! (2) for t-dt, s-ds (with the opposite fluctuation) 
    354                      DO jsmp = 1, nn_sto_eos*2 
    355                         jdof   = (jsmp + 1) / 2 
    356                         zh     = pdep(ji,jj,jk) * r1_Z0                                  ! depth 
    357                         zt     = (pts (ji,jj,jk,jp_tem) + pts_ran(ji,jj,jk,jp_tem,jdof) * zsign(jsmp)) * r1_T0    ! temperature 
    358                         zstemp = pts  (ji,jj,jk,jp_sal) + pts_ran(ji,jj,jk,jp_sal,jdof) * zsign(jsmp) 
    359                         zs     = SQRT( ABS( zstemp + rdeltaS ) * r1_S0 )   ! square root salinity 
    360                         ztm    = tmask(ji,jj,jk)                                         ! tmask 
    361                         ! 
    362                         zn3 = EOS013*zt   & 
    363                            &   + EOS103*zs+EOS003 
    364                            ! 
    365                         zn2 = (EOS022*zt   & 
    366                            &   + EOS112*zs+EOS012)*zt   & 
    367                            &   + (EOS202*zs+EOS102)*zs+EOS002 
    368                            ! 
    369                         zn1 = (((EOS041*zt   & 
    370                            &   + EOS131*zs+EOS031)*zt   & 
    371                            &   + (EOS221*zs+EOS121)*zs+EOS021)*zt   & 
    372                            &   + ((EOS311*zs+EOS211)*zs+EOS111)*zs+EOS011)*zt   & 
    373                            &   + (((EOS401*zs+EOS301)*zs+EOS201)*zs+EOS101)*zs+EOS001 
    374                            ! 
    375                         zn0_sto(jsmp) = (((((EOS060*zt   & 
    376                            &   + EOS150*zs+EOS050)*zt   & 
    377                            &   + (EOS240*zs+EOS140)*zs+EOS040)*zt   & 
    378                            &   + ((EOS330*zs+EOS230)*zs+EOS130)*zs+EOS030)*zt   & 
    379                            &   + (((EOS420*zs+EOS320)*zs+EOS220)*zs+EOS120)*zs+EOS020)*zt   & 
    380                            &   + ((((EOS510*zs+EOS410)*zs+EOS310)*zs+EOS210)*zs+EOS110)*zs+EOS010)*zt   & 
    381                            &   + (((((EOS600*zs+EOS500)*zs+EOS400)*zs+EOS300)*zs+EOS200)*zs+EOS100)*zs+EOS000 
    382                            ! 
    383                         zn_sto(jsmp)  = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0_sto(jsmp) 
    384                      END DO 
    385                      ! 
    386                      ! compute stochastic density as the mean of the (2*nn_sto_eos) densities 
    387                      prhop(ji,jj,jk) = 0._wp ; prd(ji,jj,jk) = 0._wp 
    388                      DO jsmp = 1, nn_sto_eos*2 
    389                         prhop(ji,jj,jk) = prhop(ji,jj,jk) + zn0_sto(jsmp)                      ! potential density referenced at the surface 
    390                         ! 
    391                         prd(ji,jj,jk) = prd(ji,jj,jk) + (  zn_sto(jsmp) * r1_rau0 - 1._wp  )   ! density anomaly (masked) 
    392                      END DO 
    393                      prhop(ji,jj,jk) = 0.5_wp * prhop(ji,jj,jk) * ztm / nn_sto_eos 
    394                      prd  (ji,jj,jk) = 0.5_wp * prd  (ji,jj,jk) * ztm / nn_sto_eos 
    395                   END DO 
    396                END DO 
    397             END DO 
    398             DEALLOCATE(zn0_sto,zn_sto,zsign) 
    399          ! Non-stochastic equation of state 
    400          ELSE 
    401335            DO jk = 1, jpkm1 
    402336               DO jj = 1, jpj 
     
    437371               END DO 
    438372            END DO 
    439          ENDIF 
    440373          
    441374      CASE( np_seos )                !==  simplified EOS  ==! 
  • branches/2016/dev_r6409_SIMPLIF_2_usrdef_tools/NEMOGCM/NEMO/OPA_SRC/nemogcm.F90

    r6152 r6827  
    4848   USE mppini         ! shared/distributed memory setting (mpp_init routine) 
    4949   USE domain         ! domain initialization             (dom_init routine) 
    50 #if defined key_nemocice_decomp 
    51    USE ice_domain_size, only: nx_global, ny_global 
    52 #endif 
    53    USE tideini        ! tidal components initialization   (tide_ini routine) 
    54    USE bdyini         ! open boundary cond. setting       (bdy_init routine) 
    55    USE bdydta         ! open boundary cond. setting   (bdy_dta_init routine) 
    56    USE bdytides       ! open boundary cond. setting   (bdytide_init routine) 
    57    USE sbctide, ONLY  : lk_tide 
    58    USE istate         ! initial state setting          (istate_init routine) 
    59    USE ldfdyn         ! lateral viscosity setting      (ldfdyn_init routine) 
    60    USE ldftra         ! lateral diffusivity setting    (ldftra_init routine) 
    61    USE zdfini         ! vertical physics setting          (zdf_init routine) 
    6250   USE phycst         ! physical constant                  (par_cst routine) 
    63    USE trdini         ! dyn/tra trends initialization     (trd_init routine) 
    64    USE asminc         ! assimilation increments      
    65    USE asmbkg         ! writing out state trajectory 
    66    USE diaptr         ! poleward transports           (dia_ptr_init routine) 
    67    USE diadct         ! sections transports           (dia_dct_init routine) 
    68    USE diaobs         ! Observation diagnostics       (dia_obs_init routine) 
    69    USE diacfl         ! CFL diagnostics               (dia_cfl_init routine) 
    7051   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 
    71    USE step           ! NEMO time-stepping                 (stp     routine) 
    72    USE icbini         ! handle bergs, initialisation 
    73    USE icbstp         ! handle bergs, calving, themodynamics and transport 
    74    USE cpl_oasis3     ! OASIS3 coupling 
    75    USE c1d            ! 1D configuration 
    76    USE step_c1d       ! Time stepping loop for the 1D configuration 
    77    USE dyndmp         ! Momentum damping 
    78    USE stopar         ! Stochastic param.: ??? 
    79    USE stopts         ! Stochastic param.: ??? 
    80 #if defined key_top 
    81    USE trcini         ! passive tracer initialisation 
    82 #endif 
    8352   USE lib_mpp        ! distributed memory computing 
    84    USE diurnal_bulk    ! diurnal bulk SST  
    85    USE step_diu        ! diurnal bulk SST timestepping (called from here if run offline) 
    8653#if defined key_iomput 
    8754   USE xios           ! xIOserver 
    8855#endif 
    89    USE crsini         ! initialise grid coarsening utility 
    9056   USE lbcnfd , ONLY  : isendto, nsndto, nfsloop, nfeloop   ! Setup of north fold exchanges  
    91    USE sbc_oce, ONLY  : lk_oasis 
    92    USE diatmb          ! Top,middle,bottom output 
    93    USE dia25h          ! 25h mean output 
    94    USE wet_dry         ! Wetting and drying setting   (wad_init routine) 
    9557 
    9658   IMPLICIT NONE 
     
    13799      CALL Agrif_Declare_Var_dom   ! AGRIF: set the meshes for DOM 
    138100      CALL Agrif_Declare_Var       !  "      "   "   "      "  DYN/TRA  
    139 # if defined key_top 
    140       CALL Agrif_Declare_Var_top   !  "      "   "   "      "  TOP 
    141 # endif 
    142 # if defined key_lim2 
    143       CALL Agrif_Declare_Var_lim2  !  "      "   "   "      "  LIM 
    144 # endif 
    145101#endif 
    146102      ! check that all process are still there... If some process have an error, 
     
    151107 
    152108      !                            !-----------------------! 
    153       !                            !==   time stepping   ==! 
    154       !                            !-----------------------! 
    155       istp = nit000 
    156 #if defined key_c1d 
    157          DO WHILE ( istp <= nitend .AND. nstop == 0 ) 
    158             CALL stp_c1d( istp ) 
    159             istp = istp + 1 
    160          END DO 
    161 #else 
    162           IF( lk_asminc ) THEN 
    163              IF( ln_bkgwri ) CALL asm_bkg_wri( nit000 - 1 )    ! Output background fields 
    164              IF( ln_asmdin ) THEN                        ! Direct initialization 
    165                 IF( ln_trainc ) CALL tra_asm_inc( nit000 - 1 )    ! Tracers 
    166                 IF( ln_dyninc ) CALL dyn_asm_inc( nit000 - 1 )    ! Dynamics 
    167                 IF( ln_sshinc ) CALL ssh_asm_inc( nit000 - 1 )    ! SSH 
    168              ENDIF 
    169           ENDIF 
    170  
    171 #if defined key_agrif 
    172           CALL Agrif_Regrid() 
    173 #endif 
    174  
    175          DO WHILE ( istp <= nitend .AND. nstop == 0 ) 
    176 #if defined key_agrif 
    177             CALL stp                         ! AGRIF: time stepping 
    178 #else 
    179             IF ( .NOT. ln_diurnal_only ) THEN  
    180                CALL stp( istp )                 ! standard time stepping  
    181             ELSE  
    182                CALL stp_diurnal( istp )        ! time step only the diurnal SST  
    183             ENDIF  
    184 #endif 
    185             istp = istp + 1 
    186             IF( lk_mpp )   CALL mpp_max( nstop ) 
    187          END DO 
    188 #endif 
    189  
    190       IF( ln_diaobs   )   CALL dia_obs_wri 
    191       ! 
    192       IF( ln_icebergs )   CALL icb_end( nitend ) 
    193  
    194       !                            !------------------------! 
    195109      !                            !==  finalize the run  ==! 
    196110      !                            !------------------------! 
    197       IF(lwp) WRITE(numout,cform_aaa)   ! Flag AAAAAAA 
    198111      ! 
    199112      IF( nstop /= 0 .AND. lwp ) THEN   ! error print 
     
    204117#if defined key_agrif 
    205118      IF( .NOT. Agrif_Root() ) THEN 
    206                          CALL Agrif_ParentGrid_To_ChildGrid() 
    207          IF( ln_diaobs ) CALL dia_obs_wri 
     119                                CALL Agrif_ParentGrid_To_ChildGrid() 
    208120         IF( nn_timing == 1 )   CALL timing_finalize 
    209121                                CALL Agrif_ChildGrid_To_ParentGrid() 
     
    216128#if defined key_iomput 
    217129      CALL xios_finalize                  ! end mpp communications with xios 
    218       IF( lk_oasis )   CALL cpl_finalize  ! end coupling and mpp communications with OASIS 
    219130#else 
    220       IF( lk_oasis ) THEN  
    221          CALL cpl_finalize              ! end coupling and mpp communications with OASIS 
    222       ELSE 
    223          IF( lk_mpp )   CALL mppstop    ! end mpp communications 
    224       ENDIF 
     131      IF( lk_mpp )   CALL mppstop    ! end mpp communications 
    225132#endif 
    226133      ! 
     
    294201#if defined key_iomput 
    295202      IF( Agrif_Root() ) THEN 
    296          IF( lk_oasis ) THEN 
    297             CALL cpl_init( "oceanx", ilocal_comm )                     ! nemo local communicator given by oasis 
    298             CALL xios_initialize( "not used",local_comm=ilocal_comm )    ! send nemo communicator to xios 
    299          ELSE 
    300203            CALL  xios_initialize( "for_xios_mpi_id",return_comm=ilocal_comm )    ! nemo local communicator given by xios 
    301          ENDIF 
    302204      ENDIF 
    303205      ! Nodes selection (control print return in cltxt) 
    304206      narea = mynode( cltxt, 'output.namelist.dyn', numnam_ref, numnam_cfg, numond , nstop, ilocal_comm ) 
    305207#else 
    306       IF( lk_oasis ) THEN 
    307          IF( Agrif_Root() ) THEN 
    308             CALL cpl_init( "oceanx", ilocal_comm )                      ! nemo local communicator given by oasis 
    309          ENDIF 
    310          ! Nodes selection (control print return in cltxt) 
    311          narea = mynode( cltxt, 'output.namelist.dyn', numnam_ref, numnam_cfg, numond , nstop, ilocal_comm ) 
    312       ELSE 
    313208         ilocal_comm = 0 
    314209         ! Nodes selection (control print return in cltxt) 
    315210         narea = mynode( cltxt, 'output.namelist.dyn', numnam_ref, numnam_cfg, numond , nstop ) 
    316       ENDIF 
    317211#endif 
    318212      narea = narea + 1                                     ! mynode return the rank of proc (0 --> jpnij -1 ) 
     
    402296                            CALL     phy_cst    ! Physical constants 
    403297                            CALL     eos_init   ! Equation of state 
    404       IF( lk_c1d        )   CALL     c1d_init   ! 1D column configuration 
    405                             CALL     wad_init   ! Wetting and drying options 
    406298                            CALL     dom_cfg    ! Domain configuration 
    407299                            CALL     dom_init   ! Domain 
    408       IF( ln_crs        )   CALL     crs_init   ! coarsened grid: domain initialization  
    409       IF( ln_nnogather )    CALL nemo_northcomms! northfold neighbour lists (must be done after the masks are defined) 
    410300      IF( ln_ctl        )   CALL prt_ctl_init   ! Print control 
    411        
    412       CALL diurnal_sst_bulk_init            ! diurnal sst 
    413       IF ( ln_diurnal ) CALL diurnal_sst_coolskin_init   ! cool skin    
    414        
    415       ! IF ln_diurnal_only, then we only want a subset of the initialisation routines 
    416       IF ( ln_diurnal_only ) THEN 
    417          CALL  istate_init   ! ocean initial state (Dynamics and tracers) 
    418          CALL     sbc_init   ! Forcings : surface module 
    419          CALL tra_qsr_init   ! penetrative solar radiation qsr 
    420          IF( ln_diaobs     ) THEN                  ! Observation & model comparison 
    421             CALL dia_obs_init            ! Initialize observational data 
    422             CALL dia_obs( nit000 - 1 )   ! Observation operator for restart 
    423          ENDIF      
    424          !                                     ! Assimilation increments 
    425          IF( lk_asminc     )   CALL asm_inc_init   ! Initialize assimilation increments 
    426                   
    427          IF(lwp) WRITE(numout,*) 'Euler time step switch is ', neuler 
    428          RETURN 
    429       ENDIF 
    430        
    431                             CALL  istate_init   ! ocean initial state (Dynamics and tracers) 
    432  
    433       !                                      ! external forcing  
    434 !!gm to be added : creation and call of sbc_apr_init 
    435       IF( lk_tide       )   CALL    tide_init   ! tidal harmonics 
    436                             CALL     sbc_init   ! surface boundary conditions (including sea-ice) 
    437 !!gm ==>> bdy_init should call bdy_dta_init and bdytide_init  NOT in nemogcm !!! 
    438       IF( lk_bdy        )   CALL     bdy_init   ! Open boundaries initialisation 
    439       IF( lk_bdy        )   CALL bdy_dta_init   ! Open boundaries initialisation of external data arrays 
    440       IF( lk_bdy .AND. lk_tide )   & 
    441          &                  CALL bdytide_init   ! Open boundaries initialisation of tidal harmonic forcing 
    442           
    443       !                                      ! Ocean physics 
    444       !                                         ! Vertical physics 
    445                             CALL     zdf_init      ! namelist read 
    446                             CALL zdf_bfr_init      ! bottom friction 
    447       IF( lk_zdfric     )   CALL zdf_ric_init      ! Richardson number dependent Kz 
    448       IF( lk_zdftke     )   CALL zdf_tke_init      ! TKE closure scheme 
    449       IF( lk_zdfgls     )   CALL zdf_gls_init      ! GLS closure scheme 
    450       IF( lk_zdftmx     )   CALL zdf_tmx_init      ! tidal vertical mixing 
    451       IF( lk_zdfddm     )   CALL zdf_ddm_init      ! double diffusive mixing 
    452           
    453       !                                         ! Lateral physics 
    454                             CALL ldf_tra_init      ! Lateral ocean tracer physics 
    455                             CALL ldf_eiv_init      ! eddy induced velocity param. 
    456                             CALL ldf_dyn_init      ! Lateral ocean momentum physics 
    457  
    458       !                                         ! Active tracers 
    459                             CALL tra_qsr_init      ! penetrative solar radiation qsr 
    460                             CALL tra_bbc_init      ! bottom heat flux 
    461       IF( lk_trabbl     )   CALL tra_bbl_init      ! advective (and/or diffusive) bottom boundary layer scheme 
    462                             CALL tra_dmp_init      ! internal tracer damping 
    463                             CALL tra_adv_init      ! horizontal & vertical advection 
    464                             CALL tra_ldf_init      ! lateral mixing 
    465                             CALL tra_zdf_init      ! vertical mixing and after tracer fields 
    466  
    467       !                                         ! Dynamics 
    468       IF( lk_c1d        )   CALL dyn_dmp_init      ! internal momentum damping 
    469                             CALL dyn_adv_init      ! advection (vector or flux form) 
    470                             CALL dyn_vor_init      ! vorticity term including Coriolis 
    471                             CALL dyn_ldf_init      ! lateral mixing 
    472                             CALL dyn_hpg_init      ! horizontal gradient of Hydrostatic pressure 
    473                             CALL dyn_zdf_init      ! vertical diffusion 
    474                             CALL dyn_spg_init      ! surface pressure gradient 
    475  
    476 #if defined key_top 
    477       !                                      ! Passive tracers 
    478                             CALL     trc_init 
    479 #endif 
    480       IF( l_ldfslp      )   CALL ldf_slp_init   ! slope of lateral mixing 
    481  
    482       !                                      ! Icebergs 
    483                             CALL icb_init( rdt, nit000)   ! initialise icebergs instance 
    484  
    485       !                                      ! Misc. options 
    486                             CALL sto_par_init   ! Stochastic parametrization 
    487       IF( ln_sto_eos     )  CALL sto_pts_init   ! RRandom T/S fluctuations 
    488       
    489       !                                      ! Diagnostics 
    490       IF( lk_floats     )   CALL     flo_init   ! drifting Floats 
    491                             CALL dia_cfl_init   ! Initialise CFL diagnostics 
    492       IF( lk_diaar5     )   CALL dia_ar5_init   ! ar5 diag 
    493                             CALL dia_ptr_init   ! Poleward TRansports initialization 
    494       IF( lk_diadct     )   CALL dia_dct_init   ! Sections tranports 
    495                             CALL dia_hsb_init   ! heat content, salt content and volume budgets 
    496                             CALL     trd_init   ! Mixed-layer/Vorticity/Integral constraints trends 
    497                             CALL dia_obs_init            ! Initialize observational data 
    498       IF( ln_diaobs     )   CALL dia_obs( nit000 - 1 )   ! Observation operator for restart 
    499  
    500       !                                         ! Assimilation increments 
    501       IF( lk_asminc     )   CALL asm_inc_init   ! Initialize assimilation increments 
    502       IF(lwp) WRITE(numout,*) 'Euler time step switch is ', neuler 
    503                             CALL dia_tmb_init  ! TMB outputs 
    504                             CALL dia_25h_init  ! 25h mean  outputs 
    505  
    506301      ! 
    507302   END SUBROUTINE nemo_init 
     
    600395      ENDIF 
    601396      ! 
    602       IF( nbench == 1 ) THEN              ! Benchmark 
    603          SELECT CASE ( cp_cfg ) 
    604          CASE ( 'gyre' )   ;   CALL ctl_warn( ' The Benchmark is activated ' ) 
    605          CASE DEFAULT      ;   CALL ctl_stop( ' The Benchmark is based on the GYRE configuration:',   & 
    606             &                                 ' cp_cfg = "gyre" in namelist &namcfg or set nbench = 0' ) 
    607          END SELECT 
    608       ENDIF 
    609       ! 
    610397      IF( 1_wp /= SIGN(1._wp,-0._wp)  )   CALL ctl_stop( 'nemo_ctl: The intrinsec SIGN function follows ',  & 
    611398         &                                               'f2003 standard. '                              ,  & 
     
    653440      !! ** Method  : 
    654441      !!---------------------------------------------------------------------- 
    655       USE diawri    , ONLY: dia_wri_alloc 
    656442      USE dom_oce   , ONLY: dom_oce_alloc 
    657       USE trc_oce   , ONLY: trc_oce_alloc 
    658 #if defined key_diadct  
    659       USE diadct    , ONLY: diadct_alloc  
    660 #endif  
    661 #if defined key_bdy 
    662       USE bdy_oce   , ONLY: bdy_oce_alloc 
    663 #endif 
    664443      ! 
    665444      INTEGER :: ierr 
     
    667446      ! 
    668447      ierr =        oce_alloc       ()          ! ocean 
    669       ierr = ierr + dia_wri_alloc   () 
    670448      ierr = ierr + dom_oce_alloc   ()          ! ocean domain 
    671       ierr = ierr + zdf_oce_alloc   ()          ! ocean vertical physics 
    672       ! 
    673       ierr = ierr + trc_oce_alloc   ()          ! shared TRC / TRA arrays 
    674       ! 
    675 #if defined key_diadct  
    676       ierr = ierr + diadct_alloc    ()          !  
    677 #endif  
    678 #if defined key_bdy 
    679       ierr = ierr + bdy_oce_alloc   ()          ! bdy masks (incl. initialization) 
    680 #endif 
    681449      ! 
    682450      IF( lk_mpp    )   CALL mpp_sum( ierr ) 
  • branches/2016/dev_r6409_SIMPLIF_2_usrdef_tools/NEMOGCM/NEMO/OPA_SRC/step_oce.F90

    r6140 r6827  
    99   USE oce              ! ocean dynamics and tracers variables 
    1010   USE dom_oce          ! ocean space and time domain variables 
    11    USE zdf_oce          ! ocean vertical physics variables 
    1211 
    1312   USE daymod           ! calendar                         (day     routine) 
    1413 
    15    USE sbc_oce          ! surface boundary condition: ocean 
    16    USE sbcmod           ! surface boundary condition       (sbc     routine) 
    17    USE sbcrnf           ! surface boundary condition: runoff variables 
    18    USE sbccpl           ! surface boundary condition: coupled formulation (call send at end of step) 
    19    USE sbcapr           ! surface boundary condition: atmospheric pressure 
    20    USE sbctide          ! Tide initialisation 
    21  
    22    USE traqsr           ! solar radiation penetration      (tra_qsr routine) 
    23    USE trasbc           ! surface boundary condition       (tra_sbc routine) 
    24    USE trabbc           ! bottom boundary condition        (tra_bbc routine) 
    25    USE trabbl           ! bottom boundary layer            (tra_bbl routine) 
    26    USE tradmp           ! internal damping                 (tra_dmp routine) 
    27    USE traadv           ! advection scheme control     (tra_adv_ctl routine) 
    28    USE traldf           ! lateral mixing                   (tra_ldf routine) 
    29    USE trazdf           ! vertical mixing                  (tra_zdf routine) 
    30    USE tranxt           ! time-stepping                    (tra_nxt routine) 
    31    USE tranpc           ! non-penetrative convection       (tra_npc routine) 
    32  
    3314   USE eosbn2           ! equation of state                (eos_bn2 routine) 
    3415 
    35    USE divhor           ! horizontal divergence            (div_hor routine) 
    36    USE dynadv           ! advection                        (dyn_adv routine) 
    37    USE dynbfr           ! Bottom friction terms            (dyn_bfr routine) 
    38    USE dynvor           ! vorticity term                   (dyn_vor routine) 
    39    USE dynhpg           ! hydrostatic pressure grad.       (dyn_hpg routine) 
    40    USE dynldf           ! lateral momentum diffusion       (dyn_ldf routine) 
    41    USE dynzdf           ! vertical diffusion               (dyn_zdf routine) 
    42    USE dynspg           ! surface pressure gradient        (dyn_spg routine) 
    43  
    44    USE dynnxt           ! time-stepping                    (dyn_nxt routine) 
    45  
    46    USE stopar           ! Stochastic parametrization       (sto_par routine) 
    47    USE stopts  
    48  
    49    USE bdy_par          ! for lk_bdy 
    50    USE bdy_oce          ! for dmp logical 
    51    USE bdydta           ! open boundary condition data     (bdy_dta routine) 
    52    USE bdytra           ! bdy cond. for tracers            (bdy_tra routine) 
    53    USE bdydyn3d         ! bdy cond. for baroclinic vel.  (bdy_dyn3d routine) 
    54  
    55    USE sshwzv           ! vertical velocity and ssh        (ssh_nxt routine) 
    56    !                                                       (ssh_swp routine) 
    57    !                                                       (wzv     routine) 
    58    USE domvvl           ! variable vertical scale factors  (dom_vvl_sf_nxt routine) 
    59    !                                                       (dom_vvl_sf_swp routine) 
    60  
    61    USE ldfslp           ! iso-neutral slopes               (ldf_slp routine) 
    62    USE ldfdyn           ! lateral eddy viscosity coef.     (ldf_dyn routine) 
    63    USE ldftra           ! lateral eddy diffusive coef.     (ldf_tra routine) 
    64  
    65    USE zdftmx           ! tide-induced vertical mixing     (zdf_tmx routine) 
    66    USE zdfbfr           ! bottom friction                  (zdf_bfr routine) 
    67    USE zdftke           ! TKE vertical mixing              (zdf_tke routine) 
    68    USE zdfgls           ! GLS vertical mixing              (zdf_gls routine) 
    69    USE zdfddm           ! double diffusion mixing          (zdf_ddm routine) 
    70    USE zdfevd           ! enhanced vertical diffusion      (zdf_evd routine) 
    71    USE zdfric           ! Richardson vertical mixing       (zdf_ric routine) 
    72    USE zdfmxl           ! Mixed-layer depth                (zdf_mxl routine) 
    73  
    74    USE step_diu        ! Time stepping for diurnal sst 
    75    USE diurnal_bulk    ! diurnal SST bulk routines  (diurnal_sst_takaya routine)  
    76    USE cool_skin       ! diurnal cool skin correction (diurnal_sst_coolskin routine)    
    77    USE sbc_oce         ! surface fluxes   
    78     
    79    USE zpshde           ! partial step: hor. derivative     (zps_hde routine) 
    80  
    81    USE diawri           ! Standard run outputs             (dia_wri routine) 
    82    USE diaptr           ! poleward transports              (dia_ptr routine) 
    83    USE diadct           ! sections transports              (dia_dct routine) 
    84    USE diaar5           ! AR5 diagnosics                   (dia_ar5 routine) 
    85    USE diahth           ! thermocline depth                (dia_hth routine) 
    86    USE diafwb           ! freshwater budget                (dia_fwb routine) 
    87    USE diahsb           ! heat, salt and volume budgets    (dia_hsb routine) 
    88    USE diaharm 
    89    USE diacfl 
    90    USE flo_oce          ! floats variables 
    91    USE floats           ! floats computation               (flo_stp routine) 
    92  
    93    USE crsfld           ! Standard output on coarse grid   (crs_fld routine) 
    94  
    95    USE asminc           ! assimilation increments      (tra_asm_inc routine) 
    96    !                                                   (dyn_asm_inc routine) 
    97    USE asmbkg 
    98    USE stpctl           ! time stepping control            (stp_ctl routine) 
    99    USE restart          ! ocean restart                    (rst_wri routine) 
    10016   USE prtctl           ! Print control                    (prt_ctl routine) 
    101  
    102    USE diaobs           ! Observation operator 
    10317 
    10418   USE in_out_manager   ! I/O manager 
     
    11024   USE xios 
    11125#endif 
    112 #if defined key_agrif 
    113    USE agrif_opa_sponge ! Momemtum and tracers sponges 
    114    USE agrif_opa_update ! Update (2-way nesting) 
    115 #endif 
    116 #if defined key_top 
    117    USE trcstp           ! passive tracer time-stepping      (trc_stp routine) 
    118 #endif 
    11926   !!---------------------------------------------------------------------- 
    12027   !! NEMO/OPA 3.7 , NEMO Consortium (2014) 
Note: See TracChangeset for help on using the changeset viewer.