source: NEMO/branches/2019/fix_sn_cfctl_ticket2328/src/OFF/dtadyn.F90 @ 11869

Last change on this file since 11869 was 11869, checked in by acc, 11 months ago

Branch 2019/fix_sn_cfctl_ticket2328. Changes to enable correct functionality for the sn_cfctl%l_mppout and sn_cfctl%l_mpptop options. These changes also introduce a sn_cfctl%l_oasout option to toggle the OASIS setup information (sbccpl.F90, only) which was yet another misuse of ln_ctl. The next step may be to remove most references to ln_ctl altogether and provide a single control mechanism. TBD. See ticket #2328

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