source: NEMO/branches/2019/dev_r11613_ENHANCE-04_namelists_as_internalfiles/src/OFF/dtadyn.F90

Last change on this file was 11671, checked in by acc, 16 months ago

Branch 2019/dev_r11613_ENHANCE-04_namelists_as_internalfiles. Final, non-substantive changes to complete this branch. These changes remove all REWIND statements on the old namelist fortran units (now character variables for internal files). These changes have been left until last since they are easily repeated via a script and it may be preferable to use the previous revision for merge purposes and reapply these last changes separately. This branch has been fully SETTE tested.

  • Property svn:keywords set to Id
File size: 40.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   !!----------------------------------------------------------------------
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) 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      READ  ( numnam_ref, namdta_dyn, IOSTAT = ios, ERR = 901)
228901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namdta_dyn in reference namelist' )
229      READ  ( numnam_cfg, namdta_dyn, IOSTAT = ios, ERR = 902 )
230902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namdta_dyn in configuration namelist' )
231      IF(lwm) WRITE ( numond, namdta_dyn )
232      !                                         ! store namelist information in an array
233      !                                         ! Control print
234      IF(lwp) THEN
235         WRITE(numout,*)
236         WRITE(numout,*) 'dta_dyn : offline dynamics '
237         WRITE(numout,*) '~~~~~~~ '
238         WRITE(numout,*) '   Namelist namdta_dyn'
239         WRITE(numout,*) '      runoffs option enabled (T) or not (F)            ln_dynrnf        = ', ln_dynrnf
240         WRITE(numout,*) '      runoffs is spread in vertical                    ln_dynrnf_depth  = ', ln_dynrnf_depth
241         WRITE(numout,*) '      annual global mean of empmr for ssh correction   fwbcorr          = ', fwbcorr
242         WRITE(numout,*)
243      ENDIF
244      !
245      jf_uwd  = 1     ;   jf_vwd  = 2    ;   jf_wwd = 3    ;   jf_emp = 4    ;   jf_avt = 5
246      jf_tem  = 6     ;   jf_sal  = 7    ;   jf_mld = 8    ;   jf_qsr = 9
247      jf_wnd  = 10    ;   jf_ice  = 11   ;   jf_fmf = 12   ;   jfld   = jf_fmf
248      !
249      slf_d(jf_uwd)  = sn_uwd    ;   slf_d(jf_vwd)  = sn_vwd   ;   slf_d(jf_wwd) = sn_wwd
250      slf_d(jf_emp)  = sn_emp    ;   slf_d(jf_avt)  = sn_avt
251      slf_d(jf_tem)  = sn_tem    ;   slf_d(jf_sal)  = sn_sal   ;   slf_d(jf_mld) = sn_mld
252      slf_d(jf_qsr)  = sn_qsr    ;   slf_d(jf_wnd)  = sn_wnd   ;   slf_d(jf_ice) = sn_ice
253      slf_d(jf_fmf)  = sn_fmf
254      !
255      IF( .NOT.ln_linssh ) THEN
256               jf_div  = jfld + 1   ;         jf_empb  = jfld + 2    ;   jfld = jf_empb
257         slf_d(jf_div) = sn_div     ;   slf_d(jf_empb) = sn_empb
258      ENDIF
259      !
260      IF( ln_trabbl ) THEN
261               jf_ubl  = jfld + 1   ;         jf_vbl  = jfld + 2     ;   jfld = jf_vbl
262         slf_d(jf_ubl) = sn_ubl     ;   slf_d(jf_vbl) = sn_vbl
263      ENDIF
264      !
265      IF( ln_dynrnf ) THEN
266               jf_rnf  = jfld + 1   ;     jfld  = jf_rnf
267         slf_d(jf_rnf) = sn_rnf   
268      ELSE
269         rnf(:,:) = 0._wp
270      ENDIF
271
272      ALLOCATE( sf_dyn(jfld), STAT=ierr )         ! set sf structure
273      IF( ierr > 0 )  THEN
274         CALL ctl_stop( 'dta_dyn: unable to allocate sf structure' )   ;   RETURN
275      ENDIF
276      !                                         ! fill sf with slf_i and control print
277      CALL fld_fill( sf_dyn, slf_d, cn_dir, 'dta_dyn_init', 'Data in file', 'namdta_dyn' )
278      !
279      ! Open file for each variable to get his number of dimension
280      DO ifpr = 1, jfld
281         CALL fld_clopn( sf_dyn(ifpr), nyear, nmonth, nday )
282         idv   = iom_varid( sf_dyn(ifpr)%num , slf_d(ifpr)%clvar )        ! id of the variable sdjf%clvar
283         idimv = iom_file ( sf_dyn(ifpr)%num )%ndims(idv)                 ! number of dimension for variable sdjf%clvar
284         IF( sf_dyn(ifpr)%num /= 0 )   CALL iom_close( sf_dyn(ifpr)%num ) ! close file if already open
285         ierr1=0
286         IF( idimv == 3 ) THEN    ! 2D variable
287                                      ALLOCATE( sf_dyn(ifpr)%fnow(jpi,jpj,1)    , STAT=ierr0 )
288            IF( slf_d(ifpr)%ln_tint ) ALLOCATE( sf_dyn(ifpr)%fdta(jpi,jpj,1,2)  , STAT=ierr1 )
289         ELSE                     ! 3D variable
290                                      ALLOCATE( sf_dyn(ifpr)%fnow(jpi,jpj,jpk)  , STAT=ierr0 )
291            IF( slf_d(ifpr)%ln_tint ) ALLOCATE( sf_dyn(ifpr)%fdta(jpi,jpj,jpk,2), STAT=ierr1 )
292         ENDIF
293         IF( ierr0 + ierr1 > 0 ) THEN
294            CALL ctl_stop( 'dta_dyn_init : unable to allocate sf_dyn array structure' )   ;   RETURN
295         ENDIF
296      END DO
297      !
298      IF( l_ldfslp .AND. .NOT.lk_c1d ) THEN                  ! slopes
299         IF( sf_dyn(jf_tem)%ln_tint ) THEN      ! time interpolation
300            ALLOCATE( uslpdta (jpi,jpj,jpk,2), vslpdta (jpi,jpj,jpk,2),    &
301            &         wslpidta(jpi,jpj,jpk,2), wslpjdta(jpi,jpj,jpk,2), STAT=ierr2 )
302            !
303            IF( ierr2 > 0 )  THEN
304               CALL ctl_stop( 'dta_dyn_init : unable to allocate slope arrays' )   ;   RETURN
305            ENDIF
306         ENDIF
307      ENDIF
308      !
309      IF( .NOT.ln_linssh ) THEN
310        IF( .NOT. sf_dyn(jf_uwd)%ln_clim .AND. ln_rsttr .AND.    &                     ! Restart: read in restart file
311           iom_varid( numrtr, 'sshn', ldstop = .FALSE. ) > 0 ) THEN
312           IF(lwp) WRITE(numout,*) ' sshn forcing fields read in the restart file for initialisation'
313           CALL iom_get( numrtr, jpdom_autoglo, 'sshn', sshn(:,:)   )
314           CALL iom_get( numrtr, jpdom_autoglo, 'sshb', sshb(:,:)   )
315        ELSE
316           IF(lwp) WRITE(numout,*) ' sshn forcing fields read in the restart file for initialisation'
317           CALL iom_open( 'restart', inum )
318           CALL iom_get( inum, jpdom_autoglo, 'sshn', sshn(:,:)   )
319           CALL iom_get( inum, jpdom_autoglo, 'sshb', sshb(:,:)   )
320           CALL iom_close( inum )                                        ! close file
321        ENDIF
322        !
323        DO jk = 1, jpkm1
324           e3t_n(:,:,jk) = e3t_0(:,:,jk) * ( 1._wp + sshn(:,:) * tmask(:,:,1) / ( ht_0(:,:) + 1.0 - tmask(:,:,1) ) )
325        ENDDO
326        e3t_a(:,:,jpk) = e3t_0(:,:,jpk)
327
328        ! Horizontal scale factor interpolations
329        ! --------------------------------------
330        CALL dom_vvl_interpol( e3t_n(:,:,:), e3u_n(:,:,:), 'U' )
331        CALL dom_vvl_interpol( e3t_n(:,:,:), e3v_n(:,:,:), 'V' )
332
333        ! Vertical scale factor interpolations
334        ! ------------------------------------
335        CALL dom_vvl_interpol( e3t_n(:,:,:), e3w_n(:,:,:), 'W' )
336 
337        e3t_b(:,:,:)  = e3t_n(:,:,:)
338        e3u_b(:,:,:)  = e3u_n(:,:,:)
339        e3v_b(:,:,:)  = e3v_n(:,:,:)
340
341        ! t- and w- points depth
342        ! ----------------------
343        gdept_n(:,:,1) = 0.5_wp * e3w_n(:,:,1)
344        gdepw_n(:,:,1) = 0.0_wp
345
346        DO jk = 2, jpk
347           DO jj = 1,jpj
348              DO ji = 1,jpi
349                !    zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk))   ! 0 everywhere
350                !    tmask = wmask, ie everywhere expect at jk = mikt
351                                                                   ! 1 for jk =
352                                                                   ! mikt
353                 zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk))
354                 gdepw_n(ji,jj,jk) = gdepw_n(ji,jj,jk-1) + e3t_n(ji,jj,jk-1)
355                 gdept_n(ji,jj,jk) =      zcoef  * ( gdepw_n(ji,jj,jk  ) + 0.5 * e3w_n(ji,jj,jk))  &
356                     &                + (1-zcoef) * ( gdept_n(ji,jj,jk-1) + e3w_n(ji,jj,jk))
357              END DO
358           END DO
359        END DO
360
361        gdept_b(:,:,:) = gdept_n(:,:,:)
362        gdepw_b(:,:,:) = gdepw_n(:,:,:)
363        !
364      ENDIF
365      !
366      IF( ln_dynrnf .AND. ln_dynrnf_depth ) THEN       ! read depht over which runoffs are distributed
367         IF(lwp) WRITE(numout,*) 
368         IF(lwp) WRITE(numout,*) ' read in the file depht over which runoffs are distributed'
369         CALL iom_open ( "runoffs", inum )                           ! open file
370         CALL iom_get  ( inum, jpdom_data, 'rodepth', h_rnf )   ! read the river mouth array
371         CALL iom_close( inum )                                        ! close file
372         !
373         nk_rnf(:,:) = 0                               ! set the number of level over which river runoffs are applied
374         DO jj = 1, jpj
375            DO ji = 1, jpi
376               IF( h_rnf(ji,jj) > 0._wp ) THEN
377                  jk = 2
378                  DO WHILE ( jk /= mbkt(ji,jj) .AND. gdept_0(ji,jj,jk) < h_rnf(ji,jj) ) ;  jk = jk + 1
379                  END DO
380                  nk_rnf(ji,jj) = jk
381               ELSEIF( h_rnf(ji,jj) == -1._wp   ) THEN   ;  nk_rnf(ji,jj) = 1
382               ELSEIF( h_rnf(ji,jj) == -999._wp ) THEN   ;  nk_rnf(ji,jj) = mbkt(ji,jj)
383               ELSE
384                  CALL ctl_stop( 'sbc_rnf_init: runoff depth not positive, and not -999 or -1, rnf value in file fort.999'  )
385                  WRITE(999,*) 'ji, jj, h_rnf(ji,jj) :', ji, jj, h_rnf(ji,jj)
386               ENDIF
387            END DO
388         END DO
389         DO jj = 1, jpj                                ! set the associated depth
390            DO ji = 1, jpi
391               h_rnf(ji,jj) = 0._wp
392               DO jk = 1, nk_rnf(ji,jj)
393                  h_rnf(ji,jj) = h_rnf(ji,jj) + e3t_n(ji,jj,jk)
394               END DO
395            END DO
396         END DO
397      ELSE                                       ! runoffs applied at the surface
398         nk_rnf(:,:) = 1
399         h_rnf (:,:) = e3t_n(:,:,1)
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 )
412      !
413   END SUBROUTINE dta_dyn_init
414
415   SUBROUTINE dta_dyn_sed( kt )
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      !
428      !!----------------------------------------------------------------------
429      !
430      IF( ln_timing )   CALL timing_start( 'dta_dyn_sed')
431      !
432      nsecdyn = nsec_year + nsec1jan000   ! number of seconds between Jan. 1st 00h of nit000 year and the middle of time step
433      !
434      IF( kt == nit000 ) THEN    ;    nprevrec = 0
435      ELSE                       ;    nprevrec = sf_dyn(jf_tem)%nrec_a(2)
436      ENDIF
437      CALL fld_read( kt, 1, sf_dyn )      !=  read data at kt time step   ==!
438      !
439      tsn(:,:,:,jp_tem) = sf_dyn(jf_tem)%fnow(:,:,:)  * tmask(:,:,:)    ! temperature
440      tsn(:,:,:,jp_sal) = sf_dyn(jf_sal)%fnow(:,:,:)  * tmask(:,:,:)    ! salinity
441      !
442      CALL eos    ( tsn, rhd, rhop, gdept_0(:,:,:) ) ! In any case, we need rhop
443
444      IF(ln_ctl) THEN                  ! print control
445         CALL prt_ctl(tab3d_1=tsn(:,:,:,jp_tem), clinfo1=' tn      - : ', mask1=tmask,  kdim=jpk   )
446         CALL prt_ctl(tab3d_1=tsn(:,:,:,jp_sal), clinfo1=' sn      - : ', mask1=tmask,  kdim=jpk   )
447      ENDIF
448      !
449      IF( ln_timing )   CALL timing_stop( 'dta_dyn_sed')
450      !
451   END SUBROUTINE dta_dyn_sed
452
453
454   SUBROUTINE dta_dyn_sed_init
455      !!----------------------------------------------------------------------
456      !!                  ***  ROUTINE dta_dyn_init  ***
457      !!
458      !! ** Purpose :   Initialisation of the dynamical data
459      !! ** Method  : - read the data namdta_dyn namelist
460      !!----------------------------------------------------------------------
461      INTEGER  :: ierr, ierr0, ierr1, ierr2, ierr3   ! return error code
462      INTEGER  :: ifpr                               ! dummy loop indice
463      INTEGER  :: jfld                               ! dummy loop arguments
464      INTEGER  :: inum, idv, idimv                   ! local integer
465      INTEGER  :: ios                                ! Local integer output status for namelist read
466      !!
467      CHARACTER(len=100)            ::  cn_dir        !   Root directory for location of core files
468      TYPE(FLD_N), DIMENSION(2) ::  slf_d         ! array of namelist informations on the fields to read
469      TYPE(FLD_N) :: sn_tem , sn_sal   !   "                 "
470      !!
471      NAMELIST/namdta_dyn/cn_dir, ln_dynrnf, ln_dynrnf_depth,  fwbcorr, &
472         &                sn_tem, sn_sal
473      !!----------------------------------------------------------------------
474      !
475      READ  ( numnam_ref, namdta_dyn, IOSTAT = ios, ERR = 901)
476901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namdta_dyn in reference namelist' )
477      READ  ( numnam_cfg, namdta_dyn, IOSTAT = ios, ERR = 902 )
478902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namdta_dyn in configuration namelist' )
479      IF(lwm) WRITE ( numond, namdta_dyn )
480      !                                         ! store namelist information in an array
481      !                                         ! Control print
482      IF(lwp) THEN
483         WRITE(numout,*)
484         WRITE(numout,*) 'dta_dyn : offline dynamics '
485         WRITE(numout,*) '~~~~~~~ '
486         WRITE(numout,*) '   Namelist namdta_dyn'
487         WRITE(numout,*) '      runoffs option enabled (T) or not (F)            ln_dynrnf        = ', ln_dynrnf
488         WRITE(numout,*) '      runoffs is spread in vertical                    ln_dynrnf_depth  = ', ln_dynrnf_depth
489         WRITE(numout,*) '      annual global mean of empmr for ssh correction   fwbcorr          = ', fwbcorr
490         WRITE(numout,*)
491      ENDIF
492      !
493      jf_tem  = 1     ;   jf_sal  = 2    ;   jfld   = jf_sal
494      !
495      slf_d(jf_tem)  = sn_tem    ;   slf_d(jf_sal)  = sn_sal
496      !
497      ALLOCATE( sf_dyn(jfld), STAT=ierr )         ! set sf structure
498      IF( ierr > 0 )  THEN
499         CALL ctl_stop( 'dta_dyn: unable to allocate sf structure' )   ;   RETURN
500      ENDIF
501      !                                         ! fill sf with slf_i and control print
502      CALL fld_fill( sf_dyn, slf_d, cn_dir, 'dta_dyn_init', 'Data in file', 'namdta_dyn' )
503      !
504      ! Open file for each variable to get his number of dimension
505      DO ifpr = 1, jfld
506         CALL fld_clopn( sf_dyn(ifpr), nyear, nmonth, nday )
507         idv   = iom_varid( sf_dyn(ifpr)%num , slf_d(ifpr)%clvar )        ! id of the variable sdjf%clvar
508         idimv = iom_file ( sf_dyn(ifpr)%num )%ndims(idv)                 ! number of dimension for variable sdjf%clvar
509         IF( sf_dyn(ifpr)%num /= 0 )   CALL iom_close( sf_dyn(ifpr)%num ) ! close file if already open
510         ierr1=0
511         IF( idimv == 3 ) THEN    ! 2D variable
512                                      ALLOCATE( sf_dyn(ifpr)%fnow(jpi,jpj,1)    , STAT=ierr0 )
513            IF( slf_d(ifpr)%ln_tint ) ALLOCATE( sf_dyn(ifpr)%fdta(jpi,jpj,1,2)  , STAT=ierr1 )
514         ELSE                     ! 3D variable
515                                      ALLOCATE( sf_dyn(ifpr)%fnow(jpi,jpj,jpk)  , STAT=ierr0 )
516            IF( slf_d(ifpr)%ln_tint ) ALLOCATE( sf_dyn(ifpr)%fdta(jpi,jpj,jpk,2), STAT=ierr1 )
517         ENDIF
518         IF( ierr0 + ierr1 > 0 ) THEN
519            CALL ctl_stop( 'dta_dyn_init : unable to allocate sf_dyn array structure' )   ;   RETURN
520         ENDIF
521      END DO
522      !
523      CALL dta_dyn_sed( nit000 )
524      !
525   END SUBROUTINE dta_dyn_sed_init
526
527   SUBROUTINE dta_dyn_swp( kt )
528     !!---------------------------------------------------------------------
529      !!                    ***  ROUTINE dta_dyn_swp  ***
530      !!
531      !! ** Purpose :   Swap and the data and compute the vertical scale factor
532      !!              at U/V/W pointand the depht
533      !!---------------------------------------------------------------------
534      INTEGER, INTENT(in) :: kt       ! time step
535      !
536      INTEGER             :: ji, jj, jk
537      REAL(wp)            :: zcoef
538      !!---------------------------------------------------------------------
539
540      IF( kt == nit000 ) THEN
541         IF(lwp) WRITE(numout,*)
542         IF(lwp) WRITE(numout,*) 'ssh_swp : Asselin time filter and swap of sea surface height'
543         IF(lwp) WRITE(numout,*) '~~~~~~~ '
544      ENDIF
545
546      sshb(:,:) = sshn(:,:) + atfp * ( sshb(:,:) - 2 * sshn(:,:) + ssha(:,:))  ! before <-- now filtered
547      sshn(:,:) = ssha(:,:)
548
549      e3t_n(:,:,:) = e3t_a(:,:,:)
550
551      ! Reconstruction of all vertical scale factors at now and before time steps
552      ! =============================================================================
553
554      ! Horizontal scale factor interpolations
555      ! --------------------------------------
556      CALL dom_vvl_interpol( e3t_n(:,:,:), e3u_n(:,:,:), 'U' )
557      CALL dom_vvl_interpol( e3t_n(:,:,:), e3v_n(:,:,:), 'V' )
558
559      ! Vertical scale factor interpolations
560      ! ------------------------------------
561      CALL dom_vvl_interpol( e3t_n(:,:,:), e3w_n (:,:,:), 'W' )
562
563      e3t_b(:,:,:)  = e3t_n(:,:,:)
564      e3u_b(:,:,:)  = e3u_n(:,:,:)
565      e3v_b(:,:,:)  = e3v_n(:,:,:)
566
567      ! t- and w- points depth
568      ! ----------------------
569      gdept_n(:,:,1) = 0.5_wp * e3w_n(:,:,1)
570      gdepw_n(:,:,1) = 0.0_wp
571      !
572      DO jk = 2, jpk
573         DO jj = 1,jpj
574            DO ji = 1,jpi
575               zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk))
576               gdepw_n(ji,jj,jk) = gdepw_n(ji,jj,jk-1) + e3t_n(ji,jj,jk-1)
577               gdept_n(ji,jj,jk) =      zcoef  * ( gdepw_n(ji,jj,jk  ) + 0.5 * e3w_n(ji,jj,jk))  &
578                  &                + (1-zcoef) * ( gdept_n(ji,jj,jk-1) + e3w_n(ji,jj,jk))
579            END DO
580         END DO
581      END DO
582      !
583      gdept_b(:,:,:) = gdept_n(:,:,:)
584      gdepw_b(:,:,:) = gdepw_n(:,:,:)
585      !
586   END SUBROUTINE dta_dyn_swp
587   
588
589   SUBROUTINE dta_dyn_ssh( kt, phdivtr, psshb,  pemp, pssha, pe3ta )
590      !!----------------------------------------------------------------------
591      !!                ***  ROUTINE dta_dyn_wzv  ***
592      !!                   
593      !! ** Purpose :   compute the after ssh (ssha) and the now vertical velocity
594      !!
595      !! ** Method  : Using the incompressibility hypothesis,
596      !!        - the ssh increment is computed by integrating the horizontal divergence
597      !!          and multiply by the time step.
598      !!
599      !!        - compute the after scale factor : repartition of ssh INCREMENT proportionnaly
600      !!                                           to the level thickness ( z-star case )
601      !!
602      !!        - the vertical velocity is computed by integrating the horizontal divergence 
603      !!          from the bottom to the surface minus the scale factor evolution.
604      !!          The boundary conditions are w=0 at the bottom (no flux)
605      !!
606      !! ** action  :   ssha / e3t_a / wn
607      !!
608      !! Reference  : Leclair, M., and G. Madec, 2009, Ocean Modelling.
609      !!----------------------------------------------------------------------
610      INTEGER,                                   INTENT(in )    :: kt        !  time-step
611      REAL(wp), DIMENSION(jpi,jpj,jpk)          , INTENT(in )   :: phdivtr   ! horizontal divergence transport
612      REAL(wp), DIMENSION(jpi,jpj)    , OPTIONAL, INTENT(in )   :: psshb     ! now ssh
613      REAL(wp), DIMENSION(jpi,jpj)    , OPTIONAL, INTENT(in )   :: pemp      ! evaporation minus precipitation
614      REAL(wp), DIMENSION(jpi,jpj)    , OPTIONAL, INTENT(inout) :: pssha     ! after ssh
615      REAL(wp), DIMENSION(jpi,jpj,jpk), OPTIONAL, INTENT(out)   :: pe3ta     ! after vertical scale factor
616      !
617      INTEGER                       :: jk
618      REAL(wp), DIMENSION(jpi,jpj)  :: zhdiv 
619      REAL(wp)  :: z2dt 
620      !!----------------------------------------------------------------------
621      !
622      z2dt = 2._wp * rdt
623      !
624      zhdiv(:,:) = 0._wp
625      DO jk = 1, jpkm1
626         zhdiv(:,:) = zhdiv(:,:) +  phdivtr(:,:,jk) * tmask(:,:,jk)
627      END DO
628      !                                                ! Sea surface  elevation time-stepping
629      pssha(:,:) = ( psshb(:,:) - z2dt * ( r1_rau0 * pemp(:,:)  + zhdiv(:,:) ) ) * ssmask(:,:)
630      !                                                 !
631      !                                                 ! After acale factors at t-points ( z_star coordinate )
632      DO jk = 1, jpkm1
633        pe3ta(:,:,jk) = e3t_0(:,:,jk) * ( 1._wp + pssha(:,:) * tmask(:,:,1) / ( ht_0(:,:) + 1.0 - tmask(:,:,1) ) )
634      END DO
635      !
636   END SUBROUTINE dta_dyn_ssh
637
638
639   SUBROUTINE dta_dyn_hrnf
640      !!----------------------------------------------------------------------
641      !!                  ***  ROUTINE sbc_rnf  ***
642      !!
643      !! ** Purpose :   update the horizontal divergence with the runoff inflow
644      !!
645      !! ** Method  :
646      !!                CAUTION : rnf is positive (inflow) decreasing the
647      !!                          divergence and expressed in m/s
648      !!
649      !! ** Action  :   phdivn   decreased by the runoff inflow
650      !!----------------------------------------------------------------------
651      !!
652      INTEGER  ::   ji, jj, jk   ! dummy loop indices
653      !!----------------------------------------------------------------------
654      !
655      DO jj = 1, jpj                   ! update the depth over which runoffs are distributed
656         DO ji = 1, jpi
657            h_rnf(ji,jj) = 0._wp
658            DO jk = 1, nk_rnf(ji,jj)                           ! recalculates h_rnf to be the depth in metres
659                h_rnf(ji,jj) = h_rnf(ji,jj) + e3t_n(ji,jj,jk)   ! to the bottom of the relevant grid box
660            END DO
661        END DO
662      END DO
663      !
664   END SUBROUTINE dta_dyn_hrnf
665
666
667
668   SUBROUTINE dta_dyn_slp( kt )
669      !!---------------------------------------------------------------------
670      !!                    ***  ROUTINE dta_dyn_slp  ***
671      !!
672      !! ** Purpose : Computation of slope
673      !!
674      !!---------------------------------------------------------------------
675      INTEGER,  INTENT(in) :: kt       ! time step
676      !
677      INTEGER  ::   ji, jj     ! dummy loop indices
678      REAL(wp) ::   ztinta     ! ratio applied to after  records when doing time interpolation
679      REAL(wp) ::   ztintb     ! ratio applied to before records when doing time interpolation
680      INTEGER  ::   iswap 
681      REAL(wp), DIMENSION(jpi,jpj,jpk)      ::   zuslp, zvslp, zwslpi, zwslpj
682      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts) ::   zts
683      !!---------------------------------------------------------------------
684      !
685      IF( sf_dyn(jf_tem)%ln_tint ) THEN    ! Computes slopes (here avt is used as workspace)                       
686         IF( kt == nit000 ) THEN
687            IF(lwp) WRITE(numout,*) ' Compute new slopes at kt = ', kt
688            zts(:,:,:,jp_tem) = sf_dyn(jf_tem)%fdta(:,:,:,1) * tmask(:,:,:)   ! temperature
689            zts(:,:,:,jp_sal) = sf_dyn(jf_sal)%fdta(:,:,:,1) * tmask(:,:,:)   ! salinity
690            avt(:,:,:)        = sf_dyn(jf_avt)%fdta(:,:,:,1) * tmask(:,:,:)   ! vertical diffusive coef.
691            CALL compute_slopes( kt, zts, zuslp, zvslp, zwslpi, zwslpj )
692            uslpdta (:,:,:,1) = zuslp (:,:,:) 
693            vslpdta (:,:,:,1) = zvslp (:,:,:) 
694            wslpidta(:,:,:,1) = zwslpi(:,:,:) 
695            wslpjdta(:,:,:,1) = zwslpj(:,:,:) 
696            !
697            zts(:,:,:,jp_tem) = sf_dyn(jf_tem)%fdta(:,:,:,2) * tmask(:,:,:)   ! temperature
698            zts(:,:,:,jp_sal) = sf_dyn(jf_sal)%fdta(:,:,:,2) * tmask(:,:,:)   ! salinity
699            avt(:,:,:)        = sf_dyn(jf_avt)%fdta(:,:,:,2) * tmask(:,:,:)   ! vertical diffusive coef.
700            CALL compute_slopes( kt, zts, zuslp, zvslp, zwslpi, zwslpj )
701            uslpdta (:,:,:,2) = zuslp (:,:,:) 
702            vslpdta (:,:,:,2) = zvslp (:,:,:) 
703            wslpidta(:,:,:,2) = zwslpi(:,:,:) 
704            wslpjdta(:,:,:,2) = zwslpj(:,:,:) 
705         ELSE
706           !
707           iswap = 0
708           IF( sf_dyn(jf_tem)%nrec_a(2) - nprevrec /= 0 )  iswap = 1
709           IF( nsecdyn > sf_dyn(jf_tem)%nrec_b(2) .AND. iswap == 1 )  THEN    ! read/update the after data
710              IF(lwp) WRITE(numout,*) ' Compute new slopes at kt = ', kt
711              uslpdta (:,:,:,1) =  uslpdta (:,:,:,2)         ! swap the data
712              vslpdta (:,:,:,1) =  vslpdta (:,:,:,2) 
713              wslpidta(:,:,:,1) =  wslpidta(:,:,:,2) 
714              wslpjdta(:,:,:,1) =  wslpjdta(:,:,:,2) 
715              !
716              zts(:,:,:,jp_tem) = sf_dyn(jf_tem)%fdta(:,:,:,2) * tmask(:,:,:)   ! temperature
717              zts(:,:,:,jp_sal) = sf_dyn(jf_sal)%fdta(:,:,:,2) * tmask(:,:,:)   ! salinity
718              avt(:,:,:)        = sf_dyn(jf_avt)%fdta(:,:,:,2) * tmask(:,:,:)   ! vertical diffusive coef.
719              CALL compute_slopes( kt, zts, zuslp, zvslp, zwslpi, zwslpj )
720              !
721              uslpdta (:,:,:,2) = zuslp (:,:,:) 
722              vslpdta (:,:,:,2) = zvslp (:,:,:) 
723              wslpidta(:,:,:,2) = zwslpi(:,:,:) 
724              wslpjdta(:,:,:,2) = zwslpj(:,:,:) 
725            ENDIF
726         ENDIF
727      ENDIF
728      !
729      IF( sf_dyn(jf_tem)%ln_tint )  THEN
730         ztinta =  REAL( nsecdyn - sf_dyn(jf_tem)%nrec_b(2), wp )  &
731            &    / REAL( sf_dyn(jf_tem)%nrec_a(2) - sf_dyn(jf_tem)%nrec_b(2), wp )
732         ztintb =  1. - ztinta
733         IF( l_ldfslp .AND. .NOT.lk_c1d ) THEN    ! Computes slopes (here avt is used as workspace)
734            uslp (:,:,:) = ztintb * uslpdta (:,:,:,1)  + ztinta * uslpdta (:,:,:,2) 
735            vslp (:,:,:) = ztintb * vslpdta (:,:,:,1)  + ztinta * vslpdta (:,:,:,2) 
736            wslpi(:,:,:) = ztintb * wslpidta(:,:,:,1)  + ztinta * wslpidta(:,:,:,2) 
737            wslpj(:,:,:) = ztintb * wslpjdta(:,:,:,1)  + ztinta * wslpjdta(:,:,:,2) 
738         ENDIF
739      ELSE
740         zts(:,:,:,jp_tem) = sf_dyn(jf_tem)%fnow(:,:,:) * tmask(:,:,:)   ! temperature
741         zts(:,:,:,jp_sal) = sf_dyn(jf_sal)%fnow(:,:,:) * tmask(:,:,:)   ! salinity
742         avt(:,:,:)        = sf_dyn(jf_avt)%fnow(:,:,:) * tmask(:,:,:)   ! vertical diffusive coef.
743         CALL compute_slopes( kt, zts, zuslp, zvslp, zwslpi, zwslpj )
744         !
745         IF( l_ldfslp .AND. .NOT.lk_c1d ) THEN    ! Computes slopes (here avt is used as workspace)
746            uslp (:,:,:) = zuslp (:,:,:)
747            vslp (:,:,:) = zvslp (:,:,:)
748            wslpi(:,:,:) = zwslpi(:,:,:)
749            wslpj(:,:,:) = zwslpj(:,:,:)
750         ENDIF
751      ENDIF
752      !
753   END SUBROUTINE dta_dyn_slp
754
755
756   SUBROUTINE compute_slopes( kt, pts, puslp, pvslp, pwslpi, pwslpj )
757      !!---------------------------------------------------------------------
758      !!                    ***  ROUTINE dta_dyn_slp  ***
759      !!
760      !! ** Purpose :   Computation of slope
761      !!---------------------------------------------------------------------
762      INTEGER ,                              INTENT(in ) :: kt       ! time step
763      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in ) :: pts      ! temperature/salinity
764      REAL(wp), DIMENSION(jpi,jpj,jpk)     , INTENT(out) :: puslp    ! zonal isopycnal slopes
765      REAL(wp), DIMENSION(jpi,jpj,jpk)     , INTENT(out) :: pvslp    ! meridional isopycnal slopes
766      REAL(wp), DIMENSION(jpi,jpj,jpk)     , INTENT(out) :: pwslpi   ! zonal diapycnal slopes
767      REAL(wp), DIMENSION(jpi,jpj,jpk)     , INTENT(out) :: pwslpj   ! meridional diapycnal slopes
768      !!---------------------------------------------------------------------
769      !
770      IF( l_ldfslp .AND. .NOT.lk_c1d ) THEN    ! Computes slopes (here avt is used as workspace)
771         CALL eos    ( pts, rhd, rhop, gdept_0(:,:,:) )
772         CALL eos_rab( pts, rab_n )       ! now local thermal/haline expension ratio at T-points
773         CALL bn2    ( pts, rab_n, rn2  ) ! now    Brunt-Vaisala
774
775      ! Partial steps: before Horizontal DErivative
776      IF( ln_zps  .AND. .NOT. ln_isfcav)                            &
777         &            CALL zps_hde    ( kt, jpts, pts, gtsu, gtsv,  &  ! Partial steps: before horizontal gradient
778         &                                        rhd, gru , grv    )  ! of t, s, rd at the last ocean level
779      IF( ln_zps .AND.        ln_isfcav)                            &
780         &            CALL zps_hde_isf( kt, jpts, pts, gtsu, gtsv, gtui, gtvi, &  ! Partial steps for top cell (ISF)
781         &                                        rhd, gru , grv , grui, grvi )  ! of t, s, rd at the first ocean level
782
783         rn2b(:,:,:) = rn2(:,:,:)         ! need for zdfmxl
784         CALL zdf_mxl( kt )            ! mixed layer depth
785         CALL ldf_slp( kt, rhd, rn2 )  ! slopes
786         puslp (:,:,:) = uslp (:,:,:)
787         pvslp (:,:,:) = vslp (:,:,:)
788         pwslpi(:,:,:) = wslpi(:,:,:)
789         pwslpj(:,:,:) = wslpj(:,:,:)
790     ELSE
791         puslp (:,:,:) = 0.            ! to avoid warning when compiling
792         pvslp (:,:,:) = 0.
793         pwslpi(:,:,:) = 0.
794         pwslpj(:,:,:) = 0.
795     ENDIF
796      !
797   END SUBROUTINE compute_slopes
798
799   !!======================================================================
800END MODULE dtadyn
Note: See TracBrowser for help on using the repository browser.