Context Navigation

source: CONFIG/UNIFORM/v6/NEMO_v6.5/SOURCES/dtadyn.F90 @ 6230

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

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

Download in other formats: