New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
Changeset 7646 for trunk/NEMOGCM/NEMO/OFF_SRC/dtadyn.F90 – NEMO

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

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

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

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

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/NEMOGCM/NEMO/OFF_SRC/dtadyn.F90

    r6140 r7646  
    2222   USE c1d             ! 1D configuration: lk_c1d 
    2323   USE dom_oce         ! ocean domain: variables 
     24   USE domvvl          ! variable volume 
    2425   USE zdf_oce         ! ocean vertical physics: variables 
    2526   USE sbc_oce         ! surface module: variables 
    2627   USE trc_oce         ! share ocean/biogeo variables 
    2728   USE phycst          ! physical constants 
    28    USE ldftra          ! lateral diffusivity coefficients 
    2929   USE trabbl          ! active tracer: bottom boundary layer 
    3030   USE ldfslp          ! lateral diffusion: iso-neutral slopes 
     31   USE sbcrnf          ! river runoffs 
     32   USE ldftra          ! ocean tracer   lateral physics 
    3133   USE zdfmxl          ! vertical physics: mixed layer depth 
    3234   USE eosbn2          ! equation of state - Brunt Vaisala frequency 
    … …  
    3840   USE prtctl          ! print control 
    3941   USE fldread         ! read input fields  
     42   USE wrk_nemo        ! Memory allocation  
    4043   USE timing          ! Timing 
    41    USE wrk_nemo 
     44   USE trc, ONLY : ln_rsttr, numrtr, numrtw, lrst_trc 
    4245 
    4346   IMPLICIT NONE 
    … …  
    4649   PUBLIC   dta_dyn_init   ! called by opa.F90 
    4750   PUBLIC   dta_dyn        ! called by step.F90 
    48  
    49    CHARACTER(len=100) ::   cn_dir       !: Root directory for location of ssr files 
    50    LOGICAL            ::   ln_dynwzv    !: vertical velocity read in a file (T) or computed from u/v (F) 
    51    LOGICAL            ::   ln_dynbbl    !: bbl coef read in a file (T) or computed (F) 
    52    LOGICAL            ::   ln_dynrnf    !: read runoff data in file (T) or set to zero (F) 
    53  
    54    INTEGER  , PARAMETER ::   jpfld = 15     ! maximum number of fields to read 
     51   PUBLIC   dta_dyn_swp   ! called by step.F90 
     52 
     53   CHARACTER(len=100) ::   cn_dir          !: Root directory for location of ssr files 
     54   LOGICAL            ::   ln_dynrnf       !: read runoff data in file (T) or set to zero (F) 
     55   LOGICAL            ::   ln_dynrnf_depth       !: read runoff data in file (T) or set to zero (F) 
     56   REAL(wp)           ::   fwbcorr 
     57 
     58 
     59   INTEGER  , PARAMETER ::   jpfld = 20     ! maximum number of fields to read 
    5560   INTEGER  , SAVE      ::   jf_tem         ! index of temperature 
    5661   INTEGER  , SAVE      ::   jf_sal         ! index of salinity 
    57    INTEGER  , SAVE      ::   jf_uwd         ! index of u-wind 
    58    INTEGER  , SAVE      ::   jf_vwd         ! index of v-wind 
    59    INTEGER  , SAVE      ::   jf_wwd         ! index of w-wind 
     62   INTEGER  , SAVE      ::   jf_uwd         ! index of u-transport 
     63   INTEGER  , SAVE      ::   jf_vwd         ! index of v-transport 
     64   INTEGER  , SAVE      ::   jf_wwd         ! index of v-transport 
    6065   INTEGER  , SAVE      ::   jf_avt         ! index of Kz 
    6166   INTEGER  , SAVE      ::   jf_mld         ! index of mixed layer deptht 
    6267   INTEGER  , SAVE      ::   jf_emp         ! index of water flux 
     68   INTEGER  , SAVE      ::   jf_empb        ! index of water flux 
    6369   INTEGER  , SAVE      ::   jf_qsr         ! index of solar radiation 
    6470   INTEGER  , SAVE      ::   jf_wnd         ! index of wind speed 
    6571   INTEGER  , SAVE      ::   jf_ice         ! index of sea ice cover 
    6672   INTEGER  , SAVE      ::   jf_rnf         ! index of river runoff 
     73   INTEGER  , SAVE      ::   jf_fmf         ! index of downward salt flux 
    6774   INTEGER  , SAVE      ::   jf_ubl         ! index of u-bbl coef 
    6875   INTEGER  , SAVE      ::   jf_vbl         ! index of v-bbl coef 
    69    INTEGER  , SAVE      ::   jf_fmf         ! index of downward salt flux 
    70  
    71    TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: sf_dyn  ! structure of input fields (file informations, fields read) 
     76   INTEGER  , SAVE      ::   jf_div         ! index of e3t 
     77 
     78 
     79   TYPE(FLD), ALLOCATABLE, SAVE, DIMENSION(:) :: sf_dyn  ! structure of input fields (file informations, fields read) 
    7280   !                                               !  
    73    REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: wdta       ! vertical velocity at 2 time step 
    74    REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:,:  ) :: wnow       ! vertical velocity at 2 time step 
    7581   REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: uslpdta    ! zonal isopycnal slopes 
    7682   REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: vslpdta    ! meridional isopycnal slopes 
    7783   REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: wslpidta   ! zonal diapycnal slopes 
    7884   REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: wslpjdta   ! meridional diapycnal slopes 
    79    REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:,:)   :: uslpnow    ! zonal isopycnal slopes 
    80    REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:,:)   :: vslpnow    ! meridional isopycnal slopes 
    81    REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:,:)   :: wslpinow   ! zonal diapycnal slopes 
    82    REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:,:)   :: wslpjnow   ! meridional diapycnal slopes 
    83  
    84    INTEGER :: nrecprev_tem , nrecprev_uwd 
    85  
    86    !! * Substitutions 
    87 #  include "vectopt_loop_substitute.h90" 
     85 
     86   INTEGER, SAVE  :: nprevrec, nsecdyn 
     87 
     88 
    8889   !!---------------------------------------------------------------------- 
    8990   !! NEMO/OFF 3.3 , NEMO Consortium (2010) 
    … …  
    104105      !!             - interpolates data if needed 
    105106      !!---------------------------------------------------------------------- 
    106       USE oce, ONLY:  zts    => tsa 
    107       USE oce, ONLY:  zuslp  => ua   , zvslp  => va 
    108       USE oce, ONLY:  zu     => ub   , zv     => vb,  zw => rke 
    109       ! 
     107      ! 
     108      USE oce, ONLY:  zhdivtr => ua 
    110109      INTEGER, INTENT(in) ::   kt   ! ocean time-step index 
    111       ! 
    112        REAL(wp), DIMENSION(jpi,jpj,jpk     )  :: zwslpi, zwslpj 
    113 !      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts)  :: zts 
    114 !      REAL(wp), DIMENSION(jpi,jpj,jpk     )  :: zuslp, zvslp, zwslpi, zwslpj 
    115 !      REAL(wp), DIMENSION(jpi,jpj,jpk     )  :: zu, zv, zw 
    116       ! 
    117       ! 
    118       INTEGER  ::   ji, jj     ! dummy loop indices 
    119       INTEGER  ::   isecsbc    ! number of seconds between Jan. 1st 00h of nit000 year and the middle of time step 
    120       REAL(wp) ::   ztinta     ! ratio applied to after  records when doing time interpolation 
    121       REAL(wp) ::   ztintb     ! ratio applied to before records when doing time interpolation 
    122       INTEGER  ::   iswap_tem, iswap_uwd    !  
     110      INTEGER             ::   ji, jj, jk 
     111      REAL(wp), POINTER, DIMENSION(:,:)   :: zemp 
     112      ! 
    123113      !!---------------------------------------------------------------------- 
    124114       
    … …  
    126116      IF( nn_timing == 1 )  CALL timing_start( 'dta_dyn') 
    127117      ! 
    128       isecsbc = nsec_year + nsec1jan000  
    129       ! 
    130       IF( kt == nit000 ) THEN 
    131          nrecprev_tem = 0 
    132          nrecprev_uwd = 0 
    133          ! 
    134          CALL fld_read( kt, 1, sf_dyn )      !==   read data at kt time step   ==! 
    135          ! 
    136          IF( l_ldfslp .AND. .NOT.lk_c1d .AND. sf_dyn(jf_tem)%ln_tint ) THEN    ! Computes slopes (here avt is used as workspace)                        
    137             zts(:,:,:,jp_tem) = sf_dyn(jf_tem)%fdta(:,:,:,1) * tmask(:,:,:)   ! temperature 
    138             zts(:,:,:,jp_sal) = sf_dyn(jf_sal)%fdta(:,:,:,1) * tmask(:,:,:)   ! salinity  
    139             avt(:,:,:)        = sf_dyn(jf_avt)%fdta(:,:,:,1) * tmask(:,:,:)   ! vertical diffusive coef. 
    140             CALL dta_dyn_slp( kt, zts, zuslp, zvslp, zwslpi, zwslpj ) 
    141             uslpdta (:,:,:,1) = zuslp (:,:,:)  
    142             vslpdta (:,:,:,1) = zvslp (:,:,:)  
    143             wslpidta(:,:,:,1) = zwslpi(:,:,:)  
    144             wslpjdta(:,:,:,1) = zwslpj(:,:,:)  
    145          ENDIF 
    146          IF( ln_dynwzv .AND. sf_dyn(jf_uwd)%ln_tint )  THEN    ! compute vertical velocity from u/v 
    147             zu(:,:,:) = sf_dyn(jf_uwd)%fdta(:,:,:,1) 
    148             zv(:,:,:) = sf_dyn(jf_vwd)%fdta(:,:,:,1) 
    149             CALL dta_dyn_wzv( zu, zv, zw ) 
    150             wdta(:,:,:,1) = zw(:,:,:) * tmask(:,:,:) 
    151          ENDIF 
    152       ELSE 
    153          nrecprev_tem = sf_dyn(jf_tem)%nrec_a(2) 
    154          nrecprev_uwd = sf_dyn(jf_uwd)%nrec_a(2) 
    155          ! 
    156          CALL fld_read( kt, 1, sf_dyn )      !==   read data at kt time step   ==! 
    157          ! 
    158       ENDIF 
    159       !  
    160       IF( l_ldfslp .AND. .NOT.lk_c1d ) THEN    ! Computes slopes (here avt is used as workspace)                        
    161          iswap_tem = 0 
    162          IF(  kt /= nit000 .AND. ( sf_dyn(jf_tem)%nrec_a(2) - nrecprev_tem ) /= 0 )  iswap_tem = 1 
    163          IF( ( isecsbc > sf_dyn(jf_tem)%nrec_b(2) .AND. iswap_tem == 1 ) .OR. kt == nit000 )  THEN    ! read/update the after data 
    164             IF(lwp) WRITE(numout,*) ' Compute new slopes at kt = ', kt 
    165             IF( sf_dyn(jf_tem)%ln_tint ) THEN                 ! time interpolation of data 
    166                IF( kt /= nit000 ) THEN 
    167                   uslpdta (:,:,:,1) =  uslpdta (:,:,:,2)         ! swap the data 
    168                   vslpdta (:,:,:,1) =  vslpdta (:,:,:,2)   
    169                   wslpidta(:,:,:,1) =  wslpidta(:,:,:,2)  
    170                   wslpjdta(:,:,:,1) =  wslpjdta(:,:,:,2)  
    171                ENDIF 
    172                ! 
    173                zts(:,:,:,jp_tem) = sf_dyn(jf_tem)%fdta(:,:,:,2) * tmask(:,:,:)   ! temperature 
    174                zts(:,:,:,jp_sal) = sf_dyn(jf_sal)%fdta(:,:,:,2) * tmask(:,:,:)   ! salinity  
    175                avt(:,:,:)        = sf_dyn(jf_avt)%fdta(:,:,:,2) * tmask(:,:,:)   ! vertical diffusive coef. 
    176                CALL dta_dyn_slp( kt, zts, zuslp, zvslp, zwslpi, zwslpj ) 
    177                ! 
    178                uslpdta (:,:,:,2) = zuslp (:,:,:)  
    179                vslpdta (:,:,:,2) = zvslp (:,:,:)  
    180                wslpidta(:,:,:,2) = zwslpi(:,:,:)  
    181                wslpjdta(:,:,:,2) = zwslpj(:,:,:)  
    182             ELSE 
    183                zts(:,:,:,jp_tem) = sf_dyn(jf_tem)%fnow(:,:,:) * tmask(:,:,:) 
    184                zts(:,:,:,jp_sal) = sf_dyn(jf_sal)%fnow(:,:,:) * tmask(:,:,:) 
    185                avt(:,:,:)        = sf_dyn(jf_avt)%fnow(:,:,:) * tmask(:,:,:) 
    186                CALL dta_dyn_slp( kt, zts, zuslp, zvslp, zwslpi, zwslpj ) 
    187                uslpnow (:,:,:)   = zuslp (:,:,:)  
    188                vslpnow (:,:,:)   = zvslp (:,:,:)  
    189                wslpinow(:,:,:)   = zwslpi(:,:,:)  
    190                wslpjnow(:,:,:)   = zwslpj(:,:,:)  
    191             ENDIF 
    192          ENDIF 
    193          IF( sf_dyn(jf_tem)%ln_tint )  THEN 
    194             ztinta =  REAL( isecsbc - sf_dyn(jf_tem)%nrec_b(2), wp )  & 
    195                &    / REAL( sf_dyn(jf_tem)%nrec_a(2) - sf_dyn(jf_tem)%nrec_b(2), wp ) 
    196             ztintb =  1. - ztinta 
    197             uslp (:,:,:) = ztintb * uslpdta (:,:,:,1)  + ztinta * uslpdta (:,:,:,2)   
    198             vslp (:,:,:) = ztintb * vslpdta (:,:,:,1)  + ztinta * vslpdta (:,:,:,2)   
    199             wslpi(:,:,:) = ztintb * wslpidta(:,:,:,1)  + ztinta * wslpidta(:,:,:,2)   
    200             wslpj(:,:,:) = ztintb * wslpjdta(:,:,:,1)  + ztinta * wslpjdta(:,:,:,2)   
    201          ELSE 
    202             uslp (:,:,:) = uslpnow (:,:,:) 
    203             vslp (:,:,:) = vslpnow (:,:,:) 
    204             wslpi(:,:,:) = wslpinow(:,:,:) 
    205             wslpj(:,:,:) = wslpjnow(:,:,:) 
    206          ENDIF 
    207       ENDIF 
    208       ! 
    209       IF( ln_dynwzv )  THEN    ! compute vertical velocity from u/v 
    210          iswap_uwd = 0 
    211          IF(  kt /= nit000 .AND. ( sf_dyn(jf_uwd)%nrec_a(2) - nrecprev_uwd ) /= 0 )  iswap_uwd = 1 
    212          IF( ( isecsbc > sf_dyn(jf_uwd)%nrec_b(2) .AND. iswap_uwd == 1 ) .OR. kt == nit000 )  THEN    ! read/update the after data 
    213             IF(lwp) WRITE(numout,*) ' Compute new vertical velocity at kt = ', kt 
    214             IF(lwp) WRITE(numout,*) 
    215             IF( sf_dyn(jf_uwd)%ln_tint ) THEN                 ! time interpolation of data 
    216                IF( kt /= nit000 )  THEN 
    217                   wdta(:,:,:,1) =  wdta(:,:,:,2)     ! swap the data for initialisation 
    218                ENDIF 
    219                zu(:,:,:) = sf_dyn(jf_uwd)%fdta(:,:,:,2) 
    220                zv(:,:,:) = sf_dyn(jf_vwd)%fdta(:,:,:,2) 
    221                CALL dta_dyn_wzv( zu, zv, zw ) 
    222                wdta(:,:,:,2) = zw(:,:,:) * tmask(:,:,:) 
    223             ELSE 
    224                zu(:,:,:) = sf_dyn(jf_uwd)%fnow(:,:,:)  
    225                zv(:,:,:) = sf_dyn(jf_vwd)%fnow(:,:,:) 
    226                CALL dta_dyn_wzv( zu, zv, zw ) 
    227                wnow(:,:,:)  = zw(:,:,:) * tmask(:,:,:) 
    228             ENDIF 
    229          ENDIF 
    230          IF( sf_dyn(jf_uwd)%ln_tint )  THEN 
    231             ztinta =  REAL( isecsbc - sf_dyn(jf_uwd)%nrec_b(2), wp )  & 
    232                &    / REAL( sf_dyn(jf_uwd)%nrec_a(2) - sf_dyn(jf_uwd)%nrec_b(2), wp ) 
    233             ztintb =  1. - ztinta 
    234             wn(:,:,:) = ztintb * wdta(:,:,:,1)  + ztinta * wdta(:,:,:,2)   
    235          ELSE 
    236             wn(:,:,:) = wnow(:,:,:) 
    237          ENDIF 
    238       ENDIF 
     118      nsecdyn = nsec_year + nsec1jan000   ! number of seconds between Jan. 1st 00h of nit000 year and the middle of time step 
     119      ! 
     120      IF( kt == nit000 ) THEN    ;    nprevrec = 0 
     121      ELSE                       ;    nprevrec = sf_dyn(jf_tem)%nrec_a(2) 
     122      ENDIF 
     123      CALL fld_read( kt, 1, sf_dyn )      !=  read data at kt time step   ==! 
     124      ! 
     125      IF( l_ldfslp .AND. .NOT.lk_c1d )   CALL  dta_dyn_slp( kt )    ! Computation of slopes 
    239126      ! 
    240127      tsn(:,:,:,jp_tem) = sf_dyn(jf_tem)%fnow(:,:,:)  * tmask(:,:,:)    ! temperature 
    241128      tsn(:,:,:,jp_sal) = sf_dyn(jf_sal)%fnow(:,:,:)  * tmask(:,:,:)    ! salinity 
    242       ! 
     129      wndm(:,:)         = sf_dyn(jf_wnd)%fnow(:,:,1)  * tmask(:,:,1)    ! wind speed - needed for gas exchange 
     130      fmmflx(:,:)       = sf_dyn(jf_fmf)%fnow(:,:,1)  * tmask(:,:,1)    ! downward salt flux (v3.5+) 
     131      fr_i(:,:)         = sf_dyn(jf_ice)%fnow(:,:,1)  * tmask(:,:,1)    ! Sea-ice fraction 
     132      qsr (:,:)         = sf_dyn(jf_qsr)%fnow(:,:,1)  * tmask(:,:,1)    ! solar radiation 
     133      emp (:,:)         = sf_dyn(jf_emp)%fnow(:,:,1)  * tmask(:,:,1)    ! E-P 
     134      IF( ln_dynrnf ) THEN  
     135         rnf (:,:)      = sf_dyn(jf_rnf)%fnow(:,:,1) * tmask(:,:,1)    ! E-P 
     136         IF( ln_dynrnf_depth .AND. .NOT. ln_linssh )    CALL  dta_dyn_hrnf 
     137      ENDIF 
     138      ! 
     139      un(:,:,:)        = sf_dyn(jf_uwd)%fnow(:,:,:) * umask(:,:,:)    ! effective u-transport 
     140      vn(:,:,:)        = sf_dyn(jf_vwd)%fnow(:,:,:) * vmask(:,:,:)    ! effective v-transport 
     141      wn(:,:,:)        = sf_dyn(jf_wwd)%fnow(:,:,:) * tmask(:,:,:)    ! effective v-transport 
     142      ! 
     143      IF( .NOT.ln_linssh ) THEN 
     144         CALL wrk_alloc(jpi, jpj, zemp ) 
     145         zhdivtr(:,:,:) = sf_dyn(jf_div)%fnow(:,:,:) * tmask(:,:,:)    ! effective u-transport 
     146         emp_b (:,:)    = sf_dyn(jf_empb)%fnow(:,:,1) * tmask(:,:,1)    ! E-P 
     147         zemp   (:,:)   = 0.5_wp * ( emp(:,:) + emp_b(:,:) ) + rnf(:,:) + fwbcorr * tmask(:,:,1) 
     148         CALL dta_dyn_ssh( kt, zhdivtr, sshb, zemp, ssha, e3t_a(:,:,:) )  !=  ssh, vertical scale factor & vertical transport 
     149         CALL wrk_dealloc(jpi, jpj, zemp ) 
     150         !                                           Write in the tracer restart file 
     151         !                                          ******************************* 
     152         IF( lrst_trc ) THEN 
     153            IF(lwp) WRITE(numout,*) 
     154            IF(lwp) WRITE(numout,*) 'dta_dyn_ssh : ssh field written in tracer restart file ',   & 
     155               &                    'at it= ', kt,' date= ', ndastp 
     156            IF(lwp) WRITE(numout,*) '~~~~' 
     157            CALL iom_rstput( kt, nitrst, numrtw, 'sshn', ssha ) 
     158            CALL iom_rstput( kt, nitrst, numrtw, 'sshb', sshn ) 
     159         ENDIF 
     160      ENDIF 
    243161      ! 
    244162      CALL eos    ( tsn, rhd, rhop, gdept_0(:,:,:) ) ! In any case, we need rhop 
    … …  
    247165 
    248166      rn2b(:,:,:) = rn2(:,:,:)         ! need for zdfmxl 
    249       CALL zdf_mxl( kt )                                                   ! In any case, we need mxl  
    250       ! 
    251       avt(:,:,:)       = sf_dyn(jf_avt)%fnow(:,:,:)  * tmask(:,:,:)    ! vertical diffusive coefficient  
    252       un (:,:,:)       = sf_dyn(jf_uwd)%fnow(:,:,:)  * umask(:,:,:)    ! u-velocity 
    253       vn (:,:,:)       = sf_dyn(jf_vwd)%fnow(:,:,:)  * vmask(:,:,:)    ! v-velocity  
    254       IF( .NOT.ln_dynwzv ) &                                          ! w-velocity read in file  
    255          wn (:,:,:)    = sf_dyn(jf_wwd)%fnow(:,:,:) * tmask(:,:,:)     
    256       hmld(:,:)        = sf_dyn(jf_mld)%fnow(:,:,1) * tmask(:,:,1)    ! mixed layer depht 
    257       wndm(:,:)        = sf_dyn(jf_wnd)%fnow(:,:,1) * tmask(:,:,1)    ! wind speed - needed for gas exchange 
    258       emp (:,:)        = sf_dyn(jf_emp)%fnow(:,:,1) * tmask(:,:,1)    ! E-P 
    259       fmmflx(:,:)      = sf_dyn(jf_fmf)%fnow(:,:,1) * tmask(:,:,1)    ! downward salt flux (v3.5+) 
    260       fr_i(:,:)        = sf_dyn(jf_ice)%fnow(:,:,1) * tmask(:,:,1)    ! Sea-ice fraction 
    261       qsr (:,:)        = sf_dyn(jf_qsr)%fnow(:,:,1) * tmask(:,:,1)    ! solar radiation 
    262       IF( ln_dynrnf ) & 
    263       rnf (:,:)        = sf_dyn(jf_rnf)%fnow(:,:,1) * tmask(:,:,1)    ! river runoffs  
    264  
    265       !                                               ! update eddy diffusivity coeff. and/or eiv coeff. at kt 
    266       IF( l_ldftra_time .OR. l_ldfeiv_time )   CALL ldf_tra( kt )  
    267       !                                                      ! bbl diffusive coef 
     167      CALL zdf_mxl( kt )                                                   ! In any case, we need mxl 
     168      ! 
     169      hmld(:,:)         = sf_dyn(jf_mld)%fnow(:,:,1) * tmask(:,:,1)    ! mixed layer depht 
     170      avt(:,:,:)        = sf_dyn(jf_avt)%fnow(:,:,:) * tmask(:,:,:)    ! vertical diffusive coefficient  
     171      ! 
    268172#if defined key_trabbl && ! defined key_c1d 
    269       IF( ln_dynbbl ) THEN                                        ! read in a file 
    270          ahu_bbl(:,:)  = sf_dyn(jf_ubl)%fnow(:,:,1) * umask(:,:,1) 
    271          ahv_bbl(:,:)  = sf_dyn(jf_vbl)%fnow(:,:,1) * vmask(:,:,1) 
    272       ELSE                                                        ! Compute bbl coefficients if needed 
    273          tsb(:,:,:,:) = tsn(:,:,:,:) 
    274          CALL bbl( kt, nit000, 'TRC') 
    275       END IF 
     173      ahu_bbl(:,:)      = sf_dyn(jf_ubl)%fnow(:,:,1) * umask(:,:,1)    ! bbl diffusive coef 
     174      ahv_bbl(:,:)      = sf_dyn(jf_vbl)%fnow(:,:,1) * vmask(:,:,1) 
    276175#endif 
     176      ! 
     177      ! 
     178      CALL eos( tsn, rhd, rhop, gdept_0(:,:,:) ) ! In any case, we need rhop 
    277179      ! 
    278180      IF(ln_ctl) THEN                  ! print control 
    … …  
    283185         CALL prt_ctl(tab3d_1=wn               , clinfo1=' wn      - : ', mask1=tmask, ovlap=1, kdim=jpk   ) 
    284186         CALL prt_ctl(tab3d_1=avt              , clinfo1=' kz      - : ', mask1=tmask, ovlap=1, kdim=jpk   ) 
    285          CALL prt_ctl(tab2d_1=fr_i             , clinfo1=' fr_i    - : ', mask1=tmask, ovlap=1 ) 
    286          CALL prt_ctl(tab2d_1=hmld             , clinfo1=' hmld    - : ', mask1=tmask, ovlap=1 ) 
    287          CALL prt_ctl(tab2d_1=fmmflx           , clinfo1=' fmmflx  - : ', mask1=tmask, ovlap=1 ) 
    288          CALL prt_ctl(tab2d_1=emp              , clinfo1=' emp     - : ', mask1=tmask, ovlap=1 ) 
    289          CALL prt_ctl(tab2d_1=wndm             , clinfo1=' wspd    - : ', mask1=tmask, ovlap=1 ) 
    290          CALL prt_ctl(tab2d_1=qsr              , clinfo1=' qsr     - : ', mask1=tmask, ovlap=1 ) 
     187         CALL prt_ctl(tab3d_1=uslp             , clinfo1=' slp  - u : ', tab3d_2=vslp, clinfo2=' v : ', kdim=jpk) 
     188         CALL prt_ctl(tab3d_1=wslpi            , clinfo1=' slp  - wi: ', tab3d_2=wslpj, clinfo2=' wj: ', kdim=jpk) 
     189!         CALL prt_ctl(tab2d_1=fr_i             , clinfo1=' fr_i    - : ', mask1=tmask, ovlap=1 ) 
     190!         CALL prt_ctl(tab2d_1=hmld             , clinfo1=' hmld    - : ', mask1=tmask, ovlap=1 ) 
     191!         CALL prt_ctl(tab2d_1=fmmflx           , clinfo1=' fmmflx  - : ', mask1=tmask, ovlap=1 ) 
     192!         CALL prt_ctl(tab2d_1=emp              , clinfo1=' emp     - : ', mask1=tmask, ovlap=1 ) 
     193!         CALL prt_ctl(tab2d_1=wndm             , clinfo1=' wspd    - : ', mask1=tmask, ovlap=1 ) 
     194!         CALL prt_ctl(tab2d_1=qsr              , clinfo1=' qsr     - : ', mask1=tmask, ovlap=1 ) 
    291195      ENDIF 
    292196      ! 
    … …  
    310214      INTEGER  :: inum, idv, idimv                   ! local integer 
    311215      INTEGER  :: ios                                ! Local integer output status for namelist read 
    312       !! 
    313       CHARACTER(len=100)            ::  cn_dir   !   Root directory for location of core files 
    314       TYPE(FLD_N), DIMENSION(jpfld) ::  slf_d    ! array of namelist informations on the fields to read 
    315       TYPE(FLD_N) :: sn_tem, sn_sal, sn_mld, sn_emp, sn_ice, sn_qsr, sn_wnd, sn_rnf  ! informations about the fields to be read 
    316       TYPE(FLD_N) :: sn_uwd, sn_vwd, sn_wwd, sn_avt, sn_ubl, sn_vbl, sn_fmf          !   "                                 " 
    317       NAMELIST/namdta_dyn/cn_dir, ln_dynwzv, ln_dynbbl, ln_dynrnf,    & 
    318          &                sn_tem, sn_sal, sn_mld, sn_emp, sn_ice, sn_qsr, sn_wnd, sn_rnf,  & 
    319          &                sn_uwd, sn_vwd, sn_wwd, sn_avt, sn_ubl, sn_vbl, sn_fmf   
    320       !!---------------------------------------------------------------------- 
     216      INTEGER  :: ji, jj, jk 
     217      REAL(wp) :: zcoef 
     218      INTEGER  :: nkrnf_max 
     219      REAL(wp) :: hrnf_max 
     220      !! 
     221      CHARACTER(len=100)            ::  cn_dir        !   Root directory for location of core files 
     222      TYPE(FLD_N), DIMENSION(jpfld) ::  slf_d         ! array of namelist informations on the fields to read 
     223      TYPE(FLD_N) :: sn_uwd, sn_vwd, sn_wwd, sn_empb, sn_emp  ! informations about the fields to be read 
     224      TYPE(FLD_N) :: sn_tem , sn_sal , sn_avt   !   "                 " 
     225      TYPE(FLD_N) :: sn_mld, sn_qsr, sn_wnd , sn_ice , sn_fmf   !   "               " 
     226      TYPE(FLD_N) :: sn_ubl, sn_vbl, sn_rnf    !   "              " 
     227      TYPE(FLD_N) :: sn_div  ! informations about the fields to be read 
     228 
     229      !!---------------------------------------------------------------------- 
     230      ! 
     231      NAMELIST/namdta_dyn/cn_dir, ln_dynrnf, ln_dynrnf_depth,  fwbcorr, & 
     232         &                sn_uwd, sn_vwd, sn_wwd, sn_emp,    & 
     233         &                sn_avt, sn_tem, sn_sal, sn_mld , sn_qsr ,   & 
     234         &                sn_wnd, sn_ice, sn_fmf,                    & 
     235         &                sn_ubl, sn_vbl, sn_rnf,                   & 
     236         &                sn_empb, sn_div  
    321237      ! 
    322238      REWIND( numnam_ref )              ! Namelist namdta_dyn in reference namelist : Offline: init. of dynamical data 
    … …  
    335251         WRITE(numout,*) '~~~~~~~ ' 
    336252         WRITE(numout,*) '   Namelist namdta_dyn' 
    337          WRITE(numout,*) '      vertical velocity read from file (T) or computed (F) ln_dynwzv  = ', ln_dynwzv 
    338          WRITE(numout,*) '      bbl coef read from file (T) or computed (F)          ln_dynbbl  = ', ln_dynbbl 
    339          WRITE(numout,*) '      river runoff option enabled (T) or not (F)           ln_dynrnf  = ', ln_dynrnf 
     253         WRITE(numout,*) '      runoffs option enabled (T) or not (F)            ln_dynrnf        = ', ln_dynrnf 
     254         WRITE(numout,*) '      runoffs is spread in vertical                    ln_dynrnf_depth  = ', ln_dynrnf_depth 
     255         WRITE(numout,*) '      annual global mean of empmr for ssh correction   fwbcorr          = ', fwbcorr 
    340256         WRITE(numout,*) 
    341257      ENDIF 
    342258      !  
    343       IF( ln_dynbbl .AND. ( .NOT.lk_trabbl .OR. lk_c1d ) ) THEN 
    344          CALL ctl_warn( 'dta_dyn_init: bbl option requires key_trabbl activated ; force ln_dynbbl to false' ) 
    345          ln_dynbbl = .FALSE. 
    346       ENDIF 
    347  
    348       jf_tem = 1   ;   jf_sal = 2   ;  jf_mld = 3   ;  jf_emp = 4   ;   jf_fmf  = 5   ;  jf_ice = 6   ;   jf_qsr = 7 
    349       jf_wnd = 8   ;   jf_uwd = 9   ;  jf_vwd = 10  ;  jf_wwd = 11  ;   jf_avt  = 12  ;  jfld  = jf_avt 
    350       ! 
    351       slf_d(jf_tem) = sn_tem   ;   slf_d(jf_sal)  = sn_sal   ;   slf_d(jf_mld) = sn_mld 
    352       slf_d(jf_emp) = sn_emp   ;   slf_d(jf_fmf ) = sn_fmf   ;   slf_d(jf_ice) = sn_ice  
    353       slf_d(jf_qsr) = sn_qsr   ;   slf_d(jf_wnd)  = sn_wnd   ;   slf_d(jf_avt) = sn_avt  
    354       slf_d(jf_uwd) = sn_uwd   ;   slf_d(jf_vwd)  = sn_vwd   ;   slf_d(jf_wwd) = sn_wwd 
    355  
     259 
     260      jf_uwd  = 1     ;   jf_vwd  = 2    ;   jf_wwd = 3    ;   jf_emp = 4    ;   jf_avt = 5 
     261      jf_tem  = 6     ;   jf_sal  = 7    ;   jf_mld = 8    ;   jf_qsr = 9 
     262      jf_wnd  = 10    ;   jf_ice  = 11   ;   jf_fmf = 12   ;   jfld   = jf_fmf 
     263 
     264      ! 
     265      slf_d(jf_uwd)  = sn_uwd    ;   slf_d(jf_vwd)  = sn_vwd   ;   slf_d(jf_wwd) = sn_wwd 
     266      slf_d(jf_emp)  = sn_emp    ;   slf_d(jf_avt)  = sn_avt 
     267      slf_d(jf_tem)  = sn_tem    ;   slf_d(jf_sal)  = sn_sal   ;   slf_d(jf_mld) = sn_mld 
     268      slf_d(jf_qsr)  = sn_qsr    ;   slf_d(jf_wnd)  = sn_wnd   ;   slf_d(jf_ice) = sn_ice 
     269      slf_d(jf_fmf)  = sn_fmf 
     270 
     271      ! 
     272      IF( .NOT.ln_linssh ) THEN 
     273                 jf_div  = jfld + 1    ;         jf_empb  = jfld + 2      ;      jfld = jf_empb 
     274           slf_d(jf_div) = sn_div      ;   slf_d(jf_empb) = sn_empb 
     275      ENDIF 
     276      ! 
     277      IF( lk_trabbl ) THEN 
     278                 jf_ubl  = jfld + 1    ;         jf_vbl  = jfld + 2     ;      jfld = jf_vbl 
     279           slf_d(jf_ubl) = sn_ubl      ;   slf_d(jf_vbl) = sn_vbl 
     280      ENDIF 
    356281      ! 
    357282      IF( ln_dynrnf ) THEN 
    358                 jf_rnf = jfld + 1  ;  jfld  = jf_rnf 
    359          slf_d(jf_rnf) = sn_rnf 
    360          ! Activate runoff key of sbc_oce 
    361          ln_rnf = .true. 
    362          WRITE(numout,*) 'dta_dyn : Activate the runoff data structure from ocean core ( force ln_rnf = .true.) ' 
    363          WRITE(numout,*) 
     283                jf_rnf  = jfld + 1     ;     jfld  = jf_rnf 
     284          slf_d(jf_rnf) = sn_rnf     
    364285      ELSE 
    365          rnf (:,:) = 0._wp 
    366       ENDIF 
    367  
    368       IF( ln_dynbbl ) THEN         ! eiv & bbl 
    369                  jf_ubl  = jfld + 1 ;        jf_vbl  = jfld + 2 ;  jfld = jf_vbl 
    370            slf_d(jf_ubl) = sn_ubl  ;   slf_d(jf_vbl) = sn_vbl 
    371       ENDIF 
    372  
    373  
     286         rnf(:,:) = 0._wp 
     287      ENDIF 
     288 
     289   
    374290      ALLOCATE( sf_dyn(jfld), STAT=ierr )         ! set sf structure 
    375       IF( ierr > 0 ) THEN 
     291      IF( ierr > 0 )  THEN 
    376292         CALL ctl_stop( 'dta_dyn: unable to allocate sf structure' )   ;   RETURN 
    377293      ENDIF 
    378294      !                                         ! fill sf with slf_i and control print 
    379295      CALL fld_fill( sf_dyn, slf_d, cn_dir, 'dta_dyn_init', 'Data in file', 'namdta_dyn' ) 
     296      ! 
    380297      ! Open file for each variable to get his number of dimension 
    381298      DO ifpr = 1, jfld 
    … …  
    401318            ALLOCATE( uslpdta (jpi,jpj,jpk,2), vslpdta (jpi,jpj,jpk,2),    & 
    402319            &         wslpidta(jpi,jpj,jpk,2), wslpjdta(jpi,jpj,jpk,2), STAT=ierr2 ) 
    403          ELSE 
    404             ALLOCATE( uslpnow (jpi,jpj,jpk)  , vslpnow (jpi,jpj,jpk)  ,    & 
    405             &         wslpinow(jpi,jpj,jpk)  , wslpjnow(jpi,jpj,jpk)  , STAT=ierr2 ) 
    406          ENDIF  
    407          IF( ierr2 > 0 ) THEN 
    408             CALL ctl_stop( 'dta_dyn_init : unable to allocate slope arrays' )   ;   RETURN 
     320            ! 
     321            IF( ierr2 > 0 )  THEN 
     322               CALL ctl_stop( 'dta_dyn_init : unable to allocate slope arrays' )   ;   RETURN 
     323            ENDIF 
    409324         ENDIF 
    410325      ENDIF 
    411       IF( ln_dynwzv ) THEN                  ! slopes  
    412          IF( sf_dyn(jf_uwd)%ln_tint ) THEN      ! time interpolation 
    413             ALLOCATE( wdta(jpi,jpj,jpk,2), STAT=ierr3 ) 
    414          ELSE 
    415             ALLOCATE( wnow(jpi,jpj,jpk)  , STAT=ierr3 ) 
    416          ENDIF  
    417          IF( ierr3 > 0 ) THEN 
    418             CALL ctl_stop( 'dta_dyn_init : unable to allocate wdta arrays' )   ;   RETURN 
    419          ENDIF 
    420       ENDIF 
    421       ! 
    422       CALL dta_dyn( nit000 ) 
    423       ! 
    424    END SUBROUTINE dta_dyn_init 
    425  
    426  
    427    SUBROUTINE dta_dyn_wzv( pu, pv, pw ) 
    428       !!---------------------------------------------------------------------- 
    429       !!                    ***  ROUTINE wzv  *** 
    430       !! 
    431       !! ** Purpose :   Compute the now vertical velocity after the array swap 
    432       !! 
    433       !! ** Method  : - compute the now divergence given by : 
    434       !!         * z-coordinate ONLY !!!! 
    435       !!         hdiv = 1/(e1t*e2t) [ di(e2u  u) + dj(e1v  v) ] 
    436       !!     - Using the incompressibility hypothesis, the vertical 
    437       !!      velocity is computed by integrating the horizontal divergence 
    438       !!      from the bottom to the surface. 
    439       !!        The boundary conditions are w=0 at the bottom (no flux). 
    440       !!---------------------------------------------------------------------- 
    441       USE oce, ONLY:  zhdiv => hdivn 
    442       ! 
    443       REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) :: pu, pv    !:  horizontal velocities 
    444       REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(  out) :: pw        !:  vertical velocity 
    445       !! 
    446       INTEGER  ::  ji, jj, jk 
    447       REAL(wp) ::  zu, zu1, zv, zv1, zet 
    448       !!---------------------------------------------------------------------- 
    449       ! 
    450       ! Computation of vertical velocity using horizontal divergence 
    451       zhdiv(:,:,:) = 0._wp 
    452       DO jk = 1, jpkm1 
    453          DO jj = 2, jpjm1 
    454             DO ji = fs_2, fs_jpim1   ! vector opt. 
    455                zu  = pu(ji  ,jj  ,jk) * umask(ji  ,jj  ,jk) * e2u(ji  ,jj  ) * e3u_n(ji  ,jj  ,jk) 
    456                zu1 = pu(ji-1,jj  ,jk) * umask(ji-1,jj  ,jk) * e2u(ji-1,jj  ) * e3u_n(ji-1,jj  ,jk) 
    457                zv  = pv(ji  ,jj  ,jk) * vmask(ji  ,jj  ,jk) * e1v(ji  ,jj  ) * e3v_n(ji  ,jj  ,jk) 
    458                zv1 = pv(ji  ,jj-1,jk) * vmask(ji  ,jj-1,jk) * e1v(ji  ,jj-1) * e3v_n(ji  ,jj-1,jk) 
    459                zhdiv(ji,jj,jk) = ( zu - zu1 + zv - zv1 ) * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk)  
     326      ! 
     327      IF( .NOT.ln_linssh ) THEN 
     328        IF( .NOT. sf_dyn(jf_uwd)%ln_clim .AND. ln_rsttr .AND.    &                     ! Restart: read in restart file 
     329           iom_varid( numrtr, 'sshn', ldstop = .FALSE. ) > 0 ) THEN 
     330           IF(lwp) WRITE(numout,*) ' sshn forcing fields read in the restart file for initialisation' 
     331           CALL iom_get( numrtr, jpdom_autoglo, 'sshn', sshn(:,:)   ) 
     332           CALL iom_get( numrtr, jpdom_autoglo, 'sshb', sshb(:,:)   ) 
     333        ELSE 
     334           IF(lwp) WRITE(numout,*) ' sshn forcing fields read in the restart file for initialisation' 
     335           CALL iom_open( 'restart', inum ) 
     336           CALL iom_get( inum, jpdom_autoglo, 'sshn', sshn(:,:)   ) 
     337           CALL iom_get( inum, jpdom_autoglo, 'sshb', sshb(:,:)   ) 
     338           CALL iom_close( inum )                                        ! close file 
     339        ENDIF 
     340        ! 
     341        DO jk = 1, jpkm1 
     342           e3t_n(:,:,jk) = e3t_0(:,:,jk) * ( 1._wp + sshn(:,:) * tmask(:,:,1) / ( ht_0(:,:) + 1.0 - tmask(:,:,1) ) ) 
     343        ENDDO 
     344        e3t_a(:,:,jpk) = e3t_0(:,:,jpk) 
     345 
     346        ! Horizontal scale factor interpolations 
     347        ! -------------------------------------- 
     348        CALL dom_vvl_interpol( e3t_n(:,:,:), e3u_n(:,:,:), 'U' ) 
     349        CALL dom_vvl_interpol( e3t_n(:,:,:), e3v_n(:,:,:), 'V' ) 
     350 
     351        ! Vertical scale factor interpolations 
     352        ! ------------------------------------ 
     353        CALL dom_vvl_interpol( e3t_n(:,:,:), e3w_n(:,:,:), 'W' ) 
     354   
     355        e3t_b(:,:,:)  = e3t_n(:,:,:) 
     356        e3u_b(:,:,:)  = e3u_n(:,:,:) 
     357        e3v_b(:,:,:)  = e3v_n(:,:,:) 
     358 
     359        ! t- and w- points depth 
     360        ! ---------------------- 
     361        gdept_n(:,:,1) = 0.5_wp * e3w_n(:,:,1) 
     362        gdepw_n(:,:,1) = 0.0_wp 
     363 
     364        DO jk = 2, jpk 
     365           DO jj = 1,jpj 
     366              DO ji = 1,jpi 
     367                !    zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk))   ! 0 everywhere 
     368                !    tmask = wmask, ie everywhere expect at jk = mikt 
     369                                                                   ! 1 for jk = 
     370                                                                   ! mikt 
     371                 zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk)) 
     372                 gdepw_n(ji,jj,jk) = gdepw_n(ji,jj,jk-1) + e3t_n(ji,jj,jk-1) 
     373                 gdept_n(ji,jj,jk) =      zcoef  * ( gdepw_n(ji,jj,jk  ) + 0.5 * e3w_n(ji,jj,jk))  & 
     374                     &                + (1-zcoef) * ( gdept_n(ji,jj,jk-1) + e3w_n(ji,jj,jk)) 
     375              END DO 
     376           END DO 
     377        END DO 
     378 
     379        gdept_b(:,:,:) = gdept_n(:,:,:) 
     380        gdepw_b(:,:,:) = gdepw_n(:,:,:) 
     381        ! 
     382      ENDIF 
     383      ! 
     384      IF( ln_dynrnf .AND. ln_dynrnf_depth ) THEN       ! read depht over which runoffs are distributed 
     385         IF(lwp) WRITE(numout,*)  
     386         IF(lwp) WRITE(numout,*) ' read in the file depht over which runoffs are distributed' 
     387         CALL iom_open ( "runoffs", inum )                           ! open file 
     388         CALL iom_get  ( inum, jpdom_data, 'rodepth', h_rnf )   ! read the river mouth array 
     389         CALL iom_close( inum )                                        ! close file 
     390         ! 
     391         nk_rnf(:,:) = 0                               ! set the number of level over which river runoffs are applied 
     392         DO jj = 1, jpj 
     393            DO ji = 1, jpi 
     394               IF( h_rnf(ji,jj) > 0._wp ) THEN 
     395                  jk = 2 
     396                  DO WHILE ( jk /= mbkt(ji,jj) .AND. gdept_0(ji,jj,jk) < h_rnf(ji,jj) ) ;  jk = jk + 1 
     397                  END DO 
     398                  nk_rnf(ji,jj) = jk 
     399               ELSEIF( h_rnf(ji,jj) == -1._wp   ) THEN   ;  nk_rnf(ji,jj) = 1 
     400               ELSEIF( h_rnf(ji,jj) == -999._wp ) THEN   ;  nk_rnf(ji,jj) = mbkt(ji,jj) 
     401               ELSE 
     402                  CALL ctl_stop( 'sbc_rnf_init: runoff depth not positive, and not -999 or -1, rnf value in file fort.999'  ) 
     403                  WRITE(999,*) 'ji, jj, h_rnf(ji,jj) :', ji, jj, h_rnf(ji,jj) 
     404               ENDIF 
    460405            END DO 
    461406         END DO 
     407         DO jj = 1, jpj                                ! set the associated depth 
     408            DO ji = 1, jpi 
     409               h_rnf(ji,jj) = 0._wp 
     410               DO jk = 1, nk_rnf(ji,jj) 
     411                  h_rnf(ji,jj) = h_rnf(ji,jj) + e3t_n(ji,jj,jk) 
     412               END DO 
     413            END DO 
     414         END DO 
     415      ELSE                                       ! runoffs applied at the surface 
     416         nk_rnf(:,:) = 1 
     417         h_rnf (:,:) = e3t_n(:,:,1) 
     418      ENDIF 
     419      nkrnf_max = MAXVAL( nk_rnf(:,:) ) 
     420      hrnf_max = MAXVAL( h_rnf(:,:) ) 
     421      IF( lk_mpp )  THEN 
     422         CALL mpp_max( nkrnf_max )                 ! max over the  global domain 
     423         CALL mpp_max( hrnf_max )                 ! max over the  global domain 
     424      ENDIF 
     425      IF(lwp) WRITE(numout,*) ' ' 
     426      IF(lwp) WRITE(numout,*) ' max depht of runoff : ', hrnf_max,'    max level  : ', nkrnf_max 
     427      IF(lwp) WRITE(numout,*) ' ' 
     428      ! 
     429      CALL dta_dyn( nit000 ) 
     430      ! 
     431   END SUBROUTINE dta_dyn_init 
     432 
     433   SUBROUTINE dta_dyn_swp( kt ) 
     434     !!--------------------------------------------------------------------- 
     435      !!                    ***  ROUTINE dta_dyn_swp  *** 
     436      !! 
     437      !! ** Purpose : Swap and the data and compute the vertical scale factor at U/V/W point 
     438      !!              and the depht 
     439      !! 
     440      !!--------------------------------------------------------------------- 
     441      INTEGER, INTENT(in) :: kt       ! time step 
     442      INTEGER             :: ji, jj, jk 
     443      REAL(wp)            :: zcoef 
     444      ! 
     445      !!--------------------------------------------------------------------- 
     446 
     447      IF( kt == nit000 ) THEN 
     448         IF(lwp) WRITE(numout,*) 
     449         IF(lwp) WRITE(numout,*) 'ssh_swp : Asselin time filter and swap of sea surface height' 
     450         IF(lwp) WRITE(numout,*) '~~~~~~~ ' 
     451      ENDIF 
     452 
     453      sshb(:,:) = sshn(:,:) + atfp * ( sshb(:,:) - 2 * sshn(:,:) + ssha(:,:))  ! before <-- now filtered 
     454      sshn(:,:) = ssha(:,:) 
     455 
     456      e3t_n(:,:,:) = e3t_a(:,:,:) 
     457 
     458      ! Reconstruction of all vertical scale factors at now and before time steps 
     459      ! ============================================================================= 
     460 
     461      ! Horizontal scale factor interpolations 
     462      ! -------------------------------------- 
     463      CALL dom_vvl_interpol( e3t_n(:,:,:), e3u_n(:,:,:), 'U' ) 
     464      CALL dom_vvl_interpol( e3t_n(:,:,:), e3v_n(:,:,:), 'V' ) 
     465 
     466      ! Vertical scale factor interpolations 
     467      ! ------------------------------------ 
     468      CALL dom_vvl_interpol( e3t_n(:,:,:), e3w_n (:,:,:), 'W' ) 
     469 
     470      e3t_b(:,:,:)  = e3t_n(:,:,:) 
     471      e3u_b(:,:,:)  = e3u_n(:,:,:) 
     472      e3v_b(:,:,:)  = e3v_n(:,:,:) 
     473 
     474      ! t- and w- points depth 
     475      ! ---------------------- 
     476      gdept_n(:,:,1) = 0.5_wp * e3w_n(:,:,1) 
     477      gdepw_n(:,:,1) = 0.0_wp 
     478 
     479      DO jk = 2, jpk 
     480         DO jj = 1,jpj 
     481            DO ji = 1,jpi 
     482                 zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk)) 
     483                 gdepw_n(ji,jj,jk) = gdepw_n(ji,jj,jk-1) + e3t_n(ji,jj,jk-1) 
     484                 gdept_n(ji,jj,jk) =      zcoef  * ( gdepw_n(ji,jj,jk  ) + 0.5 * e3w_n(ji,jj,jk))  & 
     485                     &                + (1-zcoef) * ( gdept_n(ji,jj,jk-1) + e3w_n(ji,jj,jk)) 
     486              END DO 
     487           END DO 
     488        END DO 
     489 
     490      gdept_b(:,:,:) = gdept_n(:,:,:) 
     491      gdepw_b(:,:,:) = gdepw_n(:,:,:) 
     492 
     493      ! 
     494   END SUBROUTINE dta_dyn_swp 
     495 
     496   SUBROUTINE dta_dyn_ssh( kt, phdivtr, psshb,  pemp, pssha, pe3ta ) 
     497      !!---------------------------------------------------------------------- 
     498      !!                ***  ROUTINE dta_dyn_wzv  *** 
     499      !!                    
     500      !! ** Purpose :   compute the after ssh (ssha) and the now vertical velocity 
     501      !! 
     502      !! ** Method  : Using the incompressibility hypothesis,  
     503      !!        - the ssh increment is computed by integrating the horizontal divergence  
     504      !!          and multiply by the time step. 
     505      !! 
     506      !!        - compute the after scale factor : repartition of ssh INCREMENT proportionnaly 
     507      !!                                           to the level thickness ( z-star case ) 
     508      !! 
     509      !!        - the vertical velocity is computed by integrating the horizontal divergence   
     510      !!          from the bottom to the surface minus the scale factor evolution. 
     511      !!          The boundary conditions are w=0 at the bottom (no flux) 
     512      !! 
     513      !! ** action  :   ssha / e3t_a / wn 
     514      !! 
     515      !! Reference  : Leclair, M., and G. Madec, 2009, Ocean Modelling. 
     516      !!---------------------------------------------------------------------- 
     517      !! * Arguments 
     518      INTEGER,                                   INTENT(in )    :: kt        !  time-step 
     519      REAL(wp), DIMENSION(jpi,jpj,jpk)          , INTENT(in )   :: phdivtr   ! horizontal divergence transport 
     520      REAL(wp), DIMENSION(jpi,jpj)    , OPTIONAL, INTENT(in )   :: psshb     ! now ssh 
     521      REAL(wp), DIMENSION(jpi,jpj)    , OPTIONAL, INTENT(in )   :: pemp      ! evaporation minus precipitation 
     522      REAL(wp), DIMENSION(jpi,jpj)    , OPTIONAL, INTENT(inout) :: pssha     ! after ssh 
     523      REAL(wp), DIMENSION(jpi,jpj,jpk), OPTIONAL, INTENT(out)   :: pe3ta     ! after vertical scale factor 
     524      !! * Local declarations 
     525      INTEGER                       :: jk 
     526      REAL(wp), DIMENSION(jpi,jpj)  :: zhdiv   
     527      REAL(wp)  :: z2dt   
     528      !!---------------------------------------------------------------------- 
     529       
     530      ! 
     531      z2dt = 2._wp * rdt 
     532      ! 
     533      zhdiv(:,:) = 0._wp 
     534      DO jk = 1, jpkm1 
     535         zhdiv(:,:) = zhdiv(:,:) +  phdivtr(:,:,jk) * tmask(:,:,jk) 
    462536      END DO 
    463       !                              !  update the horizontal divergence with the runoff inflow 
    464       IF( ln_dynrnf )   zhdiv(:,:,1) = zhdiv(:,:,1) - rnf(:,:) * r1_rau0 / e3t_n(:,:,1) 
    465       ! 
    466       CALL lbc_lnk( zhdiv, 'T', 1. )      ! Lateral boundary conditions on zhdiv 
    467       ! computation of vertical velocity from the bottom 
    468       pw(:,:,jpk) = 0._wp 
    469       DO jk = jpkm1, 1, -1 
    470          pw(:,:,jk) = pw(:,:,jk+1) - e3t_n(:,:,jk) * zhdiv(:,:,jk) 
     537      !                                                ! Sea surface  elevation time-stepping 
     538      pssha(:,:) = ( psshb(:,:) - z2dt * ( r1_rau0 * pemp(:,:)  + zhdiv(:,:) ) ) * ssmask(:,:) 
     539      !                                                 !  
     540      !                                                 ! After acale factors at t-points ( z_star coordinate ) 
     541      DO jk = 1, jpkm1 
     542        pe3ta(:,:,jk) = e3t_0(:,:,jk) * ( 1._wp + pssha(:,:) * tmask(:,:,1) / ( ht_0(:,:) + 1.0 - tmask(:,:,1) ) ) 
    471543      END DO 
    472544      ! 
    473    END SUBROUTINE dta_dyn_wzv 
    474  
    475    SUBROUTINE dta_dyn_slp( kt, pts, puslp, pvslp, pwslpi, pwslpj ) 
     545   END SUBROUTINE dta_dyn_ssh 
     546 
     547 
     548   SUBROUTINE dta_dyn_hrnf 
     549      !!---------------------------------------------------------------------- 
     550      !!                  ***  ROUTINE sbc_rnf  *** 
     551      !! 
     552      !! ** Purpose :   update the horizontal divergence with the runoff inflow 
     553      !! 
     554      !! ** Method  : 
     555      !!                CAUTION : rnf is positive (inflow) decreasing the 
     556      !!                          divergence and expressed in m/s 
     557      !! 
     558      !! ** Action  :   phdivn   decreased by the runoff inflow 
     559      !!---------------------------------------------------------------------- 
     560      !! 
     561      INTEGER  ::   ji, jj, jk   ! dummy loop indices 
     562      !!---------------------------------------------------------------------- 
     563      ! 
     564      DO jj = 1, jpj                   ! update the depth over which runoffs are distributed 
     565         DO ji = 1, jpi 
     566            h_rnf(ji,jj) = 0._wp 
     567            DO jk = 1, nk_rnf(ji,jj)                           ! recalculates h_rnf to be the depth in metres 
     568                h_rnf(ji,jj) = h_rnf(ji,jj) + e3t_n(ji,jj,jk)   ! to the bottom of the relevant grid box 
     569            END DO 
     570        END DO 
     571      END DO 
     572      ! 
     573   END SUBROUTINE dta_dyn_hrnf 
     574 
     575 
     576 
     577   SUBROUTINE dta_dyn_slp( kt ) 
     578      !!--------------------------------------------------------------------- 
     579      !!                    ***  ROUTINE dta_dyn_slp  *** 
     580      !! 
     581      !! ** Purpose : Computation of slope 
     582      !! 
     583      !!--------------------------------------------------------------------- 
     584      USE oce, ONLY:  zts => tsa  
     585      ! 
     586      INTEGER,  INTENT(in) :: kt       ! time step 
     587      ! 
     588      INTEGER  ::   ji, jj     ! dummy loop indices 
     589      REAL(wp) ::   ztinta     ! ratio applied to after  records when doing time interpolation 
     590      REAL(wp) ::   ztintb     ! ratio applied to before records when doing time interpolation 
     591      INTEGER  ::   iswap  
     592      REAL(wp), POINTER, DIMENSION(:,:,:) :: zuslp, zvslp, zwslpi, zwslpj 
     593      !!--------------------------------------------------------------------- 
     594      ! 
     595      CALL wrk_alloc(jpi, jpj, jpk, zuslp, zvslp, zwslpi, zwslpj ) 
     596      ! 
     597      IF( sf_dyn(jf_tem)%ln_tint ) THEN    ! Computes slopes (here avt is used as workspace)                        
     598         IF( kt == nit000 ) THEN 
     599            IF(lwp) WRITE(numout,*) ' Compute new slopes at kt = ', kt 
     600            zts(:,:,:,jp_tem) = sf_dyn(jf_tem)%fdta(:,:,:,1) * tmask(:,:,:)   ! temperature 
     601            zts(:,:,:,jp_sal) = sf_dyn(jf_sal)%fdta(:,:,:,1) * tmask(:,:,:)   ! salinity  
     602            avt(:,:,:)        = sf_dyn(jf_avt)%fdta(:,:,:,1) * tmask(:,:,:)   ! vertical diffusive coef. 
     603            CALL compute_slopes( kt, zts, zuslp, zvslp, zwslpi, zwslpj ) 
     604            uslpdta (:,:,:,1) = zuslp (:,:,:)  
     605            vslpdta (:,:,:,1) = zvslp (:,:,:)  
     606            wslpidta(:,:,:,1) = zwslpi(:,:,:)  
     607            wslpjdta(:,:,:,1) = zwslpj(:,:,:)  
     608            ! 
     609            zts(:,:,:,jp_tem) = sf_dyn(jf_tem)%fdta(:,:,:,2) * tmask(:,:,:)   ! temperature 
     610            zts(:,:,:,jp_sal) = sf_dyn(jf_sal)%fdta(:,:,:,2) * tmask(:,:,:)   ! salinity  
     611            avt(:,:,:)        = sf_dyn(jf_avt)%fdta(:,:,:,2) * tmask(:,:,:)   ! vertical diffusive coef. 
     612            CALL compute_slopes( kt, zts, zuslp, zvslp, zwslpi, zwslpj ) 
     613            uslpdta (:,:,:,2) = zuslp (:,:,:)  
     614            vslpdta (:,:,:,2) = zvslp (:,:,:)  
     615            wslpidta(:,:,:,2) = zwslpi(:,:,:)  
     616            wslpjdta(:,:,:,2) = zwslpj(:,:,:)  
     617         ELSE 
     618           !  
     619           iswap = 0 
     620           IF( sf_dyn(jf_tem)%nrec_a(2) - nprevrec /= 0 )  iswap = 1 
     621           IF( nsecdyn > sf_dyn(jf_tem)%nrec_b(2) .AND. iswap == 1 )  THEN    ! read/update the after data 
     622              IF(lwp) WRITE(numout,*) ' Compute new slopes at kt = ', kt 
     623              uslpdta (:,:,:,1) =  uslpdta (:,:,:,2)         ! swap the data 
     624              vslpdta (:,:,:,1) =  vslpdta (:,:,:,2)   
     625              wslpidta(:,:,:,1) =  wslpidta(:,:,:,2)  
     626              wslpjdta(:,:,:,1) =  wslpjdta(:,:,:,2)  
     627              ! 
     628              zts(:,:,:,jp_tem) = sf_dyn(jf_tem)%fdta(:,:,:,2) * tmask(:,:,:)   ! temperature 
     629              zts(:,:,:,jp_sal) = sf_dyn(jf_sal)%fdta(:,:,:,2) * tmask(:,:,:)   ! salinity  
     630              avt(:,:,:)        = sf_dyn(jf_avt)%fdta(:,:,:,2) * tmask(:,:,:)   ! vertical diffusive coef. 
     631              CALL compute_slopes( kt, zts, zuslp, zvslp, zwslpi, zwslpj ) 
     632              ! 
     633              uslpdta (:,:,:,2) = zuslp (:,:,:)  
     634              vslpdta (:,:,:,2) = zvslp (:,:,:)  
     635              wslpidta(:,:,:,2) = zwslpi(:,:,:)  
     636              wslpjdta(:,:,:,2) = zwslpj(:,:,:)  
     637            ENDIF 
     638         ENDIF 
     639      ENDIF 
     640      ! 
     641      IF( sf_dyn(jf_tem)%ln_tint )  THEN 
     642         ztinta =  REAL( nsecdyn - sf_dyn(jf_tem)%nrec_b(2), wp )  & 
     643            &    / REAL( sf_dyn(jf_tem)%nrec_a(2) - sf_dyn(jf_tem)%nrec_b(2), wp ) 
     644         ztintb =  1. - ztinta 
     645         IF( l_ldfslp .AND. .NOT.lk_c1d ) THEN    ! Computes slopes (here avt is used as workspace) 
     646            uslp (:,:,:) = ztintb * uslpdta (:,:,:,1)  + ztinta * uslpdta (:,:,:,2)   
     647            vslp (:,:,:) = ztintb * vslpdta (:,:,:,1)  + ztinta * vslpdta (:,:,:,2)   
     648            wslpi(:,:,:) = ztintb * wslpidta(:,:,:,1)  + ztinta * wslpidta(:,:,:,2)   
     649            wslpj(:,:,:) = ztintb * wslpjdta(:,:,:,1)  + ztinta * wslpjdta(:,:,:,2)   
     650         ENDIF 
     651      ELSE 
     652         zts(:,:,:,jp_tem) = sf_dyn(jf_tem)%fnow(:,:,:) * tmask(:,:,:)   ! temperature 
     653         zts(:,:,:,jp_sal) = sf_dyn(jf_sal)%fnow(:,:,:) * tmask(:,:,:)   ! salinity  
     654         avt(:,:,:)        = sf_dyn(jf_avt)%fnow(:,:,:) * tmask(:,:,:)   ! vertical diffusive coef. 
     655         CALL compute_slopes( kt, zts, zuslp, zvslp, zwslpi, zwslpj ) 
     656         ! 
     657         IF( l_ldfslp .AND. .NOT.lk_c1d ) THEN    ! Computes slopes (here avt is used as workspace) 
     658            uslp (:,:,:) = zuslp (:,:,:) 
     659            vslp (:,:,:) = zvslp (:,:,:) 
     660            wslpi(:,:,:) = zwslpi(:,:,:) 
     661            wslpj(:,:,:) = zwslpj(:,:,:) 
     662         ENDIF 
     663      ENDIF 
     664      ! 
     665      CALL wrk_dealloc(jpi, jpj, jpk, zuslp, zvslp, zwslpi, zwslpj ) 
     666      ! 
     667   END SUBROUTINE dta_dyn_slp 
     668 
     669   SUBROUTINE compute_slopes( kt, pts, puslp, pvslp, pwslpi, pwslpj ) 
    476670      !!--------------------------------------------------------------------- 
    477671      !!                    ***  ROUTINE dta_dyn_slp  *** 
    … …  
    487681      REAL(wp), DIMENSION(jpi,jpj,jpk)     , INTENT(out) :: pwslpj   ! meridional diapycnal slopes 
    488682      !!--------------------------------------------------------------------- 
    489       IF( l_ldfslp .AND. .NOT.lk_c1d ) THEN    ! Computes slopes (here avt is used as workspace)                        
     683      IF( l_ldfslp .AND. .NOT.lk_c1d ) THEN    ! Computes slopes (here avt is used as workspace) 
    490684         CALL eos    ( pts, rhd, rhop, gdept_0(:,:,:) ) 
    491685         CALL eos_rab( pts, rab_n )       ! now local thermal/haline expension ratio at T-points 
    … …  
    497691         &                                        rhd, gru , grv    )  ! of t, s, rd at the last ocean level 
    498692      IF( ln_zps .AND.        ln_isfcav)                            & 
    499          &            CALL zps_hde_isf( kt, jpts, pts, gtsu, gtsv, gtui, gtvi,  &  ! Partial steps for top cell (ISF) 
    500          &                                        rhd, gru , grv , grui, grvi   )  ! of t, s, rd at the first ocean level 
     693         &            CALL zps_hde_isf( kt, jpts, pts, gtsu, gtsv, gtui, gtvi, &  ! Partial steps for top cell (ISF) 
     694         &                                        rhd, gru , grv , grui, grvi )  ! of t, s, rd at the first ocean level 
    501695 
    502696         rn2b(:,:,:) = rn2(:,:,:)         ! need for zdfmxl 
    503697         CALL zdf_mxl( kt )            ! mixed layer depth 
    504698         CALL ldf_slp( kt, rhd, rn2 )  ! slopes 
    505          puslp (:,:,:) = uslp (:,:,:)  
    506          pvslp (:,:,:) = vslp (:,:,:)  
    507          pwslpi(:,:,:) = wslpi(:,:,:)  
    508          pwslpj(:,:,:) = wslpj(:,:,:)  
     699         puslp (:,:,:) = uslp (:,:,:) 
     700         pvslp (:,:,:) = vslp (:,:,:) 
     701         pwslpi(:,:,:) = wslpi(:,:,:) 
     702         pwslpj(:,:,:) = wslpj(:,:,:) 
    509703     ELSE 
    510704         puslp (:,:,:) = 0.            ! to avoid warning when compiling 
    … …  
    514708     ENDIF 
    515709      ! 
    516    END SUBROUTINE dta_dyn_slp 
     710   END SUBROUTINE compute_slopes 
    517711   !!====================================================================== 
    518712END MODULE dtadyn 
Note: See TracChangeset for help on using the changeset viewer.