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.
dtadyn.F90 in NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OFF – NEMO

source: NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OFF/dtadyn.F90 @ 12779

Last change on this file since 12779 was 12779, checked in by techene, 4 years ago

modify TOP to be able to run with key_qco i.e. remove e3. gdep. and use a substitute instead (# include "domzgr_substitute.h90")

  • Property svn:keywords set to Id
File size: 42.2 KB
RevLine 
[325]1MODULE dtadyn
2   !!======================================================================
3   !!                       ***  MODULE  dtadyn  ***
[2528]4   !! Off-line : interpolation of the physical fields
5   !!======================================================================
6   !! History :   OPA  ! 1992-01 (M. Imbard) Original code
7   !!             8.0  ! 1998-04 (L.Bopp MA Foujols) slopes for isopyc.
8   !!              -   ! 1998-05 (L. Bopp) read output of coupled run
9   !!             8.2  ! 2001-01 (M. Levy et M. Benjelloul) add netcdf FORMAT
10   !!   NEMO      1.0  ! 2005-03 (O. Aumont and A. El Moussaoui) F90
11   !!              -   ! 2005-12 (C. Ethe) Adapted for DEGINT
12   !!             3.0  ! 2007-06 (C. Ethe) use of iom module
13   !!             3.3  ! 2010-11 (C. Ethe) Full reorganization of the off-line: phasing with the on-line
[3294]14   !!             3.4  ! 2011-05 (C. Ethe) Use of fldread
[12377]15   !!             4.1  ! 2019-08 (A. Coward, D. Storkey) split dta_dyn_sf_swp -> dta_dyn_sf_atf and dta_dyn_sf_interp
[2528]16   !!----------------------------------------------------------------------
[325]17
18   !!----------------------------------------------------------------------
[3294]19   !!   dta_dyn_init : initialization, namelist read, and SAVEs control
[325]20   !!   dta_dyn      : Interpolation of the fields
21   !!----------------------------------------------------------------------
22   USE oce             ! ocean dynamics and tracers variables
[2528]23   USE c1d             ! 1D configuration: lk_c1d
24   USE dom_oce         ! ocean domain: variables
[12779]25#if ! defined key_qco
[7646]26   USE domvvl          ! variable volume
[12779]27#else
28   USE domqco
29#endif
[2528]30   USE zdf_oce         ! ocean vertical physics: variables
31   USE sbc_oce         ! surface module: variables
[3294]32   USE trc_oce         ! share ocean/biogeo variables
[325]33   USE phycst          ! physical constants
[2528]34   USE trabbl          ! active tracer: bottom boundary layer
35   USE ldfslp          ! lateral diffusion: iso-neutral slopes
[7646]36   USE sbcrnf          ! river runoffs
37   USE ldftra          ! ocean tracer   lateral physics
[2528]38   USE zdfmxl          ! vertical physics: mixed layer depth
39   USE eosbn2          ! equation of state - Brunt Vaisala frequency
[325]40   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
[2528]41   USE zpshde          ! z-coord. with partial steps: horizontal derivatives
42   USE in_out_manager  ! I/O manager
43   USE iom             ! I/O library
[325]44   USE lib_mpp         ! distributed memory computing library
[3294]45   USE prtctl          ! print control
46   USE fldread         ! read input fields
47   USE timing          ! Timing
[7646]48   USE trc, ONLY : ln_rsttr, numrtr, numrtw, lrst_trc
[325]49
50   IMPLICIT NONE
51   PRIVATE
52
[12377]53   PUBLIC   dta_dyn_init       ! called by nemo_init
54   PUBLIC   dta_dyn            ! called by nemo_gcm
55   PUBLIC   dta_dyn_sed_init   ! called by nemo_init
56   PUBLIC   dta_dyn_sed        ! called by nemo_gcm
57   PUBLIC   dta_dyn_atf        ! called by nemo_gcm
[12779]58#if ! defined key_qco
[12377]59   PUBLIC   dta_dyn_sf_interp  ! called by nemo_gcm
[12779]60#endif
[325]61
[7646]62   CHARACTER(len=100) ::   cn_dir          !: Root directory for location of ssr files
63   LOGICAL            ::   ln_dynrnf       !: read runoff data in file (T) or set to zero (F)
64   LOGICAL            ::   ln_dynrnf_depth       !: read runoff data in file (T) or set to zero (F)
65   REAL(wp)           ::   fwbcorr
[325]66
[7646]67
68   INTEGER  , PARAMETER ::   jpfld = 20     ! maximum number of fields to read
[3294]69   INTEGER  , SAVE      ::   jf_tem         ! index of temperature
70   INTEGER  , SAVE      ::   jf_sal         ! index of salinity
[7646]71   INTEGER  , SAVE      ::   jf_uwd         ! index of u-transport
72   INTEGER  , SAVE      ::   jf_vwd         ! index of v-transport
73   INTEGER  , SAVE      ::   jf_wwd         ! index of v-transport
[3294]74   INTEGER  , SAVE      ::   jf_avt         ! index of Kz
75   INTEGER  , SAVE      ::   jf_mld         ! index of mixed layer deptht
76   INTEGER  , SAVE      ::   jf_emp         ! index of water flux
[7646]77   INTEGER  , SAVE      ::   jf_empb        ! index of water flux
[3294]78   INTEGER  , SAVE      ::   jf_qsr         ! index of solar radiation
79   INTEGER  , SAVE      ::   jf_wnd         ! index of wind speed
80   INTEGER  , SAVE      ::   jf_ice         ! index of sea ice cover
[4570]81   INTEGER  , SAVE      ::   jf_rnf         ! index of river runoff
[7646]82   INTEGER  , SAVE      ::   jf_fmf         ! index of downward salt flux
[3294]83   INTEGER  , SAVE      ::   jf_ubl         ! index of u-bbl coef
84   INTEGER  , SAVE      ::   jf_vbl         ! index of v-bbl coef
[7646]85   INTEGER  , SAVE      ::   jf_div         ! index of e3t
[325]86
[7646]87
88   TYPE(FLD), ALLOCATABLE, SAVE, DIMENSION(:) :: sf_dyn  ! structure of input fields (file informations, fields read)
[3294]89   !                                               !
90   REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: uslpdta    ! zonal isopycnal slopes
91   REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: vslpdta    ! meridional isopycnal slopes
92   REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: wslpidta   ! zonal diapycnal slopes
93   REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: wslpjdta   ! meridional diapycnal slopes
[1735]94
[7646]95   INTEGER, SAVE  :: nprevrec, nsecdyn
[325]96
[12377]97   !! * Substitutions
98#  include "do_loop_substitute.h90"
[343]99   !!----------------------------------------------------------------------
[9598]100   !! NEMO/OFF 4.0 , NEMO Consortium (2018)
[2528]101   !! $Id$
[10068]102   !! Software governed by the CeCILL license (see ./LICENSE)
[343]103   !!----------------------------------------------------------------------
[325]104CONTAINS
105
[12377]106   SUBROUTINE dta_dyn( kt, Kbb, Kmm, Kaa )
[325]107      !!----------------------------------------------------------------------
108      !!                  ***  ROUTINE dta_dyn  ***
109      !!
[3294]110      !! ** Purpose :  Prepares dynamics and physics fields from a NEMO run
111      !!               for an off-line simulation of passive tracers
[325]112      !!
[3294]113      !! ** Method : calculates the position of data
114      !!             - computes slopes if needed
115      !!             - interpolates data if needed
[2528]116      !!----------------------------------------------------------------------
[12377]117      INTEGER, INTENT(in) ::   kt             ! ocean time-step index
118      INTEGER, INTENT(in) ::   Kbb, Kmm, Kaa  ! ocean time level indices
[9212]119      !
[7646]120      INTEGER             ::   ji, jj, jk
[9212]121      REAL(wp), ALLOCATABLE, DIMENSION(:,:)   ::   zemp
122      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   zhdivtr
[325]123      !!----------------------------------------------------------------------
[3294]124      !
[9124]125      IF( ln_timing )   CALL timing_start( 'dta_dyn')
[3294]126      !
[7646]127      nsecdyn = nsec_year + nsec1jan000   ! number of seconds between Jan. 1st 00h of nit000 year and the middle of time step
[3294]128      !
[7646]129      IF( kt == nit000 ) THEN    ;    nprevrec = 0
130      ELSE                       ;    nprevrec = sf_dyn(jf_tem)%nrec_a(2)
[325]131      ENDIF
[7646]132      CALL fld_read( kt, 1, sf_dyn )      !=  read data at kt time step   ==!
133      !
[12377]134      IF( l_ldfslp .AND. .NOT.lk_c1d )   CALL  dta_dyn_slp( kt, Kbb, Kmm )    ! Computation of slopes
[7646]135      !
[12779]136      ts    (:,:,:,jp_tem,Kmm) = sf_dyn(jf_tem)%fnow(:,:,:) * tmask(:,:,:)    ! temperature
137      ts    (:,:,:,jp_sal,Kmm) = sf_dyn(jf_sal)%fnow(:,:,:) * tmask(:,:,:)    ! salinity
138      wndm  (:,:)              = sf_dyn(jf_wnd)%fnow(:,:,1) * tmask(:,:,1)    ! wind speed - needed for gas exchange
139      fmmflx(:,:)              = sf_dyn(jf_fmf)%fnow(:,:,1) * tmask(:,:,1)    ! downward salt flux (v3.5+)
140      fr_i  (:,:)              = sf_dyn(jf_ice)%fnow(:,:,1) * tmask(:,:,1)    ! Sea-ice fraction
141      qsr   (:,:)              = sf_dyn(jf_qsr)%fnow(:,:,1) * tmask(:,:,1)    ! solar radiation
142      emp   (:,:)              = sf_dyn(jf_emp)%fnow(:,:,1) * tmask(:,:,1)    ! E-P
[7646]143      IF( ln_dynrnf ) THEN
[12779]144         rnf(:,:)              = sf_dyn(jf_rnf)%fnow(:,:,1) * tmask(:,:,1)    ! E-P
145         IF( ln_dynrnf_depth .AND. .NOT.ln_linssh )   CALL dta_dyn_hrnf( Kmm )
[3294]146      ENDIF
[325]147      !
[12779]148      uu(:,:,:,Kmm)            = sf_dyn(jf_uwd)%fnow(:,:,:) * umask(:,:,:)    ! effective u-transport
149      vv(:,:,:,Kmm)            = sf_dyn(jf_vwd)%fnow(:,:,:) * vmask(:,:,:)    ! effective v-transport
150      ww(:,:,:)                = sf_dyn(jf_wwd)%fnow(:,:,:) * tmask(:,:,:)    ! effective v-transport
[7646]151      !
152      IF( .NOT.ln_linssh ) THEN
[9212]153         ALLOCATE( zemp(jpi,jpj) , zhdivtr(jpi,jpj,jpk) )
[12779]154         zhdivtr(:,:,:) = sf_dyn(jf_div )%fnow(:,:,:) * tmask(:,:,:)    ! effective u-transport
[9212]155         emp_b  (:,:)   = sf_dyn(jf_empb)%fnow(:,:,1) * tmask(:,:,1)    ! E-P
156         zemp   (:,:)   = ( 0.5_wp * ( emp(:,:) + emp_b(:,:) ) + rnf(:,:) + fwbcorr ) * tmask(:,:,1)
[12779]157#if defined key_qco
158         CALL dta_dyn_ssh( kt, zhdivtr, ssh(:,:,Kbb), zemp, ssh(:,:,Kaa) )
159         CALL dom_qco_r3c( ssh(:,:,Kaa), r3t(:,:,Kaa), r3u(:,:,Kaa), r3v(:,:,Kaa) )
160#else
161         CALL dta_dyn_ssh( kt, zhdivtr, ssh(:,:,Kbb), zemp, ssh(:,:,Kaa), e3t(:,:,:,Kaa) )  !=  ssh, vertical scale factor
162#endif
[9212]163         DEALLOCATE( zemp , zhdivtr )
[7646]164         !                                           Write in the tracer restart file
[9212]165         !                                          *********************************
[7646]166         IF( lrst_trc ) THEN
[3827]167            IF(lwp) WRITE(numout,*)
[9212]168            IF(lwp) WRITE(numout,*) 'dta_dyn_ssh : ssh field written in tracer restart file at it= ', kt,' date= ', ndastp
169            IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
[12377]170            CALL iom_rstput( kt, nitrst, numrtw, 'sshn', ssh(:,:,Kaa) )
171            CALL iom_rstput( kt, nitrst, numrtw, 'sshb', ssh(:,:,Kmm) )
[1501]172         ENDIF
[3294]173      ENDIF
[325]174      !
[12377]175      CALL eos    ( ts(:,:,:,:,Kmm), rhd, rhop, gdept_0(:,:,:) ) ! In any case, we need rhop
176      CALL eos_rab( ts(:,:,:,:,Kmm), rab_n, Kmm )       ! now    local thermal/haline expension ratio at T-points
177      CALL bn2    ( ts(:,:,:,:,Kmm), rab_n, rn2, Kmm )  ! before Brunt-Vaisala frequency need for zdfmxl
[5131]178
[12377]179      rn2b(:,:,:) = rn2(:,:,:)         ! needed for zdfmxl
180      CALL zdf_mxl( kt, Kmm )          ! In any case, we need mxl
[2528]181      !
[9019]182      hmld(:,:)       = sf_dyn(jf_mld)%fnow(:,:,1) * tmask(:,:,1)    ! mixed layer depht
183      avt(:,:,:)      = sf_dyn(jf_avt)%fnow(:,:,:) * tmask(:,:,:)    ! vertical diffusive coefficient
[10213]184      avs(:,:,:)      = avt(:,:,:)
[7646]185      !
[9019]186      IF( ln_trabbl .AND. .NOT.lk_c1d ) THEN       ! diffusive Bottom boundary layer param
187         ahu_bbl(:,:) = sf_dyn(jf_ubl)%fnow(:,:,1) * umask(:,:,1)    ! bbl diffusive coef
188         ahv_bbl(:,:) = sf_dyn(jf_vbl)%fnow(:,:,1) * vmask(:,:,1)
189      ENDIF
[2762]190      !
[7646]191      !
[12377]192      CALL eos( ts(:,:,:,:,Kmm), rhd, rhop, gdept_0(:,:,:) ) ! In any case, we need rhop
[7646]193      !
[12377]194      IF(sn_cfctl%l_prtctl) THEN                 ! print control
195         CALL prt_ctl(tab3d_1=ts(:,:,:,jp_tem,Kmm), clinfo1=' tn      - : ', mask1=tmask,  kdim=jpk   )
196         CALL prt_ctl(tab3d_1=ts(:,:,:,jp_sal,Kmm), clinfo1=' sn      - : ', mask1=tmask,  kdim=jpk   )
197         CALL prt_ctl(tab3d_1=uu(:,:,:,Kmm)               , clinfo1=' uu(:,:,:,Kmm)      - : ', mask1=umask,  kdim=jpk   )
198         CALL prt_ctl(tab3d_1=vv(:,:,:,Kmm)               , clinfo1=' vv(:,:,:,Kmm)      - : ', mask1=vmask,  kdim=jpk   )
199         CALL prt_ctl(tab3d_1=ww               , clinfo1=' ww      - : ', mask1=tmask,  kdim=jpk   )
[9440]200         CALL prt_ctl(tab3d_1=avt              , clinfo1=' kz      - : ', mask1=tmask,  kdim=jpk   )
[7646]201         CALL prt_ctl(tab3d_1=uslp             , clinfo1=' slp  - u : ', tab3d_2=vslp, clinfo2=' v : ', kdim=jpk)
202         CALL prt_ctl(tab3d_1=wslpi            , clinfo1=' slp  - wi: ', tab3d_2=wslpj, clinfo2=' wj: ', kdim=jpk)
[2528]203      ENDIF
204      !
[9124]205      IF( ln_timing )   CALL timing_stop( 'dta_dyn')
[3294]206      !
[325]207   END SUBROUTINE dta_dyn
208
[2528]209
[12377]210   SUBROUTINE dta_dyn_init( Kbb, Kmm, Kaa )
[325]211      !!----------------------------------------------------------------------
[3294]212      !!                  ***  ROUTINE dta_dyn_init  ***
[325]213      !!
[3294]214      !! ** Purpose :   Initialisation of the dynamical data     
215      !! ** Method  : - read the data namdta_dyn namelist
[325]216      !!----------------------------------------------------------------------
[12377]217      INTEGER, INTENT(in) ::   Kbb, Kmm, Kaa  ! ocean time level indices
218      !
[3294]219      INTEGER  :: ierr, ierr0, ierr1, ierr2, ierr3   ! return error code
220      INTEGER  :: ifpr                               ! dummy loop indice
221      INTEGER  :: jfld                               ! dummy loop arguments
222      INTEGER  :: inum, idv, idimv                   ! local integer
[4147]223      INTEGER  :: ios                                ! Local integer output status for namelist read
[7646]224      INTEGER  :: ji, jj, jk
225      REAL(wp) :: zcoef
226      INTEGER  :: nkrnf_max
227      REAL(wp) :: hrnf_max
[3294]228      !!
[7646]229      CHARACTER(len=100)            ::  cn_dir        !   Root directory for location of core files
230      TYPE(FLD_N), DIMENSION(jpfld) ::  slf_d         ! array of namelist informations on the fields to read
231      TYPE(FLD_N) :: sn_uwd, sn_vwd, sn_wwd, sn_empb, sn_emp  ! informations about the fields to be read
232      TYPE(FLD_N) :: sn_tem , sn_sal , sn_avt   !   "                 "
233      TYPE(FLD_N) :: sn_mld, sn_qsr, sn_wnd , sn_ice , sn_fmf   !   "               "
234      TYPE(FLD_N) :: sn_ubl, sn_vbl, sn_rnf    !   "              "
235      TYPE(FLD_N) :: sn_div  ! informations about the fields to be read
[9212]236      !!
[7646]237      NAMELIST/namdta_dyn/cn_dir, ln_dynrnf, ln_dynrnf_depth,  fwbcorr, &
[9212]238         &                sn_uwd, sn_vwd, sn_wwd, sn_emp,               &
239         &                sn_avt, sn_tem, sn_sal, sn_mld , sn_qsr ,     &
240         &                sn_wnd, sn_ice, sn_fmf,                       &
241         &                sn_ubl, sn_vbl, sn_rnf,                       &
[7646]242         &                sn_empb, sn_div 
[9212]243      !!----------------------------------------------------------------------
[7646]244      !
[4147]245      READ  ( numnam_ref, namdta_dyn, IOSTAT = ios, ERR = 901)
[11536]246901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namdta_dyn in reference namelist' )
[4147]247      READ  ( numnam_cfg, namdta_dyn, IOSTAT = ios, ERR = 902 )
[11536]248902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namdta_dyn in configuration namelist' )
[4624]249      IF(lwm) WRITE ( numond, namdta_dyn )
[3294]250      !                                         ! store namelist information in an array
251      !                                         ! Control print
[325]252      IF(lwp) THEN
253         WRITE(numout,*)
[3294]254         WRITE(numout,*) 'dta_dyn : offline dynamics '
255         WRITE(numout,*) '~~~~~~~ '
256         WRITE(numout,*) '   Namelist namdta_dyn'
[7646]257         WRITE(numout,*) '      runoffs option enabled (T) or not (F)            ln_dynrnf        = ', ln_dynrnf
258         WRITE(numout,*) '      runoffs is spread in vertical                    ln_dynrnf_depth  = ', ln_dynrnf_depth
259         WRITE(numout,*) '      annual global mean of empmr for ssh correction   fwbcorr          = ', fwbcorr
[325]260         WRITE(numout,*)
261      ENDIF
[3294]262      !
[7646]263      jf_uwd  = 1     ;   jf_vwd  = 2    ;   jf_wwd = 3    ;   jf_emp = 4    ;   jf_avt = 5
264      jf_tem  = 6     ;   jf_sal  = 7    ;   jf_mld = 8    ;   jf_qsr = 9
265      jf_wnd  = 10    ;   jf_ice  = 11   ;   jf_fmf = 12   ;   jfld   = jf_fmf
[3294]266      !
[7646]267      slf_d(jf_uwd)  = sn_uwd    ;   slf_d(jf_vwd)  = sn_vwd   ;   slf_d(jf_wwd) = sn_wwd
268      slf_d(jf_emp)  = sn_emp    ;   slf_d(jf_avt)  = sn_avt
269      slf_d(jf_tem)  = sn_tem    ;   slf_d(jf_sal)  = sn_sal   ;   slf_d(jf_mld) = sn_mld
270      slf_d(jf_qsr)  = sn_qsr    ;   slf_d(jf_wnd)  = sn_wnd   ;   slf_d(jf_ice) = sn_ice
271      slf_d(jf_fmf)  = sn_fmf
[3294]272      !
[7646]273      IF( .NOT.ln_linssh ) THEN
[9212]274               jf_div  = jfld + 1   ;         jf_empb  = jfld + 2    ;   jfld = jf_empb
275         slf_d(jf_div) = sn_div     ;   slf_d(jf_empb) = sn_empb
[7646]276      ENDIF
277      !
[9019]278      IF( ln_trabbl ) THEN
[9212]279               jf_ubl  = jfld + 1   ;         jf_vbl  = jfld + 2     ;   jfld = jf_vbl
280         slf_d(jf_ubl) = sn_ubl     ;   slf_d(jf_vbl) = sn_vbl
[7646]281      ENDIF
282      !
[5385]283      IF( ln_dynrnf ) THEN
[9212]284               jf_rnf  = jfld + 1   ;     jfld  = jf_rnf
285         slf_d(jf_rnf) = sn_rnf   
[4570]286      ELSE
[7646]287         rnf(:,:) = 0._wp
[4570]288      ENDIF
289
[3294]290      ALLOCATE( sf_dyn(jfld), STAT=ierr )         ! set sf structure
[7646]291      IF( ierr > 0 )  THEN
[3294]292         CALL ctl_stop( 'dta_dyn: unable to allocate sf structure' )   ;   RETURN
293      ENDIF
[5768]294      !                                         ! fill sf with slf_i and control print
295      CALL fld_fill( sf_dyn, slf_d, cn_dir, 'dta_dyn_init', 'Data in file', 'namdta_dyn' )
[7646]296      !
[3294]297      ! Open file for each variable to get his number of dimension
298      DO ifpr = 1, jfld
[12377]299         CALL fld_def( sf_dyn(ifpr) )
300         CALL iom_open( sf_dyn(ifpr)%clname, sf_dyn(ifpr)%num )
[5768]301         idv   = iom_varid( sf_dyn(ifpr)%num , slf_d(ifpr)%clvar )        ! id of the variable sdjf%clvar
302         idimv = iom_file ( sf_dyn(ifpr)%num )%ndims(idv)                 ! number of dimension for variable sdjf%clvar
[12377]303         CALL iom_close( sf_dyn(ifpr)%num )                               ! close file if already open
[5768]304         ierr1=0
[3294]305         IF( idimv == 3 ) THEN    ! 2D variable
306                                      ALLOCATE( sf_dyn(ifpr)%fnow(jpi,jpj,1)    , STAT=ierr0 )
307            IF( slf_d(ifpr)%ln_tint ) ALLOCATE( sf_dyn(ifpr)%fdta(jpi,jpj,1,2)  , STAT=ierr1 )
308         ELSE                     ! 3D variable
309                                      ALLOCATE( sf_dyn(ifpr)%fnow(jpi,jpj,jpk)  , STAT=ierr0 )
310            IF( slf_d(ifpr)%ln_tint ) ALLOCATE( sf_dyn(ifpr)%fdta(jpi,jpj,jpk,2), STAT=ierr1 )
[2528]311         ENDIF
[3294]312         IF( ierr0 + ierr1 > 0 ) THEN
313            CALL ctl_stop( 'dta_dyn_init : unable to allocate sf_dyn array structure' )   ;   RETURN
314         ENDIF
315      END DO
[325]316      !
[5836]317      IF( l_ldfslp .AND. .NOT.lk_c1d ) THEN                  ! slopes
[3294]318         IF( sf_dyn(jf_tem)%ln_tint ) THEN      ! time interpolation
319            ALLOCATE( uslpdta (jpi,jpj,jpk,2), vslpdta (jpi,jpj,jpk,2),    &
320            &         wslpidta(jpi,jpj,jpk,2), wslpjdta(jpi,jpj,jpk,2), STAT=ierr2 )
[7646]321            !
322            IF( ierr2 > 0 )  THEN
323               CALL ctl_stop( 'dta_dyn_init : unable to allocate slope arrays' )   ;   RETURN
324            ENDIF
[3294]325         ENDIF
[2528]326      ENDIF
[7646]327      !
328      IF( .NOT.ln_linssh ) THEN
329        IF( .NOT. sf_dyn(jf_uwd)%ln_clim .AND. ln_rsttr .AND.    &                     ! Restart: read in restart file
330           iom_varid( numrtr, 'sshn', ldstop = .FALSE. ) > 0 ) THEN
[12377]331           IF(lwp) WRITE(numout,*) ' ssh(:,:,Kmm) forcing fields read in the restart file for initialisation'
332           CALL iom_get( numrtr, jpdom_autoglo, 'sshn', ssh(:,:,Kmm)   )
333           CALL iom_get( numrtr, jpdom_autoglo, 'sshb', ssh(:,:,Kbb)   )
[7646]334        ELSE
[12377]335           IF(lwp) WRITE(numout,*) ' ssh(:,:,Kmm) forcing fields read in the restart file for initialisation'
[7646]336           CALL iom_open( 'restart', inum )
[12377]337           CALL iom_get( inum, jpdom_autoglo, 'sshn', ssh(:,:,Kmm)   )
338           CALL iom_get( inum, jpdom_autoglo, 'sshb', ssh(:,:,Kbb)   )
[7646]339           CALL iom_close( inum )                                        ! close file
340        ENDIF
341        !
[12779]342#if defined key_qco
343        CALL dom_qco_r3c( ssh(:,:,Kbb), r3t(:,:,Kbb), r3u(:,:,Kbb), r3v(:,:,Kbb) )
344        CALL dom_qco_r3c( ssh(:,:,Kmm), r3t(:,:,Kmm), r3u(:,:,Kmm), r3v(:,:,Kmm) )
345#else
[7646]346        DO jk = 1, jpkm1
[12779]347           e3t(:,:,jk,Kmm) = e3t_0(:,:,jk) * ( 1._wp + ssh(:,:,Kmm) * r1_ht_0(:,:) * tmask(:,:,jk) )
348        END DO
[12377]349        e3t(:,:,jpk,Kaa) = e3t_0(:,:,jpk)
[7646]350
351        ! Horizontal scale factor interpolations
352        ! --------------------------------------
[12377]353        CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3u(:,:,:,Kmm), 'U' )
354        CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3v(:,:,:,Kmm), 'V' )
[7646]355
356        ! Vertical scale factor interpolations
357        ! ------------------------------------
[12377]358        CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3w(:,:,:,Kmm), 'W' )
[12779]359!!gm this should be computed from ssh(Kbb) 
[12377]360        e3t(:,:,:,Kbb)  = e3t(:,:,:,Kmm)
361        e3u(:,:,:,Kbb)  = e3u(:,:,:,Kmm)
362        e3v(:,:,:,Kbb)  = e3v(:,:,:,Kmm)
[7646]363
364        ! t- and w- points depth
365        ! ----------------------
[12377]366        gdept(:,:,1,Kmm) = 0.5_wp * e3w(:,:,1,Kmm)
367        gdepw(:,:,1,Kmm) = 0.0_wp
[7646]368
[12377]369        DO_3D_11_11( 2, jpk )
370          !    zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk))   ! 0 everywhere
371          !    tmask = wmask, ie everywhere expect at jk = mikt
372                                                             ! 1 for jk =
373                                                             ! mikt
374           zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk))
375           gdepw(ji,jj,jk,Kmm) = gdepw(ji,jj,jk-1,Kmm) + e3t(ji,jj,jk-1,Kmm)
376           gdept(ji,jj,jk,Kmm) =      zcoef  * ( gdepw(ji,jj,jk  ,Kmm) + 0.5 * e3w(ji,jj,jk,Kmm))  &
377               &                + (1-zcoef) * ( gdept(ji,jj,jk-1,Kmm) + e3w(ji,jj,jk,Kmm))
378        END_3D
[7646]379
[12377]380        gdept(:,:,:,Kbb) = gdept(:,:,:,Kmm)
381        gdepw(:,:,:,Kbb) = gdepw(:,:,:,Kmm)
[7646]382        !
[12779]383     ENDIF
384#endif
[2715]385      !
[7646]386      IF( ln_dynrnf .AND. ln_dynrnf_depth ) THEN       ! read depht over which runoffs are distributed
387         IF(lwp) WRITE(numout,*) 
388         IF(lwp) WRITE(numout,*) ' read in the file depht over which runoffs are distributed'
389         CALL iom_open ( "runoffs", inum )                           ! open file
390         CALL iom_get  ( inum, jpdom_data, 'rodepth', h_rnf )   ! read the river mouth array
391         CALL iom_close( inum )                                        ! close file
392         !
393         nk_rnf(:,:) = 0                               ! set the number of level over which river runoffs are applied
[12377]394         DO_2D_11_11
395            IF( h_rnf(ji,jj) > 0._wp ) THEN
396               jk = 2
397               DO WHILE ( jk /= mbkt(ji,jj) .AND. gdept_0(ji,jj,jk) < h_rnf(ji,jj) ) ;  jk = jk + 1
[7646]398               END DO
[12377]399               nk_rnf(ji,jj) = jk
400            ELSEIF( h_rnf(ji,jj) == -1._wp   ) THEN   ;  nk_rnf(ji,jj) = 1
401            ELSEIF( h_rnf(ji,jj) == -999._wp ) THEN   ;  nk_rnf(ji,jj) = mbkt(ji,jj)
402            ELSE
403               CALL ctl_stop( 'sbc_rnf_init: runoff depth not positive, and not -999 or -1, rnf value in file fort.999'  )
404               WRITE(999,*) 'ji, jj, h_rnf(ji,jj) :', ji, jj, h_rnf(ji,jj)
405            ENDIF
406         END_2D
[12779]407!!st pourquoi on n'utilise pas le gde3w ici plutôt que de faire une boucle ?
[12377]408         DO_2D_11_11
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(ji,jj,jk,Kmm)
[7646]412            END DO
[12377]413         END_2D
[7646]414      ELSE                                       ! runoffs applied at the surface
415         nk_rnf(:,:) = 1
[12377]416         h_rnf (:,:) = e3t(:,:,1,Kmm)
[7646]417      ENDIF
418      nkrnf_max = MAXVAL( nk_rnf(:,:) )
419      hrnf_max = MAXVAL( h_rnf(:,:) )
420      IF( lk_mpp )  THEN
[10425]421         CALL mpp_max( 'dtadyn', nkrnf_max )                 ! max over the  global domain
422         CALL mpp_max( 'dtadyn', hrnf_max )                 ! max over the  global domain
[7646]423      ENDIF
424      IF(lwp) WRITE(numout,*) ' '
425      IF(lwp) WRITE(numout,*) ' max depht of runoff : ', hrnf_max,'    max level  : ', nkrnf_max
426      IF(lwp) WRITE(numout,*) ' '
427      !
[12377]428      CALL dta_dyn( nit000, Kbb, Kmm, Kaa )
[2528]429      !
[1501]430   END SUBROUTINE dta_dyn_init
431
[12779]432   
[12377]433   SUBROUTINE dta_dyn_sed( kt, Kmm )
[10222]434      !!----------------------------------------------------------------------
435      !!                  ***  ROUTINE dta_dyn  ***
436      !!
437      !! ** Purpose :  Prepares dynamics and physics fields from a NEMO run
438      !!               for an off-line simulation of passive tracers
439      !!
440      !! ** Method : calculates the position of data
441      !!             - computes slopes if needed
442      !!             - interpolates data if needed
443      !!----------------------------------------------------------------------
444      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
[12377]445      INTEGER, INTENT(in) ::   Kmm  ! ocean time level index
[10222]446      !
447      !!----------------------------------------------------------------------
448      !
449      IF( ln_timing )   CALL timing_start( 'dta_dyn_sed')
450      !
451      nsecdyn = nsec_year + nsec1jan000   ! number of seconds between Jan. 1st 00h of nit000 year and the middle of time step
452      !
453      IF( kt == nit000 ) THEN    ;    nprevrec = 0
454      ELSE                       ;    nprevrec = sf_dyn(jf_tem)%nrec_a(2)
455      ENDIF
456      CALL fld_read( kt, 1, sf_dyn )      !=  read data at kt time step   ==!
457      !
[12377]458      ts(:,:,:,jp_tem,Kmm) = sf_dyn(jf_tem)%fnow(:,:,:)  * tmask(:,:,:)    ! temperature
459      ts(:,:,:,jp_sal,Kmm) = sf_dyn(jf_sal)%fnow(:,:,:)  * tmask(:,:,:)    ! salinity
[10222]460      !
[12377]461      CALL eos    ( ts(:,:,:,:,Kmm), rhd, rhop, gdept_0(:,:,:) ) ! In any case, we need rhop
[9212]462
[12377]463      IF(sn_cfctl%l_prtctl) THEN                     ! print control
464         CALL prt_ctl(tab3d_1=ts(:,:,:,jp_tem,Kmm), clinfo1=' tn      - : ', mask1=tmask,  kdim=jpk   )
465         CALL prt_ctl(tab3d_1=ts(:,:,:,jp_sal,Kmm), clinfo1=' sn      - : ', mask1=tmask,  kdim=jpk   )
[10222]466      ENDIF
467      !
468      IF( ln_timing )   CALL timing_stop( 'dta_dyn_sed')
469      !
470   END SUBROUTINE dta_dyn_sed
471
472
[12377]473   SUBROUTINE dta_dyn_sed_init( Kmm )
[10222]474      !!----------------------------------------------------------------------
475      !!                  ***  ROUTINE dta_dyn_init  ***
476      !!
477      !! ** Purpose :   Initialisation of the dynamical data
478      !! ** Method  : - read the data namdta_dyn namelist
479      !!----------------------------------------------------------------------
[12377]480      INTEGER, INTENT( in ) :: Kmm                   ! ocean time level index
481      !
[10222]482      INTEGER  :: ierr, ierr0, ierr1, ierr2, ierr3   ! return error code
483      INTEGER  :: ifpr                               ! dummy loop indice
484      INTEGER  :: jfld                               ! dummy loop arguments
485      INTEGER  :: inum, idv, idimv                   ! local integer
486      INTEGER  :: ios                                ! Local integer output status for namelist read
487      !!
488      CHARACTER(len=100)            ::  cn_dir        !   Root directory for location of core files
489      TYPE(FLD_N), DIMENSION(2) ::  slf_d         ! array of namelist informations on the fields to read
490      TYPE(FLD_N) :: sn_tem , sn_sal   !   "                 "
491      !!
492      NAMELIST/namdta_dyn/cn_dir, ln_dynrnf, ln_dynrnf_depth,  fwbcorr, &
493         &                sn_tem, sn_sal
494      !!----------------------------------------------------------------------
495      !
496      READ  ( numnam_ref, namdta_dyn, IOSTAT = ios, ERR = 901)
[11536]497901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namdta_dyn in reference namelist' )
[10222]498      READ  ( numnam_cfg, namdta_dyn, IOSTAT = ios, ERR = 902 )
[11536]499902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namdta_dyn in configuration namelist' )
[10222]500      IF(lwm) WRITE ( numond, namdta_dyn )
501      !                                         ! store namelist information in an array
502      !                                         ! Control print
503      IF(lwp) THEN
504         WRITE(numout,*)
505         WRITE(numout,*) 'dta_dyn : offline dynamics '
506         WRITE(numout,*) '~~~~~~~ '
507         WRITE(numout,*) '   Namelist namdta_dyn'
508         WRITE(numout,*) '      runoffs option enabled (T) or not (F)            ln_dynrnf        = ', ln_dynrnf
509         WRITE(numout,*) '      runoffs is spread in vertical                    ln_dynrnf_depth  = ', ln_dynrnf_depth
510         WRITE(numout,*) '      annual global mean of empmr for ssh correction   fwbcorr          = ', fwbcorr
511         WRITE(numout,*)
512      ENDIF
513      !
514      jf_tem  = 1     ;   jf_sal  = 2    ;   jfld   = jf_sal
515      !
516      slf_d(jf_tem)  = sn_tem    ;   slf_d(jf_sal)  = sn_sal
517      !
518      ALLOCATE( sf_dyn(jfld), STAT=ierr )         ! set sf structure
519      IF( ierr > 0 )  THEN
520         CALL ctl_stop( 'dta_dyn: unable to allocate sf structure' )   ;   RETURN
521      ENDIF
522      !                                         ! fill sf with slf_i and control print
523      CALL fld_fill( sf_dyn, slf_d, cn_dir, 'dta_dyn_init', 'Data in file', 'namdta_dyn' )
524      !
525      ! Open file for each variable to get his number of dimension
526      DO ifpr = 1, jfld
[12377]527         CALL fld_def( sf_dyn(ifpr) )
528         CALL iom_open( sf_dyn(ifpr)%clname, sf_dyn(ifpr)%num )
[10222]529         idv   = iom_varid( sf_dyn(ifpr)%num , slf_d(ifpr)%clvar )        ! id of the variable sdjf%clvar
530         idimv = iom_file ( sf_dyn(ifpr)%num )%ndims(idv)                 ! number of dimension for variable sdjf%clvar
[12377]531         CALL iom_close( sf_dyn(ifpr)%num )                               ! close file if already open
[10222]532         ierr1=0
533         IF( idimv == 3 ) THEN    ! 2D variable
534                                      ALLOCATE( sf_dyn(ifpr)%fnow(jpi,jpj,1)    , STAT=ierr0 )
535            IF( slf_d(ifpr)%ln_tint ) ALLOCATE( sf_dyn(ifpr)%fdta(jpi,jpj,1,2)  , STAT=ierr1 )
536         ELSE                     ! 3D variable
537                                      ALLOCATE( sf_dyn(ifpr)%fnow(jpi,jpj,jpk)  , STAT=ierr0 )
538            IF( slf_d(ifpr)%ln_tint ) ALLOCATE( sf_dyn(ifpr)%fdta(jpi,jpj,jpk,2), STAT=ierr1 )
539         ENDIF
540         IF( ierr0 + ierr1 > 0 ) THEN
541            CALL ctl_stop( 'dta_dyn_init : unable to allocate sf_dyn array structure' )   ;   RETURN
542         ENDIF
543      END DO
544      !
[12377]545      CALL dta_dyn_sed( nit000, Kmm )
[10222]546      !
547   END SUBROUTINE dta_dyn_sed_init
548
[12779]549   
[12377]550   SUBROUTINE dta_dyn_atf( kt, Kbb, Kmm, Kaa )
[7646]551     !!---------------------------------------------------------------------
552      !!                    ***  ROUTINE dta_dyn_swp  ***
553      !!
[12377]554      !! ** Purpose :   Asselin time filter of now SSH
[7646]555      !!---------------------------------------------------------------------
[12377]556      INTEGER, INTENT(in) :: kt             ! time step
557      INTEGER, INTENT(in) :: Kbb, Kmm, Kaa  ! ocean time level indices
[9212]558      !
[7646]559      !!---------------------------------------------------------------------
[6140]560
[7646]561      IF( kt == nit000 ) THEN
562         IF(lwp) WRITE(numout,*)
[12377]563         IF(lwp) WRITE(numout,*) 'dta_dyn_atf : Asselin time filter of sea surface height'
564         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ '
[7646]565      ENDIF
566
[12724]567      ssh(:,:,Kmm) = ssh(:,:,Kmm) + rn_atfp * ( ssh(:,:,Kbb) - 2 * ssh(:,:,Kmm) + ssh(:,:,Kaa)) 
[7646]568
[12377]569      !! Do we also need to time filter e3t??
570      !
571   END SUBROUTINE dta_dyn_atf
[12779]572
[12377]573   
[12779]574#if ! defined key_qco   
[12377]575   SUBROUTINE dta_dyn_sf_interp( kt, Kmm )
576      !!---------------------------------------------------------------------
577      !!                    ***  ROUTINE dta_dyn_sf_interp  ***
578      !!
579      !! ** Purpose :   Calculate scale factors at U/V/W points and depths
580      !!                given the after e3t field
581      !!---------------------------------------------------------------------
582      INTEGER, INTENT(in) :: kt   ! time step
583      INTEGER, INTENT(in) :: Kmm  ! ocean time level indices
584      !
585      INTEGER             :: ji, jj, jk
586      REAL(wp)            :: zcoef
587      !!---------------------------------------------------------------------
[7646]588
589      ! Horizontal scale factor interpolations
590      ! --------------------------------------
[12377]591      CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3u(:,:,:,Kmm), 'U' )
592      CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3v(:,:,:,Kmm), 'V' )
[7646]593
594      ! Vertical scale factor interpolations
595      ! ------------------------------------
[12377]596      CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3w (:,:,:,Kmm), 'W' )
[7646]597
598      ! t- and w- points depth
599      ! ----------------------
[12377]600      gdept(:,:,1,Kmm) = 0.5_wp * e3w(:,:,1,Kmm)
601      gdepw(:,:,1,Kmm) = 0.0_wp
[9212]602      !
[12377]603      DO_3D_11_11( 2, jpk )
604         zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk))
605         gdepw(ji,jj,jk,Kmm) = gdepw(ji,jj,jk-1,Kmm) + e3t(ji,jj,jk-1,Kmm)
606         gdept(ji,jj,jk,Kmm) =      zcoef  * ( gdepw(ji,jj,jk  ,Kmm) + 0.5 * e3w(ji,jj,jk,Kmm))  &
607            &                + (1-zcoef) * ( gdept(ji,jj,jk-1,Kmm) + e3w(ji,jj,jk,Kmm))
608      END_3D
[9212]609      !
[12377]610   END SUBROUTINE dta_dyn_sf_interp
[12779]611#endif
[7646]612
[12779]613   
[7646]614   SUBROUTINE dta_dyn_ssh( kt, phdivtr, psshb,  pemp, pssha, pe3ta )
[1501]615      !!----------------------------------------------------------------------
[7646]616      !!                ***  ROUTINE dta_dyn_wzv  ***
617      !!                   
[12377]618      !! ** Purpose :   compute the after ssh (ssh(:,:,Kaa)) and the now vertical velocity
[1501]619      !!
[7646]620      !! ** Method  : Using the incompressibility hypothesis,
621      !!        - the ssh increment is computed by integrating the horizontal divergence
622      !!          and multiply by the time step.
[1501]623      !!
[7646]624      !!        - compute the after scale factor : repartition of ssh INCREMENT proportionnaly
625      !!                                           to the level thickness ( z-star case )
626      !!
627      !!        - the vertical velocity is computed by integrating the horizontal divergence 
628      !!          from the bottom to the surface minus the scale factor evolution.
629      !!          The boundary conditions are w=0 at the bottom (no flux)
630      !!
[12779]631      !! ** action  :   ssh(:,:,Kaa) / e3t(:,:,k,Kaa) / ww
[7646]632      !!
633      !! Reference  : Leclair, M., and G. Madec, 2009, Ocean Modelling.
[2528]634      !!----------------------------------------------------------------------
[7646]635      INTEGER,                                   INTENT(in )    :: kt        !  time-step
636      REAL(wp), DIMENSION(jpi,jpj,jpk)          , INTENT(in )   :: phdivtr   ! horizontal divergence transport
637      REAL(wp), DIMENSION(jpi,jpj)    , OPTIONAL, INTENT(in )   :: psshb     ! now ssh
638      REAL(wp), DIMENSION(jpi,jpj)    , OPTIONAL, INTENT(in )   :: pemp      ! evaporation minus precipitation
639      REAL(wp), DIMENSION(jpi,jpj)    , OPTIONAL, INTENT(inout) :: pssha     ! after ssh
640      REAL(wp), DIMENSION(jpi,jpj,jpk), OPTIONAL, INTENT(out)   :: pe3ta     ! after vertical scale factor
[9212]641      !
[7646]642      INTEGER                       :: jk
643      REAL(wp), DIMENSION(jpi,jpj)  :: zhdiv 
644      REAL(wp)  :: z2dt 
645      !!----------------------------------------------------------------------
[3294]646      !
[12724]647      z2dt = 2._wp * rn_Dt
[7646]648      !
649      zhdiv(:,:) = 0._wp
650      DO jk = 1, jpkm1
651         zhdiv(:,:) = zhdiv(:,:) +  phdivtr(:,:,jk) * tmask(:,:,jk)
652      END DO
653      !                                                ! Sea surface  elevation time-stepping
[12724]654      pssha(:,:) = ( psshb(:,:) - z2dt * ( r1_rho0 * pemp(:,:)  + zhdiv(:,:) ) ) * ssmask(:,:)
[7646]655      !
[12779]656      IF( PRESENT( pe3ta ) ) THEN                      ! After acale factors at t-points ( z_star coordinate )
657         DO jk = 1, jpkm1
658            pe3ta(:,:,jk) = e3t_0(:,:,jk) * ( 1._wp + pssha(:,:) * r1_ht_0(:,:) * tmask(:,:,jk) )
659         END DO
660      ENDIF
661      !
[7646]662   END SUBROUTINE dta_dyn_ssh
663
664
[12377]665   SUBROUTINE dta_dyn_hrnf( Kmm )
[7646]666      !!----------------------------------------------------------------------
667      !!                  ***  ROUTINE sbc_rnf  ***
[1501]668      !!
[7646]669      !! ** Purpose :   update the horizontal divergence with the runoff inflow
670      !!
671      !! ** Method  :
672      !!                CAUTION : rnf is positive (inflow) decreasing the
673      !!                          divergence and expressed in m/s
674      !!
675      !! ** Action  :   phdivn   decreased by the runoff inflow
[2528]676      !!----------------------------------------------------------------------
[7646]677      !!
[12377]678      INTEGER, INTENT(in) :: Kmm  ! ocean time level index
679      !
680      INTEGER  ::   ji, jj, jk    ! dummy loop indices
[7646]681      !!----------------------------------------------------------------------
[2528]682      !
[12779]683!!st code dupliqué même remarque que plus haut pourquoi ne pas utiliser gdepw ?
[12377]684      DO_2D_11_11
685         h_rnf(ji,jj) = 0._wp
686         DO jk = 1, nk_rnf(ji,jj)                           ! recalculates h_rnf to be the depth in metres
687             h_rnf(ji,jj) = h_rnf(ji,jj) + e3t(ji,jj,jk,Kmm)   ! to the bottom of the relevant grid box
688         END DO
689      END_2D
[5836]690      !
[7646]691   END SUBROUTINE dta_dyn_hrnf
692
693
694
[12377]695   SUBROUTINE dta_dyn_slp( kt, Kbb, Kmm )
[7646]696      !!---------------------------------------------------------------------
697      !!                    ***  ROUTINE dta_dyn_slp  ***
698      !!
699      !! ** Purpose : Computation of slope
700      !!
701      !!---------------------------------------------------------------------
702      INTEGER,  INTENT(in) :: kt       ! time step
[12377]703      INTEGER,  INTENT(in) :: Kbb, Kmm ! ocean time level indices
[7646]704      !
705      INTEGER  ::   ji, jj     ! dummy loop indices
706      REAL(wp) ::   ztinta     ! ratio applied to after  records when doing time interpolation
707      REAL(wp) ::   ztintb     ! ratio applied to before records when doing time interpolation
708      INTEGER  ::   iswap 
[9212]709      REAL(wp), DIMENSION(jpi,jpj,jpk)      ::   zuslp, zvslp, zwslpi, zwslpj
710      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts) ::   zts
[7646]711      !!---------------------------------------------------------------------
712      !
713      IF( sf_dyn(jf_tem)%ln_tint ) THEN    ! Computes slopes (here avt is used as workspace)                       
714         IF( kt == nit000 ) THEN
715            IF(lwp) WRITE(numout,*) ' Compute new slopes at kt = ', kt
716            zts(:,:,:,jp_tem) = sf_dyn(jf_tem)%fdta(:,:,:,1) * tmask(:,:,:)   ! temperature
717            zts(:,:,:,jp_sal) = sf_dyn(jf_sal)%fdta(:,:,:,1) * tmask(:,:,:)   ! salinity
718            avt(:,:,:)        = sf_dyn(jf_avt)%fdta(:,:,:,1) * tmask(:,:,:)   ! vertical diffusive coef.
[12377]719            CALL compute_slopes( kt, zts, zuslp, zvslp, zwslpi, zwslpj, Kbb, Kmm )
[7646]720            uslpdta (:,:,:,1) = zuslp (:,:,:) 
721            vslpdta (:,:,:,1) = zvslp (:,:,:) 
722            wslpidta(:,:,:,1) = zwslpi(:,:,:) 
723            wslpjdta(:,:,:,1) = zwslpj(:,:,:) 
724            !
725            zts(:,:,:,jp_tem) = sf_dyn(jf_tem)%fdta(:,:,:,2) * tmask(:,:,:)   ! temperature
726            zts(:,:,:,jp_sal) = sf_dyn(jf_sal)%fdta(:,:,:,2) * tmask(:,:,:)   ! salinity
727            avt(:,:,:)        = sf_dyn(jf_avt)%fdta(:,:,:,2) * tmask(:,:,:)   ! vertical diffusive coef.
[12377]728            CALL compute_slopes( kt, zts, zuslp, zvslp, zwslpi, zwslpj, Kbb, Kmm )
[7646]729            uslpdta (:,:,:,2) = zuslp (:,:,:) 
730            vslpdta (:,:,:,2) = zvslp (:,:,:) 
731            wslpidta(:,:,:,2) = zwslpi(:,:,:) 
732            wslpjdta(:,:,:,2) = zwslpj(:,:,:) 
733         ELSE
734           !
735           iswap = 0
736           IF( sf_dyn(jf_tem)%nrec_a(2) - nprevrec /= 0 )  iswap = 1
737           IF( nsecdyn > sf_dyn(jf_tem)%nrec_b(2) .AND. iswap == 1 )  THEN    ! read/update the after data
738              IF(lwp) WRITE(numout,*) ' Compute new slopes at kt = ', kt
739              uslpdta (:,:,:,1) =  uslpdta (:,:,:,2)         ! swap the data
740              vslpdta (:,:,:,1) =  vslpdta (:,:,:,2) 
741              wslpidta(:,:,:,1) =  wslpidta(:,:,:,2) 
742              wslpjdta(:,:,:,1) =  wslpjdta(:,:,:,2) 
743              !
744              zts(:,:,:,jp_tem) = sf_dyn(jf_tem)%fdta(:,:,:,2) * tmask(:,:,:)   ! temperature
745              zts(:,:,:,jp_sal) = sf_dyn(jf_sal)%fdta(:,:,:,2) * tmask(:,:,:)   ! salinity
746              avt(:,:,:)        = sf_dyn(jf_avt)%fdta(:,:,:,2) * tmask(:,:,:)   ! vertical diffusive coef.
[12377]747              CALL compute_slopes( kt, zts, zuslp, zvslp, zwslpi, zwslpj, Kbb, Kmm )
[7646]748              !
749              uslpdta (:,:,:,2) = zuslp (:,:,:) 
750              vslpdta (:,:,:,2) = zvslp (:,:,:) 
751              wslpidta(:,:,:,2) = zwslpi(:,:,:) 
752              wslpjdta(:,:,:,2) = zwslpj(:,:,:) 
753            ENDIF
754         ENDIF
755      ENDIF
756      !
757      IF( sf_dyn(jf_tem)%ln_tint )  THEN
758         ztinta =  REAL( nsecdyn - sf_dyn(jf_tem)%nrec_b(2), wp )  &
759            &    / REAL( sf_dyn(jf_tem)%nrec_a(2) - sf_dyn(jf_tem)%nrec_b(2), wp )
760         ztintb =  1. - ztinta
761         IF( l_ldfslp .AND. .NOT.lk_c1d ) THEN    ! Computes slopes (here avt is used as workspace)
762            uslp (:,:,:) = ztintb * uslpdta (:,:,:,1)  + ztinta * uslpdta (:,:,:,2) 
763            vslp (:,:,:) = ztintb * vslpdta (:,:,:,1)  + ztinta * vslpdta (:,:,:,2) 
764            wslpi(:,:,:) = ztintb * wslpidta(:,:,:,1)  + ztinta * wslpidta(:,:,:,2) 
765            wslpj(:,:,:) = ztintb * wslpjdta(:,:,:,1)  + ztinta * wslpjdta(:,:,:,2) 
766         ENDIF
767      ELSE
768         zts(:,:,:,jp_tem) = sf_dyn(jf_tem)%fnow(:,:,:) * tmask(:,:,:)   ! temperature
769         zts(:,:,:,jp_sal) = sf_dyn(jf_sal)%fnow(:,:,:) * tmask(:,:,:)   ! salinity
770         avt(:,:,:)        = sf_dyn(jf_avt)%fnow(:,:,:) * tmask(:,:,:)   ! vertical diffusive coef.
[12377]771         CALL compute_slopes( kt, zts, zuslp, zvslp, zwslpi, zwslpj, Kbb, Kmm )
[7646]772         !
773         IF( l_ldfslp .AND. .NOT.lk_c1d ) THEN    ! Computes slopes (here avt is used as workspace)
774            uslp (:,:,:) = zuslp (:,:,:)
775            vslp (:,:,:) = zvslp (:,:,:)
776            wslpi(:,:,:) = zwslpi(:,:,:)
777            wslpj(:,:,:) = zwslpj(:,:,:)
778         ENDIF
779      ENDIF
780      !
781   END SUBROUTINE dta_dyn_slp
[1501]782
[9212]783
[12377]784   SUBROUTINE compute_slopes( kt, pts, puslp, pvslp, pwslpi, pwslpj, Kbb, Kmm )
[1501]785      !!---------------------------------------------------------------------
[3294]786      !!                    ***  ROUTINE dta_dyn_slp  ***
[1501]787      !!
[9212]788      !! ** Purpose :   Computation of slope
[1501]789      !!---------------------------------------------------------------------
[3294]790      INTEGER ,                              INTENT(in ) :: kt       ! time step
791      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in ) :: pts      ! temperature/salinity
792      REAL(wp), DIMENSION(jpi,jpj,jpk)     , INTENT(out) :: puslp    ! zonal isopycnal slopes
793      REAL(wp), DIMENSION(jpi,jpj,jpk)     , INTENT(out) :: pvslp    ! meridional isopycnal slopes
794      REAL(wp), DIMENSION(jpi,jpj,jpk)     , INTENT(out) :: pwslpi   ! zonal diapycnal slopes
795      REAL(wp), DIMENSION(jpi,jpj,jpk)     , INTENT(out) :: pwslpj   ! meridional diapycnal slopes
[12377]796      INTEGER ,                              INTENT(in ) :: Kbb, Kmm ! ocean time level indices
[1501]797      !!---------------------------------------------------------------------
[9212]798      !
[7646]799      IF( l_ldfslp .AND. .NOT.lk_c1d ) THEN    ! Computes slopes (here avt is used as workspace)
[5836]800         CALL eos    ( pts, rhd, rhop, gdept_0(:,:,:) )
[12377]801         CALL eos_rab( pts, rab_n, Kmm )       ! now local thermal/haline expension ratio at T-points
802         CALL bn2    ( pts, rab_n, rn2, Kmm  ) ! now    Brunt-Vaisala
[5131]803
[6140]804      ! Partial steps: before Horizontal DErivative
805      IF( ln_zps  .AND. .NOT. ln_isfcav)                            &
[12377]806         &            CALL zps_hde    ( kt, Kmm, jpts, pts, gtsu, gtsv,  &  ! Partial steps: before horizontal gradient
[6140]807         &                                        rhd, gru , grv    )  ! of t, s, rd at the last ocean level
808      IF( ln_zps .AND.        ln_isfcav)                            &
[12377]809         &            CALL zps_hde_isf( kt, Kmm, jpts, pts, gtsu, gtsv, gtui, gtvi, &  ! Partial steps for top cell (ISF)
[7646]810         &                                        rhd, gru , grv , grui, grvi )  ! of t, s, rd at the first ocean level
[4990]811
[12377]812         rn2b(:,:,:) = rn2(:,:,:)                ! needed for zdfmxl
813         CALL zdf_mxl( kt, Kmm )                 ! mixed layer depth
814         CALL ldf_slp( kt, rhd, rn2, Kbb, Kmm )  ! slopes
[7646]815         puslp (:,:,:) = uslp (:,:,:)
816         pvslp (:,:,:) = vslp (:,:,:)
817         pwslpi(:,:,:) = wslpi(:,:,:)
818         pwslpj(:,:,:) = wslpj(:,:,:)
[5836]819     ELSE
820         puslp (:,:,:) = 0.            ! to avoid warning when compiling
821         pvslp (:,:,:) = 0.
822         pwslpi(:,:,:) = 0.
823         pwslpj(:,:,:) = 0.
824     ENDIF
[2528]825      !
[7646]826   END SUBROUTINE compute_slopes
[9212]827
[2528]828   !!======================================================================
[325]829END MODULE dtadyn
Note: See TracBrowser for help on using the repository browser.