New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
dtadyn.F90 in NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OFF – NEMO

source: NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OFF/dtadyn.F90 @ 10009

Last change on this file since 10009 was 10009, checked in by gm, 6 years ago

#1911 (ENHANCE-04): RK3 branch - step II.1 time-level dimension on ssh
  • Property svn:keywords set to Id
File size: 36.8 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_swp   ! called by step.F90
51
52   CHARACTER(len=100) ::   cn_dir          !: Root directory for location of ssr files
53   LOGICAL            ::   ln_dynrnf       !: read runoff data in file (T) or set to zero (F)
54   LOGICAL            ::   ln_dynrnf_depth       !: read runoff data in file (T) or set to zero (F)
55   REAL(wp)           ::   fwbcorr
56
57
58   INTEGER  , PARAMETER ::   jpfld = 20     ! maximum number of fields to read
59   INTEGER  , SAVE      ::   jf_tem         ! index of temperature
60   INTEGER  , SAVE      ::   jf_sal         ! index of salinity
61   INTEGER  , SAVE      ::   jf_uwd         ! index of u-transport
62   INTEGER  , SAVE      ::   jf_vwd         ! index of v-transport
63   INTEGER  , SAVE      ::   jf_wwd         ! index of v-transport
64   INTEGER  , SAVE      ::   jf_avt         ! index of Kz
65   INTEGER  , SAVE      ::   jf_mld         ! index of mixed layer deptht
66   INTEGER  , SAVE      ::   jf_emp         ! index of water flux
67   INTEGER  , SAVE      ::   jf_empb        ! index of water flux
68   INTEGER  , SAVE      ::   jf_qsr         ! index of solar radiation
69   INTEGER  , SAVE      ::   jf_wnd         ! index of wind speed
70   INTEGER  , SAVE      ::   jf_ice         ! index of sea ice cover
71   INTEGER  , SAVE      ::   jf_rnf         ! index of river runoff
72   INTEGER  , SAVE      ::   jf_fmf         ! index of downward salt flux
73   INTEGER  , SAVE      ::   jf_ubl         ! index of u-bbl coef
74   INTEGER  , SAVE      ::   jf_vbl         ! index of v-bbl coef
75   INTEGER  , SAVE      ::   jf_div         ! index of e3t
76
77
78   TYPE(FLD), ALLOCATABLE, SAVE, DIMENSION(:) :: sf_dyn  ! structure of input fields (file informations, fields read)
79   !                                               !
80   REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: uslpdta    ! zonal isopycnal slopes
81   REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: vslpdta    ! meridional isopycnal slopes
82   REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: wslpidta   ! zonal diapycnal slopes
83   REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: wslpjdta   ! meridional diapycnal slopes
84
85   INTEGER, SAVE  :: nprevrec, nsecdyn
86
87   !!----------------------------------------------------------------------
88   !! NEMO/OFF 4.0 , NEMO Consortium (2018)
89   !! $Id$
90   !! Software governed by the CeCILL licence     (./LICENSE)
91   !!----------------------------------------------------------------------
92CONTAINS
93
94   SUBROUTINE dta_dyn( kt )
95      !!----------------------------------------------------------------------
96      !!                  ***  ROUTINE dta_dyn  ***
97      !!
98      !! ** Purpose :  Prepares dynamics and physics fields from a NEMO run
99      !!               for an off-line simulation of passive tracers
100      !!
101      !! ** Method : calculates the position of data
102      !!             - computes slopes if needed
103      !!             - interpolates data if needed
104      !!----------------------------------------------------------------------
105      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
106      !
107      INTEGER             ::   ji, jj, jk
108      REAL(wp), ALLOCATABLE, DIMENSION(:,:)   ::   zemp
109      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   zhdivtr
110      !!----------------------------------------------------------------------
111      !
112      IF( ln_timing )   CALL timing_start( 'dta_dyn')
113      !
114      nsecdyn = nsec_year + nsec1jan000   ! number of seconds between Jan. 1st 00h of nit000 year and the middle of time step
115      !
116      IF( kt == nit000 ) THEN    ;    nprevrec = 0
117      ELSE                       ;    nprevrec = sf_dyn(jf_tem)%nrec_a(2)
118      ENDIF
119      CALL fld_read( kt, 1, sf_dyn )      !=  read data at kt time step   ==!
120      !
121      IF( l_ldfslp .AND. .NOT.lk_c1d )   CALL  dta_dyn_slp( kt )    ! Computation of slopes
122      !
123      tsn(:,:,:,jp_tem) = sf_dyn(jf_tem)%fnow(:,:,:)  * tmask(:,:,:)    ! temperature
124      tsn(:,:,:,jp_sal) = sf_dyn(jf_sal)%fnow(:,:,:)  * tmask(:,:,:)    ! salinity
125      wndm(:,:)         = sf_dyn(jf_wnd)%fnow(:,:,1)  * tmask(:,:,1)    ! wind speed - needed for gas exchange
126      fmmflx(:,:)       = sf_dyn(jf_fmf)%fnow(:,:,1)  * tmask(:,:,1)    ! downward salt flux (v3.5+)
127      fr_i(:,:)         = sf_dyn(jf_ice)%fnow(:,:,1)  * tmask(:,:,1)    ! Sea-ice fraction
128      qsr (:,:)         = sf_dyn(jf_qsr)%fnow(:,:,1)  * tmask(:,:,1)    ! solar radiation
129      emp (:,:)         = sf_dyn(jf_emp)%fnow(:,:,1)  * tmask(:,:,1)    ! E-P
130      IF( ln_dynrnf ) THEN
131         rnf (:,:)      = sf_dyn(jf_rnf)%fnow(:,:,1) * tmask(:,:,1)    ! E-P
132         IF( ln_dynrnf_depth .AND. .NOT. ln_linssh )    CALL  dta_dyn_hrnf
133      ENDIF
134      !
135      un(:,:,:)        = sf_dyn(jf_uwd)%fnow(:,:,:) * umask(:,:,:)    ! effective u-transport
136      vn(:,:,:)        = sf_dyn(jf_vwd)%fnow(:,:,:) * vmask(:,:,:)    ! effective v-transport
137      wn(:,:,:)        = sf_dyn(jf_wwd)%fnow(:,:,:) * tmask(:,:,:)    ! effective v-transport
138      !
139      IF( .NOT.ln_linssh ) THEN
140         ALLOCATE( zemp(jpi,jpj) , zhdivtr(jpi,jpj,jpk) )
141         zhdivtr(:,:,:) = sf_dyn(jf_div)%fnow(:,:,:)  * tmask(:,:,:)    ! effective u-transport
142         emp_b  (:,:)   = sf_dyn(jf_empb)%fnow(:,:,1) * tmask(:,:,1)    ! E-P
143         zemp   (:,:)   = ( 0.5_wp * ( emp(:,:) + emp_b(:,:) ) + rnf(:,:) + fwbcorr ) * tmask(:,:,1)
144         CALL dta_dyn_ssh( kt, zhdivtr, ssh(:,:,Nbb), zemp, ssh(:,:,Naa) )  !=  ssh, vertical scale factor & vertical transport
145!!
146!!gm BUG ?  ssh after computed but no swap so, not used in the restart....
147!!
148         DEALLOCATE( zemp , zhdivtr )
149         !                                           Write in the tracer restart file
150         !                                          *********************************
151         IF( lrst_trc ) THEN
152            IF(lwp) WRITE(numout,*)
153            IF(lwp) WRITE(numout,*) 'dta_dyn : ssh field written in tracer restart file at it= ', kt,' date= ', ndastp
154            IF(lwp) WRITE(numout,*) '~~~~~~~'
155            CALL iom_rstput( kt, nitrst, numrtw, 'sshn', ssh(:,:,Nnn) )
156            CALL iom_rstput( kt, nitrst, numrtw, 'sshb', ssh(:,:,Nbb) )
157         ENDIF
158      ENDIF
159      !
160      CALL eos    ( tsn, rhd, rhop, gdept_0(:,:,:) ) ! In any case, we need rhop
161      CALL eos_rab( tsn, rab_n )       ! now    local thermal/haline expension ratio at T-points
162      CALL bn2    ( tsn, rab_n, rn2 ) ! before Brunt-Vaisala frequency need for zdfmxl
163
164      rn2b(:,:,:) = rn2(:,:,:)         ! need for zdfmxl
165      CALL zdf_mxl( kt )                                                   ! In any case, we need mxl
166      !
167      hmld(:,:)       = sf_dyn(jf_mld)%fnow(:,:,1) * tmask(:,:,1)    ! mixed layer depht
168      avt(:,:,:)      = sf_dyn(jf_avt)%fnow(:,:,:) * tmask(:,:,:)    ! vertical diffusive coefficient
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!         CALL prt_ctl(tab2d_1=fr_i             , clinfo1=' fr_i    - : ', mask1=tmask )
188!         CALL prt_ctl(tab2d_1=hmld             , clinfo1=' hmld    - : ', mask1=tmask )
189!         CALL prt_ctl(tab2d_1=fmmflx           , clinfo1=' fmmflx  - : ', mask1=tmask )
190!         CALL prt_ctl(tab2d_1=emp              , clinfo1=' emp     - : ', mask1=tmask )
191!         CALL prt_ctl(tab2d_1=wndm             , clinfo1=' wspd    - : ', mask1=tmask )
192!         CALL prt_ctl(tab2d_1=qsr              , clinfo1=' qsr     - : ', mask1=tmask )
193      ENDIF
194      !
195      IF( ln_timing )   CALL timing_stop( 'dta_dyn')
196      !
197   END SUBROUTINE dta_dyn
198
199
200   SUBROUTINE dta_dyn_init
201      !!----------------------------------------------------------------------
202      !!                  ***  ROUTINE dta_dyn_init  ***
203      !!
204      !! ** Purpose :   Initialisation of the dynamical data     
205      !! ** Method  : - read the data namdta_dyn namelist
206      !!----------------------------------------------------------------------
207      INTEGER  :: ierr, ierr0, ierr1, ierr2, ierr3   ! return error code
208      INTEGER  :: ifpr                               ! dummy loop indice
209      INTEGER  :: jfld                               ! dummy loop arguments
210      INTEGER  :: inum, idv, idimv                   ! local integer
211      INTEGER  :: ios                                ! Local integer output status for namelist read
212      INTEGER  :: ji, jj, jk
213      REAL(wp) :: zcoef
214      INTEGER  :: nkrnf_max
215      REAL(wp) :: hrnf_max
216      !!
217      CHARACTER(len=100)            ::  cn_dir        !   Root directory for location of core files
218      TYPE(FLD_N), DIMENSION(jpfld) ::  slf_d         ! array of namelist informations on the fields to read
219      TYPE(FLD_N) :: sn_uwd, sn_vwd, sn_wwd, sn_empb, sn_emp  ! informations about the fields to be read
220      TYPE(FLD_N) :: sn_tem , sn_sal , sn_avt   !   "                 "
221      TYPE(FLD_N) :: sn_mld, sn_qsr, sn_wnd , sn_ice , sn_fmf   !   "               "
222      TYPE(FLD_N) :: sn_ubl, sn_vbl, sn_rnf    !   "              "
223      TYPE(FLD_N) :: sn_div  ! informations about the fields to be read
224      !!
225      NAMELIST/namdta_dyn/cn_dir, ln_dynrnf, ln_dynrnf_depth,  fwbcorr, &
226         &                sn_uwd, sn_vwd, sn_wwd, sn_emp,               &
227         &                sn_avt, sn_tem, sn_sal, sn_mld , sn_qsr ,     &
228         &                sn_wnd, sn_ice, sn_fmf,                       &
229         &                sn_ubl, sn_vbl, sn_rnf,                       &
230         &                sn_empb, sn_div 
231      !!----------------------------------------------------------------------
232      !
233      REWIND( numnam_ref )              ! Namelist namdta_dyn in reference namelist : Offline: init. of dynamical data
234      READ  ( numnam_ref, namdta_dyn, IOSTAT = ios, ERR = 901)
235901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namdta_dyn in reference namelist', lwp )
236      REWIND( numnam_cfg )              ! Namelist namdta_dyn in configuration namelist : Offline: init. of dynamical data
237      READ  ( numnam_cfg, namdta_dyn, IOSTAT = ios, ERR = 902 )
238902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namdta_dyn in configuration namelist', lwp )
239      IF(lwm) WRITE ( numond, namdta_dyn )
240      !                                         ! store namelist information in an array
241      !                                         ! Control print
242      IF(lwp) THEN
243         WRITE(numout,*)
244         WRITE(numout,*) 'dta_dyn : offline dynamics '
245         WRITE(numout,*) '~~~~~~~ '
246         WRITE(numout,*) '   Namelist namdta_dyn'
247         WRITE(numout,*) '      runoffs option enabled (T) or not (F)            ln_dynrnf        = ', ln_dynrnf
248         WRITE(numout,*) '      runoffs is spread in vertical                    ln_dynrnf_depth  = ', ln_dynrnf_depth
249         WRITE(numout,*) '      annual global mean of empmr for ssh correction   fwbcorr          = ', fwbcorr
250         WRITE(numout,*)
251      ENDIF
252      !
253      jf_uwd  = 1     ;   jf_vwd  = 2    ;   jf_wwd = 3    ;   jf_emp = 4    ;   jf_avt = 5
254      jf_tem  = 6     ;   jf_sal  = 7    ;   jf_mld = 8    ;   jf_qsr = 9
255      jf_wnd  = 10    ;   jf_ice  = 11   ;   jf_fmf = 12   ;   jfld   = jf_fmf
256      !
257      slf_d(jf_uwd)  = sn_uwd    ;   slf_d(jf_vwd)  = sn_vwd   ;   slf_d(jf_wwd) = sn_wwd
258      slf_d(jf_emp)  = sn_emp    ;   slf_d(jf_avt)  = sn_avt
259      slf_d(jf_tem)  = sn_tem    ;   slf_d(jf_sal)  = sn_sal   ;   slf_d(jf_mld) = sn_mld
260      slf_d(jf_qsr)  = sn_qsr    ;   slf_d(jf_wnd)  = sn_wnd   ;   slf_d(jf_ice) = sn_ice
261      slf_d(jf_fmf)  = sn_fmf
262      !
263      IF( .NOT.ln_linssh ) THEN
264               jf_div  = jfld + 1   ;         jf_empb  = jfld + 2    ;   jfld = jf_empb
265         slf_d(jf_div) = sn_div     ;   slf_d(jf_empb) = sn_empb
266      ENDIF
267      !
268      IF( ln_trabbl ) THEN
269               jf_ubl  = jfld + 1   ;         jf_vbl  = jfld + 2     ;   jfld = jf_vbl
270         slf_d(jf_ubl) = sn_ubl     ;   slf_d(jf_vbl) = sn_vbl
271      ENDIF
272      !
273      IF( ln_dynrnf ) THEN
274               jf_rnf  = jfld + 1   ;     jfld  = jf_rnf
275         slf_d(jf_rnf) = sn_rnf   
276      ELSE
277         rnf(:,:) = 0._wp
278      ENDIF
279
280      ALLOCATE( sf_dyn(jfld), STAT=ierr )         ! set sf structure
281      IF( ierr > 0 )  THEN
282         CALL ctl_stop( 'dta_dyn: unable to allocate sf structure' )   ;   RETURN
283      ENDIF
284      !                                         ! fill sf with slf_i and control print
285      CALL fld_fill( sf_dyn, slf_d, cn_dir, 'dta_dyn_init', 'Data in file', 'namdta_dyn' )
286      !
287      ! Open file for each variable to get his number of dimension
288      DO ifpr = 1, jfld
289         CALL fld_clopn( sf_dyn(ifpr), nyear, nmonth, nday )
290         idv   = iom_varid( sf_dyn(ifpr)%num , slf_d(ifpr)%clvar )        ! id of the variable sdjf%clvar
291         idimv = iom_file ( sf_dyn(ifpr)%num )%ndims(idv)                 ! number of dimension for variable sdjf%clvar
292         IF( sf_dyn(ifpr)%num /= 0 )   CALL iom_close( sf_dyn(ifpr)%num ) ! close file if already open
293         ierr1=0
294         IF( idimv == 3 ) THEN    ! 2D variable
295                                      ALLOCATE( sf_dyn(ifpr)%fnow(jpi,jpj,1)    , STAT=ierr0 )
296            IF( slf_d(ifpr)%ln_tint ) ALLOCATE( sf_dyn(ifpr)%fdta(jpi,jpj,1,2)  , STAT=ierr1 )
297         ELSE                     ! 3D variable
298                                      ALLOCATE( sf_dyn(ifpr)%fnow(jpi,jpj,jpk)  , STAT=ierr0 )
299            IF( slf_d(ifpr)%ln_tint ) ALLOCATE( sf_dyn(ifpr)%fdta(jpi,jpj,jpk,2), STAT=ierr1 )
300         ENDIF
301         IF( ierr0 + ierr1 > 0 ) THEN
302            CALL ctl_stop( 'dta_dyn_init : unable to allocate sf_dyn array structure' )   ;   RETURN
303         ENDIF
304      END DO
305      !
306      IF( l_ldfslp .AND. .NOT.lk_c1d ) THEN                  ! slopes
307         IF( sf_dyn(jf_tem)%ln_tint ) THEN      ! time interpolation
308            ALLOCATE( uslpdta (jpi,jpj,jpk,2), vslpdta (jpi,jpj,jpk,2),    &
309            &         wslpidta(jpi,jpj,jpk,2), wslpjdta(jpi,jpj,jpk,2), STAT=ierr2 )
310            !
311            IF( ierr2 > 0 )  THEN
312               CALL ctl_stop( 'dta_dyn_init : unable to allocate slope arrays' )   ;   RETURN
313            ENDIF
314         ENDIF
315      ENDIF
316      !
317      IF( .NOT.ln_linssh ) THEN
318         IF( .NOT. sf_dyn(jf_uwd)%ln_clim .AND. ln_rsttr .AND.    &                     ! Restart: read in restart file
319            iom_varid( numrtr, 'sshn', ldstop = .FALSE. ) > 0 ) THEN
320            IF(lwp) WRITE(numout,*) ' ssh forcing fields read in the restart file for initialisation'
321            CALL iom_get( numrtr, jpdom_autoglo, 'sshn', ssh(:,:,Nnn)   )
322            CALL iom_get( numrtr, jpdom_autoglo, 'sshb', ssh(:,:,Nbb)   )
323         ELSE
324            IF(lwp) WRITE(numout,*) ' ssh forcing fields read in the restart file for initialisation'
325            CALL iom_open( 'restart', inum )
326            CALL iom_get( inum, jpdom_autoglo, 'sshn', ssh(:,:,Nnn)   )
327            CALL iom_get( inum, jpdom_autoglo, 'sshb', ssh(:,:,Nbb)   )
328            CALL iom_close( inum )                                        ! close file
329         ENDIF
330         !
331         !                    !== Set of all other vertical mesh fields  ==!  (now and before)     
332         !
333         !                          !* BEFORE fields :
334         CALL ssh2e3_before               ! set:      hu , hv , r1_hu, r1_hv
335         !                                !      e3t, e3w, e3u, e3uw, e3v, e3vw        (from 1 to jpkm1)
336         !
337         !                                ! set jpk level one to the e3._0 values
338         e3t_b(:,:,jpk) = e3t_0(:,:,jpk)  ;   e3u_b(:,:,jpk) =  e3w_0(:,:,jpk)  ;   e3v_b(:,:,jpk) =  e3v_0(:,:,jpk)
339         e3w_b(:,:,jpk) = e3w_0(:,:,jpk)  ;  e3uw_b(:,:,jpk) = e3uw_0(:,:,jpk)  ;  e3vw_b(:,:,jpk) = e3vw_0(:,:,jpk)
340         !
341         !                          !* NOW fields :
342         CALL ssh2e3_now                  ! set: ht , hu , hv , r1_hu, r1_hv
343         !                                !      e3t, e3w, e3u, e3uw, e3v, e3vw, e3f   (from 1 to jpkm1)
344         !                                !      gdept_n, gdepw_n, gde3w_n
345!!gm issue?   gdept_n, gdepw_n, gde3w_n never defined at jpk
346         !
347         !                                ! set one for all last level to the e3._0 value
348         e3t_n(:,:,jpk) = e3t_0(:,:,jpk)  ;   e3u_n(:,:,jpk) =  e3w_0(:,:,jpk)  ;   e3v_n(:,:,jpk) =  e3v_0(:,:,jpk)
349         e3w_n(:,:,jpk) = e3w_0(:,:,jpk)  ;  e3uw_n(:,:,jpk) = e3uw_0(:,:,jpk)  ;  e3vw_n(:,:,jpk) = e3vw_0(:,:,jpk)
350         e3f_n(:,:,jpk) = e3f_0(:,:,jpk)
351         !
352         !                          !* AFTER fields : (last level for OPA, 3D required for AGRIF initialisation)
353         e3t_a(:,:,:) = e3t_n(:,:,:)   ;   e3u_a(:,:,:) = e3u_n(:,:,:)   ;   e3v_a(:,:,:) = e3v_n(:,:,:)
354         !
355      ENDIF
356      !
357      IF( ln_dynrnf .AND. ln_dynrnf_depth ) THEN       ! read depht over which runoffs are distributed
358         IF(lwp) WRITE(numout,*) 
359         IF(lwp) WRITE(numout,*) ' read in the file depht over which runoffs are distributed'
360         CALL iom_open ( "runoffs", inum )                           ! open file
361         CALL iom_get  ( inum, jpdom_data, 'rodepth', h_rnf )   ! read the river mouth array
362         CALL iom_close( inum )                                        ! close file
363         !
364         nk_rnf(:,:) = 0                               ! set the number of level over which river runoffs are applied
365         DO jj = 1, jpj
366            DO ji = 1, jpi
367               IF( h_rnf(ji,jj) > 0._wp ) THEN
368                  jk = 2
369                  DO WHILE ( jk /= mbkt(ji,jj) .AND. gdept_0(ji,jj,jk) < h_rnf(ji,jj) ) ;  jk = jk + 1
370                  END DO
371                  nk_rnf(ji,jj) = jk
372               ELSEIF( h_rnf(ji,jj) == -1._wp   ) THEN   ;  nk_rnf(ji,jj) = 1
373               ELSEIF( h_rnf(ji,jj) == -999._wp ) THEN   ;  nk_rnf(ji,jj) = mbkt(ji,jj)
374               ELSE
375                  CALL ctl_stop( 'sbc_rnf_init: runoff depth not positive, and not -999 or -1, rnf value in file fort.999'  )
376                  WRITE(999,*) 'ji, jj, h_rnf(ji,jj) :', ji, jj, h_rnf(ji,jj)
377               ENDIF
378            END DO
379         END DO
380         DO jj = 1, jpj                                ! set the associated depth
381            DO ji = 1, jpi
382               h_rnf(ji,jj) = 0._wp
383               DO jk = 1, nk_rnf(ji,jj)
384                  h_rnf(ji,jj) = h_rnf(ji,jj) + e3t_n(ji,jj,jk)
385               END DO
386            END DO
387         END DO
388      ELSE                                       ! runoffs applied at the surface
389         nk_rnf(:,:) = 1
390         h_rnf (:,:) = e3t_n(:,:,1)
391      ENDIF
392      nkrnf_max = MAXVAL( nk_rnf(:,:) )
393      hrnf_max = MAXVAL( h_rnf(:,:) )
394      IF( lk_mpp )  THEN
395         CALL mpp_max( nkrnf_max )                 ! max over the  global domain
396         CALL mpp_max( hrnf_max )                 ! max over the  global domain
397      ENDIF
398      IF(lwp) WRITE(numout,*) ' '
399      IF(lwp) WRITE(numout,*) ' max depht of runoff : ', hrnf_max,'    max level  : ', nkrnf_max
400      IF(lwp) WRITE(numout,*) ' '
401      !
402      CALL dta_dyn( nit000 )
403      !
404   END SUBROUTINE dta_dyn_init
405
406
407   SUBROUTINE dta_dyn_swp( kt )
408     !!---------------------------------------------------------------------
409      !!                    ***  ROUTINE dta_dyn_swp  ***
410      !!
411      !! ** Purpose :   Swap and the data and compute the vertical scale factor
412      !!              at U/V/W pointand the depht
413      !!---------------------------------------------------------------------
414      INTEGER, INTENT(in) :: kt       ! time step
415      !
416      INTEGER             :: ji, jj, jk
417      REAL(wp)            :: zcoef
418      REAL(wp), DIMENSION(jpi,jpj) ::   zssht_h, zsshu_h, zsshv_h
419      !!---------------------------------------------------------------------
420
421      IF( kt == nit000 ) THEN
422         IF(lwp) WRITE(numout,*)
423         IF(lwp) WRITE(numout,*) 'ssh_swp : Asselin time filter and swap of sea surface height'
424         IF(lwp) WRITE(numout,*) '~~~~~~~ '
425      ENDIF
426
427      ssh(:,:,Nbb) = ssh(:,:,Nnn) + rn_atfp * ( ssh(:,:,Nbb) - 2 * ssh(:,:,Nnn) + ssh(:,:,Naa) )   ! before <-- now filtered
428      ssh(:,:,Nnn) = ssh(:,:,Naa)
429
430      ! Reconstruction of all vertical scale factors at now and before time steps
431      ! =============================================================================
432      !
433      !                                   !==  now ssh  ==!  (u- and v-points)
434      DO jj = 2, jpjm1   ;   DO ji = 2, jpim1
435         zsshu_h(ji,jj) = 0.5_wp * ( ssh(ji,jj,Nnn) + ssh(ji+1,jj,Nnn) ) * ssumask(ji,jj)
436         zsshv_h(ji,jj) = 0.5_wp * ( ssh(ji,jj,Nnn) + ssh(ji,jj+1,Nnn) ) * ssvmask(ji,jj)
437      END DO             ;   END DO     
438      CALL lbc_lnk_multi( zsshu_h(:,:), 'U', 1._wp , zsshv_h(:,:), 'V', 1._wp )
439      !
440      !                                   !==  after depths and its inverse  ==!
441         hu_n(:,:) = hu_0(:,:) + zsshu_h(:,:)
442         hv_n(:,:) = hv_0(:,:) + zsshv_h(:,:)
443      r1_hu_n(:,:) = ssumask(:,:) / ( hu_n(:,:) + 1._wp - ssumask(:,:) )
444      r1_hv_n(:,:) = ssvmask(:,:) / ( hv_n(:,:) + 1._wp - ssvmask(:,:) )
445      !
446      !                                   !==  now scale factors  ==!  (e3t , e3u , e3v)
447      zssht_h(:,:) = ssh    (:,:,Nnn) * r1_ht_0(:,:)     ! t-point
448      zsshu_h(:,:) = zsshu_h(:,:)     * r1_hu_0(:,:)     ! u-point
449      zsshv_h(:,:) = zsshv_h(:,:)     * r1_hv_0(:,:)     ! v-point
450      DO jk = 1, jpkm1
451         e3t_n(:,:,jk) = e3t_0(:,:,jk) * ( 1._wp + zssht_h(:,:) * tmask(:,:,jk) )
452         e3u_n(:,:,jk) = e3u_0(:,:,jk) * ( 1._wp + zsshu_h(:,:) * umask(:,:,jk) )
453         e3v_n(:,:,jk) = e3v_0(:,:,jk) * ( 1._wp + zsshv_h(:,:) * vmask(:,:,jk) )
454         e3w_n(:,:,jk) = e3w_0(:,:,jk) * ( 1._wp + zssht_h(:,:) * MAX( tmask(:,:,jk) , tmask(:,:,jk+1) )
455      END DO
456      !
457      e3t_b(:,:,:)  = e3t_n(:,:,:)
458      e3u_b(:,:,:)  = e3u_n(:,:,:)
459      e3v_b(:,:,:)  = e3v_n(:,:,:)
460
461      ! t- and w- points depth
462      ! ----------------------
463      gdept_n(:,:,1) = 0.5_wp * e3w_n(:,:,1)
464      gdepw_n(:,:,1) = 0.0_wp
465      !
466      DO jk = 2, jpk
467         DO jj = 1,jpj
468            DO ji = 1,jpi
469               zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk))
470               gdepw_n(ji,jj,jk) = gdepw_n(ji,jj,jk-1) + e3t_n(ji,jj,jk-1)
471               gdept_n(ji,jj,jk) =      zcoef  * ( gdepw_n(ji,jj,jk  ) + 0.5 * e3w_n(ji,jj,jk))  &
472                  &                + (1-zcoef) * ( gdept_n(ji,jj,jk-1) + e3w_n(ji,jj,jk))
473            END DO
474         END DO
475      END DO
476      !
477      zssht_h(:,:) = 1._wp + zssht_h(:,:)               ! t-point
478      !
479      IF( ln_isfcav ) THEN    ! ISF cavities : ssh scaling not applied over the iceshelf thickness
480         DO jk = 1, jpkm1
481            gdept_n(:,:,jk) = ( gdept_0(:,:,jk) - risfdep(:,:) ) * zssht_h(:,:) + risfdep(:,:)
482            gdepw_n(:,:,jk) = ( gdepw_0(:,:,jk) - risfdep(:,:) ) * zssht_h(:,:) + risfdep(:,:)
483            gde3w_n(:,:,jk) =   gdept_n(:,:,jk) - ssh    (:,:,Nnn)
484         END DO
485      ELSE                    ! no ISF cavities
486         DO jk = 1, jpkm1
487            gdept_n(:,:,jk) = gdept_0(:,:,jk) * zssht_h(:,:)
488            gdepw_n(:,:,jk) = gdepw_0(:,:,jk) * zssht_h(:,:)
489            gde3w_n(:,:,jk) = gdept_n(:,:,jk) - ssh    (:,:,Nnn)
490         END DO
491      ENDIF
492      !
493      gdept_b(:,:,:) = gdept_n(:,:,:)
494      gdepw_b(:,:,:) = gdepw_n(:,:,:)
495      !
496   END SUBROUTINE dta_dyn_swp
497   
498
499   SUBROUTINE dta_dyn_ssh( kt, phdivtr, psshb,  pemp, pssha )
500      !!----------------------------------------------------------------------
501      !!                ***  ROUTINE dta_dyn_wzv  ***
502      !!                   
503      !! ** Purpose :   compute the after ssh (ssha) and the now vertical velocity
504      !!
505      !! ** Method  : Using the incompressibility hypothesis,
506      !!        - the ssh increment is computed by integrating the horizontal divergence
507      !!          and multiply by the time step.
508      !!
509      !!        - compute the after scale factor : repartition of ssh INCREMENT proportionnaly
510      !!                                           to the level thickness ( z-star case )
511      !!
512      !!        - the vertical velocity is computed by integrating the horizontal divergence 
513      !!          from the bottom to the surface minus the scale factor evolution.
514      !!          The boundary conditions are w=0 at the bottom (no flux)
515      !!
516      !! ** action  :   ssha / e3t_a / wn
517      !!
518      !! Reference  : Leclair, M., and G. Madec, 2009, Ocean Modelling.
519      !!----------------------------------------------------------------------
520      INTEGER,                                    INTENT(in   )    :: kt        !  time-step
521      REAL(wp), DIMENSION(jpi,jpj,jpk)          , INTENT(in   )   :: phdivtr   ! horizontal divergence transport
522      REAL(wp), DIMENSION(jpi,jpj)    , OPTIONAL, INTENT(in   )   :: psshb     ! now ssh
523      REAL(wp), DIMENSION(jpi,jpj)    , OPTIONAL, INTENT(in   )   :: pemp      ! evaporation minus precipitation
524      REAL(wp), DIMENSION(jpi,jpj)    , OPTIONAL, INTENT(inout) :: pssha     ! after ssh
525      !
526      INTEGER                       :: jk
527      REAL(wp), DIMENSION(jpi,jpj)  :: zhdiv 
528      !!----------------------------------------------------------------------
529      !
530      zhdiv(:,:) = 0._wp
531      DO jk = 1, jpkm1
532         zhdiv(:,:) = zhdiv(:,:) +  phdivtr(:,:,jk) * tmask(:,:,jk)
533      END DO
534      !                                                ! Sea surface  elevation time-stepping
535      pssha(:,:) = ( psshb(:,:) - rDt * ( r1_rho0 * pemp(:,:) + zhdiv(:,:) ) ) * ssmask(:,:)
536      !
537   END SUBROUTINE dta_dyn_ssh
538
539
540   SUBROUTINE dta_dyn_hrnf
541      !!----------------------------------------------------------------------
542      !!                  ***  ROUTINE sbc_rnf  ***
543      !!
544      !! ** Purpose :   update the horizontal divergence with the runoff inflow
545      !!
546      !! ** Method  :
547      !!                CAUTION : rnf is positive (inflow) decreasing the
548      !!                          divergence and expressed in m/s
549      !!
550      !! ** Action  :   phdivn   decreased by the runoff inflow
551      !!----------------------------------------------------------------------
552      !!
553      INTEGER  ::   ji, jj, jk   ! dummy loop indices
554      !!----------------------------------------------------------------------
555      !
556      DO jj = 1, jpj                   ! update the depth over which runoffs are distributed
557         DO ji = 1, jpi
558            h_rnf(ji,jj) = 0._wp
559            DO jk = 1, nk_rnf(ji,jj)                           ! recalculates h_rnf to be the depth in metres
560                h_rnf(ji,jj) = h_rnf(ji,jj) + e3t_n(ji,jj,jk)   ! to the bottom of the relevant grid box
561            END DO
562        END DO
563      END DO
564      !
565   END SUBROUTINE dta_dyn_hrnf
566
567
568
569   SUBROUTINE dta_dyn_slp( kt )
570      !!---------------------------------------------------------------------
571      !!                    ***  ROUTINE dta_dyn_slp  ***
572      !!
573      !! ** Purpose : Computation of slope
574      !!
575      !!---------------------------------------------------------------------
576      INTEGER,  INTENT(in) :: kt       ! time step
577      !
578      INTEGER  ::   ji, jj     ! dummy loop indices
579      REAL(wp) ::   ztinta     ! ratio applied to after  records when doing time interpolation
580      REAL(wp) ::   ztintb     ! ratio applied to before records when doing time interpolation
581      INTEGER  ::   iswap 
582      REAL(wp), DIMENSION(jpi,jpj,jpk)      ::   zuslp, zvslp, zwslpi, zwslpj
583      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts) ::   zts
584      !!---------------------------------------------------------------------
585      !
586      IF( sf_dyn(jf_tem)%ln_tint ) THEN    ! Computes slopes (here avt is used as workspace)                       
587         IF( kt == nit000 ) THEN
588            IF(lwp) WRITE(numout,*) ' Compute new slopes at kt = ', kt
589            zts(:,:,:,jp_tem) = sf_dyn(jf_tem)%fdta(:,:,:,1) * tmask(:,:,:)   ! temperature
590            zts(:,:,:,jp_sal) = sf_dyn(jf_sal)%fdta(:,:,:,1) * tmask(:,:,:)   ! salinity
591            avt(:,:,:)        = sf_dyn(jf_avt)%fdta(:,:,:,1) * tmask(:,:,:)   ! vertical diffusive coef.
592            CALL compute_slopes( kt, zts, zuslp, zvslp, zwslpi, zwslpj )
593            uslpdta (:,:,:,1) = zuslp (:,:,:) 
594            vslpdta (:,:,:,1) = zvslp (:,:,:) 
595            wslpidta(:,:,:,1) = zwslpi(:,:,:) 
596            wslpjdta(:,:,:,1) = zwslpj(:,:,:) 
597            !
598            zts(:,:,:,jp_tem) = sf_dyn(jf_tem)%fdta(:,:,:,2) * tmask(:,:,:)   ! temperature
599            zts(:,:,:,jp_sal) = sf_dyn(jf_sal)%fdta(:,:,:,2) * tmask(:,:,:)   ! salinity
600            avt(:,:,:)        = sf_dyn(jf_avt)%fdta(:,:,:,2) * tmask(:,:,:)   ! vertical diffusive coef.
601            CALL compute_slopes( kt, zts, zuslp, zvslp, zwslpi, zwslpj )
602            uslpdta (:,:,:,2) = zuslp (:,:,:) 
603            vslpdta (:,:,:,2) = zvslp (:,:,:) 
604            wslpidta(:,:,:,2) = zwslpi(:,:,:) 
605            wslpjdta(:,:,:,2) = zwslpj(:,:,:) 
606         ELSE
607           !
608           iswap = 0
609           IF( sf_dyn(jf_tem)%nrec_a(2) - nprevrec /= 0 )  iswap = 1
610           IF( nsecdyn > sf_dyn(jf_tem)%nrec_b(2) .AND. iswap == 1 )  THEN    ! read/update the after data
611              IF(lwp) WRITE(numout,*) ' Compute new slopes at kt = ', kt
612              uslpdta (:,:,:,1) =  uslpdta (:,:,:,2)         ! swap the data
613              vslpdta (:,:,:,1) =  vslpdta (:,:,:,2) 
614              wslpidta(:,:,:,1) =  wslpidta(:,:,:,2) 
615              wslpjdta(:,:,:,1) =  wslpjdta(:,:,:,2) 
616              !
617              zts(:,:,:,jp_tem) = sf_dyn(jf_tem)%fdta(:,:,:,2) * tmask(:,:,:)   ! temperature
618              zts(:,:,:,jp_sal) = sf_dyn(jf_sal)%fdta(:,:,:,2) * tmask(:,:,:)   ! salinity
619              avt(:,:,:)        = sf_dyn(jf_avt)%fdta(:,:,:,2) * tmask(:,:,:)   ! vertical diffusive coef.
620              CALL compute_slopes( kt, zts, zuslp, zvslp, zwslpi, zwslpj )
621              !
622              uslpdta (:,:,:,2) = zuslp (:,:,:) 
623              vslpdta (:,:,:,2) = zvslp (:,:,:) 
624              wslpidta(:,:,:,2) = zwslpi(:,:,:) 
625              wslpjdta(:,:,:,2) = zwslpj(:,:,:) 
626            ENDIF
627         ENDIF
628      ENDIF
629      !
630      IF( sf_dyn(jf_tem)%ln_tint )  THEN
631         ztinta =  REAL( nsecdyn - sf_dyn(jf_tem)%nrec_b(2), wp )  &
632            &    / REAL( sf_dyn(jf_tem)%nrec_a(2) - sf_dyn(jf_tem)%nrec_b(2), wp )
633         ztintb =  1. - ztinta
634         IF( l_ldfslp .AND. .NOT.lk_c1d ) THEN    ! Computes slopes (here avt is used as workspace)
635            uslp (:,:,:) = ztintb * uslpdta (:,:,:,1)  + ztinta * uslpdta (:,:,:,2) 
636            vslp (:,:,:) = ztintb * vslpdta (:,:,:,1)  + ztinta * vslpdta (:,:,:,2) 
637            wslpi(:,:,:) = ztintb * wslpidta(:,:,:,1)  + ztinta * wslpidta(:,:,:,2) 
638            wslpj(:,:,:) = ztintb * wslpjdta(:,:,:,1)  + ztinta * wslpjdta(:,:,:,2) 
639         ENDIF
640      ELSE
641         zts(:,:,:,jp_tem) = sf_dyn(jf_tem)%fnow(:,:,:) * tmask(:,:,:)   ! temperature
642         zts(:,:,:,jp_sal) = sf_dyn(jf_sal)%fnow(:,:,:) * tmask(:,:,:)   ! salinity
643         avt(:,:,:)        = sf_dyn(jf_avt)%fnow(:,:,:) * tmask(:,:,:)   ! vertical diffusive coef.
644         CALL compute_slopes( kt, zts, zuslp, zvslp, zwslpi, zwslpj )
645         !
646         IF( l_ldfslp .AND. .NOT.lk_c1d ) THEN    ! Computes slopes (here avt is used as workspace)
647            uslp (:,:,:) = zuslp (:,:,:)
648            vslp (:,:,:) = zvslp (:,:,:)
649            wslpi(:,:,:) = zwslpi(:,:,:)
650            wslpj(:,:,:) = zwslpj(:,:,:)
651         ENDIF
652      ENDIF
653      !
654   END SUBROUTINE dta_dyn_slp
655
656
657   SUBROUTINE compute_slopes( kt, pts, puslp, pvslp, pwslpi, pwslpj )
658      !!---------------------------------------------------------------------
659      !!                    ***  ROUTINE dta_dyn_slp  ***
660      !!
661      !! ** Purpose :   Computation of slope
662      !!---------------------------------------------------------------------
663      INTEGER ,                              INTENT(in ) :: kt       ! time step
664      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in ) :: pts      ! temperature/salinity
665      REAL(wp), DIMENSION(jpi,jpj,jpk)     , INTENT(out) :: puslp    ! zonal isopycnal slopes
666      REAL(wp), DIMENSION(jpi,jpj,jpk)     , INTENT(out) :: pvslp    ! meridional isopycnal slopes
667      REAL(wp), DIMENSION(jpi,jpj,jpk)     , INTENT(out) :: pwslpi   ! zonal diapycnal slopes
668      REAL(wp), DIMENSION(jpi,jpj,jpk)     , INTENT(out) :: pwslpj   ! meridional diapycnal slopes
669      !!---------------------------------------------------------------------
670      !
671      IF( l_ldfslp .AND. .NOT.lk_c1d ) THEN    ! Computes slopes (here avt is used as workspace)
672         CALL eos    ( pts, rhd, rhop, gdept_0(:,:,:) )
673         CALL eos_rab( pts, rab_n )       ! now local thermal/haline expension ratio at T-points
674         CALL bn2    ( pts, rab_n, rn2  ) ! now    Brunt-Vaisala
675
676      ! Partial steps: before Horizontal DErivative
677      IF( ln_zps  .AND. .NOT. ln_isfcav)                            &
678         &            CALL zps_hde    ( kt, jpts, pts, gtsu, gtsv,  &  ! Partial steps: before horizontal gradient
679         &                                        rhd, gru , grv    )  ! of t, s, rd at the last ocean level
680      IF( ln_zps .AND.        ln_isfcav)                            &
681         &            CALL zps_hde_isf( kt, jpts, pts, gtsu, gtsv, gtui, gtvi, &  ! Partial steps for top cell (ISF)
682         &                                        rhd, gru , grv , grui, grvi )  ! of t, s, rd at the first ocean level
683
684         rn2b(:,:,:) = rn2(:,:,:)         ! need for zdfmxl
685         CALL zdf_mxl( kt )            ! mixed layer depth
686         CALL ldf_slp( kt, rhd, rn2 )  ! slopes
687         puslp (:,:,:) = uslp (:,:,:)
688         pvslp (:,:,:) = vslp (:,:,:)
689         pwslpi(:,:,:) = wslpi(:,:,:)
690         pwslpj(:,:,:) = wslpj(:,:,:)
691     ELSE
692         puslp (:,:,:) = 0.            ! to avoid warning when compiling
693         pvslp (:,:,:) = 0.
694         pwslpi(:,:,:) = 0.
695         pwslpj(:,:,:) = 0.
696     ENDIF
697      !
698   END SUBROUTINE compute_slopes
699
700   !!======================================================================
701END MODULE dtadyn
Note: See TracBrowser for help on using the repository browser.