source: NEMO/trunk/src/OFF/dtadyn.F90 @ 12377

Last change on this file since 12377 was 12377, checked in by acc, 18 months ago

The big one. Merging all 2019 developments from the option 1 branch back onto the trunk.

This changeset reproduces 2019/dev_r11943_MERGE_2019 on the trunk using a 2-URL merge
onto a working copy of the trunk. I.e.:

svn merge —ignore-ancestry \

svn+ssh://acc@forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/NEMO/trunk \
svn+ssh://acc@forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/NEMO/branches/2019/dev_r11943_MERGE_2019 ./

The —ignore-ancestry flag avoids problems that may otherwise arise from the fact that
the merge history been trunk and branch may have been applied in a different order but
care has been taken before this step to ensure that all applicable fixes and updates
are present in the merge branch.

The trunk state just before this step has been branched to releases/release-4.0-HEAD
and that branch has been immediately tagged as releases/release-4.0.2. Any fixes
or additions in response to tickets on 4.0, 4.0.1 or 4.0.2 should be done on
releases/release-4.0-HEAD. From now on future 'point' releases (e.g. 4.0.2) will
remain unchanged with periodic releases as needs demand. Note release-4.0-HEAD is a
transitional naming convention. Future full releases, say 4.2, will have a release-4.2
branch which fulfills this role and the first point release (e.g. 4.2.0) will be made
immediately following the release branch creation.

2020 developments can be started from any trunk revision later than this one.

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