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

Last change on this file since 14053 was 14053, checked in by techene, 4 months ago

#2385 added to the trunk

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