Context Navigation

dtadyn.F90 @ 12340

Last change on this file since 12340 was 12340, checked in by acc, 4 years ago

Branch 2019/dev_r11943_MERGE_2019. This commit introduces basic do loop macro
substitution to the 2019 option 1, merge branch. These changes have been SETTE
tested. The only addition is the do_loop_substitute.h90 file in the OCE directory but
the macros defined therein are used throughout the code to replace identifiable, 2D-
and 3D- nested loop opening and closing statements with single-line alternatives. Code
indents are also adjusted accordingly.

The following explanation is taken from comments in the new header file:

This header file contains preprocessor definitions and macros used in the do-loop
substitutions introduced between version 4.0 and 4.2. The primary aim of these macros
is to assist in future applications of tiling to improve performance. This is expected
to be achieved by alternative versions of these macros in selected locations. The
initial introduction of these macros simply replaces all identifiable nested 2D- and
3D-loops with single line statements (and adjusts indenting accordingly). Do loops
are identifiable if they comform to either:

DO jk = ....

DO jj = .... DO jj = ...

DO ji = .... DO ji = ...
. OR .
. .

END DO END DO

END DO END DO

END DO

and white-space variants thereof.

Additionally, only loops with recognised jj and ji loops limits are treated; these are:
Lower limits of 1, 2 or fs_2
Upper limits of jpi, jpim1 or fs_jpim1 (for ji) or jpj, jpjm1 or fs_jpjm1 (for jj)

The macro naming convention takes the form: DO_2D_BT_LR where:

B is the Bottom offset from the PE's inner domain;
T is the Top offset from the PE's inner domain;
L is the Left offset from the PE's inner domain;
R is the Right offset from the PE's inner domain

So, given an inner domain of 2,jpim1 and 2,jpjm1, a typical example would replace:

DO jj = 2, jpj

DO ji = 1, jpim1
.
.

END DO

END DO

with:

DO_2D_01_10
.
.
END_2D

similar conventions apply to the 3D loops macros. jk loop limits are retained
through macro arguments and are not restricted. This includes the possibility of
strides for which an extra set of DO_3DS macros are defined.

In the example definition below the inner PE domain is defined by start indices of
(kIs, kJs) and end indices of (kIe, KJe)

#define DO_2D_00_00 DO jj = kJs, kJe ; DO ji = kIs, kIe
#define END_2D END DO ; END DO

TO DO:

Only conventional nested loops have been identified and replaced by this step. There are constructs such as:

DO jk = 2, jpkm1

z2d(:,:) = z2d(:,:) + e3w(:,:,jk,Kmm) * z3d(:,:,jk) * wmask(:,:,jk)

END DO

which may need to be considered.

Property svn:keywords set to Id

File size: 41.5 KB

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

Note: See TracBrowser for help on using the repository browser.

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

Context Navigation

source: NEMO/branches/2019/dev_r11943_MERGE_2019/src/OFF/dtadyn.F90 @ 12340

Download in other formats: