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/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OFF – NEMO

source: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OFF/dtadyn.F90 @ 11822

Last change on this file since 11822 was 11822, checked in by acc, 4 years ago

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