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 branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/OFF_SRC – NEMO

Context Navigation

source: branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/OFF_SRC/dtadyn.F90 @ 7522

Last change on this file since 7522 was 7522, checked in by cetlod, 7 years ago
3.6 stable : update the offline routines to be able to run passive tracers offline with linear free surface, see ticket #1775
Property svn:keywords set to `Id`
File size: 35.6 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 ldfeiv ! eddy induced velocity coef.
33	USE ldftra_oce ! 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 wrk_nemo ! Memory allocation
44	USE timing ! Timing
45	USE trc, ONLY : ln_rsttr, numrtr, numrtw, lrst_trc
46
47	IMPLICIT NONE
48	PRIVATE
49
50	PUBLIC dta_dyn_init ! called by opa.F90
51	PUBLIC dta_dyn ! 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	!! * Substitutions
90	# include "domzgr_substitute.h90"
91	# include "vectopt_loop_substitute.h90"
92	!!----------------------------------------------------------------------
93	!! NEMO/OFF 3.3 , NEMO Consortium (2010)
94	!! $Id$
95	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
96	!!----------------------------------------------------------------------
97	CONTAINS
98
99	SUBROUTINE dta_dyn( kt )
100	!!----------------------------------------------------------------------
101	!! * ROUTINE dta_dyn *
102	!!
103	!! ** Purpose : Prepares dynamics and physics fields from a NEMO run
104	!! for an off-line simulation of passive tracers
105	!!
106	!! ** Method : calculates the position of data
107	!! - computes slopes if needed
108	!! - interpolates data if needed
109	!!----------------------------------------------------------------------
110	!
111	USE oce, ONLY: zhdivtr => ua
112	INTEGER, INTENT(in) :: kt ! ocean time-step index
113	INTEGER :: ji, jj, jk
114	REAL(wp), POINTER, DIMENSION(:,:) :: zemp
115	!
116	!!----------------------------------------------------------------------
117
118	!
119	IF( nn_timing == 1 ) CALL timing_start( 'dta_dyn')
120	!
121	!
122	nsecdyn = nsec_year + nsec1jan000 ! number of seconds between Jan. 1st 00h of nit000 year and the middle of time step
123	!
124	IF( kt == nit000 ) THEN ; nprevrec = 0
125	ELSE ; nprevrec = sf_dyn(jf_tem)%nrec_a(2)
126	ENDIF
127	!
128	CALL fld_read( kt, 1, sf_dyn ) != read data at kt time step ==!
129	!
130	IF( lk_ldfslp .AND. .NOT.lk_c1d ) CALL dta_dyn_slp( kt ) ! Computation of slopes
131	!
132	tsn(:,:,:,jp_tem) = sf_dyn(jf_tem)%fnow(:,:,:) * tmask(:,:,:) ! temperature
133	tsn(:,:,:,jp_sal) = sf_dyn(jf_sal)%fnow(:,:,:) * tmask(:,:,:) ! salinity
134	wndm(:,:) = sf_dyn(jf_wnd)%fnow(:,:,1) * tmask(:,:,1) ! wind speed - needed for gas exchange
135	fmmflx(:,:) = sf_dyn(jf_fmf)%fnow(:,:,1) * tmask(:,:,1) ! downward salt flux (v3.5+)
136	fr_i(:,:) = sf_dyn(jf_ice)%fnow(:,:,1) * tmask(:,:,1) ! Sea-ice fraction
137	qsr (:,:) = sf_dyn(jf_qsr)%fnow(:,:,1) * tmask(:,:,1) ! solar radiation
138	emp (:,:) = sf_dyn(jf_emp)%fnow(:,:,1) * tmask(:,:,1) ! E-P
139	IF( ln_dynrnf ) THEN
140	rnf (:,:) = sf_dyn(jf_rnf)%fnow(:,:,1) * tmask(:,:,1) ! E-P
141	IF( ln_dynrnf_depth .AND. lk_vvl ) CALL dta_dyn_hrnf
142	ENDIF
143	!
144	un(:,:,:) = sf_dyn(jf_uwd)%fnow(:,:,:) * umask(:,:,:) ! effective u-transport
145	vn(:,:,:) = sf_dyn(jf_vwd)%fnow(:,:,:) * vmask(:,:,:) ! effective v-transport
146	wn(:,:,:) = sf_dyn(jf_wwd)%fnow(:,:,:) * tmask(:,:,:) ! effective v-transport
147	!
148	IF( lk_vvl ) THEN
149	CALL wrk_alloc(jpi, jpj, zemp )
150	zhdivtr(:,:,:) = sf_dyn(jf_div)%fnow(:,:,:) * tmask(:,:,:) ! effective u-transport
151	emp_b (:,:) = sf_dyn(jf_empb)%fnow(:,:,1) * tmask(:,:,1) ! E-P
152	zemp(:,:) = 0.5_wp * ( emp(:,:) + emp_b(:,:) ) + rnf(:,:) + fwbcorr * tmask(:,:,1)
153	CALL dta_dyn_ssh( kt, zhdivtr, sshb, zemp, ssha, fse3t_a(:,:,:) ) != ssh, vertical scale factor & vertical transport
154	CALL wrk_dealloc(jpi, jpj, zemp )
155	! Write in the tracer restart file
156	! *******************************
157	IF( lrst_trc ) THEN
158	IF(lwp) WRITE(numout,*)
159	IF(lwp) WRITE(numout,*) 'dta_dyn_ssh : ssh field written in tracer restart file ', &
160	& 'at it= ', kt,' date= ', ndastp
161	IF(lwp) WRITE(numout,*) '~~~~'
162	CALL iom_rstput( kt, nitrst, numrtw, 'sshn', ssha )
163	CALL iom_rstput( kt, nitrst, numrtw, 'sshb', sshn )
164	ENDIF
165	ENDIF
166	!
167	CALL eos ( tsn, rhd, rhop, gdept_0(:,:,:) ) ! In any case, we need rhop
168	CALL eos_rab( tsn, rab_n ) ! now local thermal/haline expension ratio at T-points
169	CALL bn2 ( tsn, rab_n, rn2 ) ! before Brunt-Vaisala frequency need for zdfmxl
170
171	rn2b(:,:,:) = rn2(:,:,:) ! need for zdfmxl
172	CALL zdf_mxl( kt ) ! In any case, we need mxl
173	!
174	hmld(:,:) = sf_dyn(jf_mld)%fnow(:,:,1) * tmask(:,:,1) ! mixed layer depht
175	avt(:,:,:) = sf_dyn(jf_avt)%fnow(:,:,:) * tmask(:,:,:) ! vertical diffusive coefficient
176	!
177	#if defined key_trabbl && ! defined key_c1d
178	ahu_bbl(:,:) = sf_dyn(jf_ubl)%fnow(:,:,1) * umask(:,:,1) ! bbl diffusive coef
179	ahv_bbl(:,:) = sf_dyn(jf_vbl)%fnow(:,:,1) * vmask(:,:,1)
180	#endif
181	!
182	!
183	CALL eos( tsn, rhd, rhop, gdept_0(:,:,:) ) ! In any case, we need rhop
184	!
185	IF(ln_ctl) THEN ! print control
186	CALL prt_ctl(tab3d_1=tsn(:,:,:,jp_tem), clinfo1=' tn - : ', mask1=tmask, ovlap=1, kdim=jpk )
187	CALL prt_ctl(tab3d_1=tsn(:,:,:,jp_sal), clinfo1=' sn - : ', mask1=tmask, ovlap=1, kdim=jpk )
188	CALL prt_ctl(tab3d_1=un , clinfo1=' un - : ', mask1=umask, ovlap=1, kdim=jpk )
189	CALL prt_ctl(tab3d_1=vn , clinfo1=' vn - : ', mask1=vmask, ovlap=1, kdim=jpk )
190	CALL prt_ctl(tab3d_1=wn , clinfo1=' wn - : ', mask1=tmask, ovlap=1, kdim=jpk )
191	CALL prt_ctl(tab3d_1=avt , clinfo1=' kz - : ', mask1=tmask, ovlap=1, kdim=jpk )
192	! CALL prt_ctl(tab2d_1=fr_i , clinfo1=' fr_i - : ', mask1=tmask, ovlap=1 )
193	! CALL prt_ctl(tab2d_1=hmld , clinfo1=' hmld - : ', mask1=tmask, ovlap=1 )
194	! CALL prt_ctl(tab2d_1=fmmflx , clinfo1=' fmmflx - : ', mask1=tmask, ovlap=1 )
195	! CALL prt_ctl(tab2d_1=emp , clinfo1=' emp - : ', mask1=tmask, ovlap=1 )
196	! CALL prt_ctl(tab2d_1=wndm , clinfo1=' wspd - : ', mask1=tmask, ovlap=1 )
197	! CALL prt_ctl(tab2d_1=qsr , clinfo1=' qsr - : ', mask1=tmask, ovlap=1 )
198	ENDIF
199	!
200	IF( nn_timing == 1 ) CALL timing_stop( 'dta_dyn')
201	!
202	END SUBROUTINE dta_dyn
203
204
205	SUBROUTINE dta_dyn_init
206	!!----------------------------------------------------------------------
207	!! * ROUTINE dta_dyn_init *
208	!!
209	!! ** Purpose : Initialisation of the dynamical data
210	!! ** Method : - read the data namdta_dyn namelist
211	!!
212	!! ** Action : - read parameters
213	!!----------------------------------------------------------------------
214	INTEGER :: ierr, ierr0, ierr1, ierr2, ierr3 ! return error code
215	INTEGER :: ifpr ! dummy loop indice
216	INTEGER :: jfld ! dummy loop arguments
217	INTEGER :: inum, idv, idimv ! local integer
218	INTEGER :: ios ! Local integer output status for namelist read
219	INTEGER :: ji, jj, jk
220	REAL(wp) :: zcoef
221	INTEGER :: nkrnf_max
222	REAL(wp) :: hrnf_max
223	!!
224	CHARACTER(len=100) :: cn_dir ! Root directory for location of core files
225	TYPE(FLD_N), DIMENSION(jpfld) :: slf_d ! array of namelist informations on the fields to read
226	TYPE(FLD_N) :: sn_uwd, sn_vwd, sn_wwd, sn_empb, sn_emp ! informations about the fields to be read
227	TYPE(FLD_N) :: sn_tem , sn_sal , sn_avt ! " "
228	TYPE(FLD_N) :: sn_mld, sn_qsr, sn_wnd , sn_ice , sn_fmf ! " "
229	TYPE(FLD_N) :: sn_ubl, sn_vbl, sn_rnf ! " "
230	TYPE(FLD_N) :: sn_div ! informations about the fields to be read
231	!!----------------------------------------------------------------------
232
233	NAMELIST/namdta_dyn/cn_dir, ln_dynrnf, ln_dynrnf_depth, fwbcorr, &
234	& sn_uwd, sn_vwd, sn_wwd, sn_emp, &
235	& sn_avt, sn_tem, sn_sal, sn_mld , sn_qsr , &
236	& sn_wnd, sn_ice, sn_fmf, &
237	& sn_ubl, sn_vbl, sn_rnf, &
238	& sn_empb, sn_div
239	!
240	REWIND( numnam_ref ) ! Namelist namdta_dyn in reference namelist : Offline: init. of dynamical data
241	READ ( numnam_ref, namdta_dyn, IOSTAT = ios, ERR = 901)
242	901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdta_dyn in reference namelist', lwp )
243
244	REWIND( numnam_cfg ) ! Namelist namdta_dyn in configuration namelist : Offline: init. of dynamical data
245	READ ( numnam_cfg, namdta_dyn, IOSTAT = ios, ERR = 902 )
246	902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdta_dyn in configuration namelist', lwp )
247	IF(lwm) WRITE ( numond, namdta_dyn )
248	! ! store namelist information in an array
249	! ! Control print
250	IF(lwp) THEN
251	WRITE(numout,*)
252	WRITE(numout,*) 'dta_dyn : offline dynamics '
253	WRITE(numout,*) '~~~~~~~ '
254	WRITE(numout,*) ' Namelist namdta_dyn'
255	WRITE(numout,*) ' runoffs option enabled (T) or not (F) ln_dynrnf = ', ln_dynrnf
256	WRITE(numout,*) ' runoffs is spread in vertical ln_dynrnf_depth = ', ln_dynrnf_depth
257	WRITE(numout,*) ' annual global mean of empmr for ssh correction fwbcorr = ', fwbcorr
258	WRITE(numout,*)
259	ENDIF
260	!
261
262	jf_uwd = 1 ; jf_vwd = 2 ; jf_wwd = 3 ; jf_emp = 4 ; jf_avt = 5
263	jf_tem = 6 ; jf_sal = 7 ; jf_mld = 8 ; jf_qsr = 9
264	jf_wnd = 10 ; jf_ice = 11 ; jf_fmf = 12 ; jfld = jf_fmf
265
266	!
267	slf_d(jf_uwd) = sn_uwd ; slf_d(jf_vwd) = sn_vwd ; slf_d(jf_wwd) = sn_wwd
268	slf_d(jf_emp) = sn_emp ; slf_d(jf_avt) = sn_avt
269	slf_d(jf_tem) = sn_tem ; slf_d(jf_sal) = sn_sal ; slf_d(jf_mld) = sn_mld
270	slf_d(jf_qsr) = sn_qsr ; slf_d(jf_wnd) = sn_wnd ; slf_d(jf_ice) = sn_ice
271	slf_d(jf_fmf) = sn_fmf
272
273
274	!
275	IF( lk_vvl ) THEN
276	jf_div = jfld + 1 ; jf_empb = jfld + 2 ; jfld = jf_empb
277	slf_d(jf_div) = sn_div ; slf_d(jf_empb) = sn_empb
278	ENDIF
279	!
280	IF( lk_trabbl ) THEN
281	jf_ubl = jfld + 1 ; jf_vbl = jfld + 2 ; jfld = jf_vbl
282	slf_d(jf_ubl) = sn_ubl ; slf_d(jf_vbl) = sn_vbl
283	ENDIF
284	!
285	IF( ln_dynrnf ) THEN
286	jf_rnf = jfld + 1 ; jfld = jf_rnf
287	slf_d(jf_rnf) = sn_rnf
288	ELSE
289	rnf(:,:) = 0._wp
290	ENDIF
291
292
293	ALLOCATE( sf_dyn(jfld), STAT=ierr ) ! set sf structure
294	IF( ierr > 0 ) THEN
295	CALL ctl_stop( 'dta_dyn: unable to allocate sf structure' ) ; RETURN
296	ENDIF
297	! ! fill sf with slf_i and control print
298	CALL fld_fill( sf_dyn, slf_d, cn_dir, 'dta_dyn_init', 'Data in file', 'namdta_dyn' )
299	!
300	! Open file for each variable to get his number of dimension
301	DO ifpr = 1, jfld
302	CALL fld_clopn( sf_dyn(ifpr), nyear, nmonth, nday )
303	idv = iom_varid( sf_dyn(ifpr)%num , slf_d(ifpr)%clvar ) ! id of the variable sdjf%clvar
304	idimv = iom_file ( sf_dyn(ifpr)%num )%ndims(idv) ! number of dimension for variable sdjf%clvar
305	IF( sf_dyn(ifpr)%num /= 0 ) CALL iom_close( sf_dyn(ifpr)%num ) ! close file if already open
306	ierr1=0
307	IF( idimv == 3 ) THEN ! 2D variable
308	ALLOCATE( sf_dyn(ifpr)%fnow(jpi,jpj,1) , STAT=ierr0 )
309	IF( slf_d(ifpr)%ln_tint ) ALLOCATE( sf_dyn(ifpr)%fdta(jpi,jpj,1,2) , STAT=ierr1 )
310	ELSE ! 3D variable
311	ALLOCATE( sf_dyn(ifpr)%fnow(jpi,jpj,jpk) , STAT=ierr0 )
312	IF( slf_d(ifpr)%ln_tint ) ALLOCATE( sf_dyn(ifpr)%fdta(jpi,jpj,jpk,2), STAT=ierr1 )
313	ENDIF
314	IF( ierr0 + ierr1 > 0 ) THEN
315	CALL ctl_stop( 'dta_dyn_init : unable to allocate sf_dyn array structure' ) ; RETURN
316	ENDIF
317	END DO
318	!
319	IF( lk_ldfslp .AND. .NOT.lk_c1d ) THEN ! slopes
320	IF( sf_dyn(jf_tem)%ln_tint ) THEN ! time interpolation
321	ALLOCATE( uslpdta (jpi,jpj,jpk,2), vslpdta (jpi,jpj,jpk,2), &
322	& wslpidta(jpi,jpj,jpk,2), wslpjdta(jpi,jpj,jpk,2), STAT=ierr2 )
323	!
324	IF( ierr2 > 0 ) THEN
325	CALL ctl_stop( 'dta_dyn_init : unable to allocate slope arrays' ) ; RETURN
326	ENDIF
327	ENDIF
328	ENDIF
329	!
330	IF( lk_vvl ) THEN
331	IF( .NOT. sf_dyn(jf_uwd)%ln_clim .AND. ln_rsttr .AND. & ! Restart: read in restart file
332	iom_varid( numrtr, 'sshn', ldstop = .FALSE. ) > 0 ) THEN
333	IF(lwp) WRITE(numout,*) ' sshn forcing fields read in the restart file for initialisation'
334	CALL iom_get( numrtr, jpdom_autoglo, 'sshn', sshn(:,:) )
335	CALL iom_get( numrtr, jpdom_autoglo, 'sshb', sshb(:,:) )
336	ELSE
337	IF(lwp) WRITE(numout,*) ' sshn forcing fields read in the restart file for initialisation'
338	CALL iom_open( 'restart', inum )
339	CALL iom_get( inum, jpdom_autoglo, 'sshn', sshn(:,:) )
340	CALL iom_get( inum, jpdom_autoglo, 'sshb', sshb(:,:) )
341	CALL iom_close( inum ) ! close file
342	ENDIF
343	!
344	DO jk = 1, jpkm1
345	fse3t_n(:,:,jk) = e3t_0(:,:,jk) * ( 1._wp + sshn(:,:) * tmask(:,:,1) / ( ht_0(:,:) + 1.0 - tmask(:,:,1) ) )
346	ENDDO
347	fse3t_a(:,:,jpk) = e3t_0(:,:,jpk)
348
349	! Horizontal scale factor interpolations
350	! --------------------------------------
351	CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3u_n(:,:,:), 'U' )
352	CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3v_n(:,:,:), 'V' )
353
354	! Vertical scale factor interpolations
355	! ------------------------------------
356	CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3w_n(:,:,:), 'W' )
357
358	fse3t_b(:,:,:) = fse3t_n(:,:,:)
359	fse3u_b(:,:,:) = fse3u_n(:,:,:)
360	fse3v_b(:,:,:) = fse3v_n(:,:,:)
361
362	! t- and w- points depth
363	! ----------------------
364	fsdept_n(:,:,1) = 0.5_wp * fse3w_n(:,:,1)
365	fsdepw_n(:,:,1) = 0.0_wp
366
367	DO jk = 2, jpk
368	DO jj = 1,jpj
369	DO ji = 1,jpi
370	! zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk)) ! 0 everywhere
371	! tmask = wmask, ie everywhere expect at jk = mikt
372	! 1 for jk =
373	! mikt
374	zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk))
375	fsdepw_n(ji,jj,jk) = fsdepw_n(ji,jj,jk-1) + fse3t_n(ji,jj,jk-1)
376	fsdept_n(ji,jj,jk) = zcoef * ( fsdepw_n(ji,jj,jk ) + 0.5 * fse3w_n(ji,jj,jk)) &
377	& + (1-zcoef) * ( fsdept_n(ji,jj,jk-1) + fse3w_n(ji,jj,jk))
378	END DO
379	END DO
380	END DO
381
382	fsdept_b(:,:,:) = fsdept_n(:,:,:)
383	fsdepw_b(:,:,:) = fsdepw_n(:,:,:)
384	!
385	ENDIF
386	!
387	IF( ln_dynrnf .AND. ln_dynrnf_depth ) THEN ! read depht over which runoffs are distributed
388	IF(lwp) WRITE(numout,*)
389	IF(lwp) WRITE(numout,*) ' read in the file depht over which runoffs are distributed'
390	CALL iom_open ( "runoffs", inum ) ! open file
391	CALL iom_get ( inum, jpdom_data, 'rodepth', h_rnf ) ! read the river mouth array
392	CALL iom_close( inum ) ! close file
393	!
394	nk_rnf(:,:) = 0 ! set the number of level over which river runoffs are applied
395	DO jj = 1, jpj
396	DO ji = 1, jpi
397	IF( h_rnf(ji,jj) > 0._wp ) THEN
398	jk = 2
399	DO WHILE ( jk /= mbkt(ji,jj) .AND. gdept_0(ji,jj,jk) < h_rnf(ji,jj) ) ; jk = jk + 1
400	END DO
401	nk_rnf(ji,jj) = jk
402	ELSEIF( h_rnf(ji,jj) == -1._wp ) THEN ; nk_rnf(ji,jj) = 1
403	ELSEIF( h_rnf(ji,jj) == -999._wp ) THEN ; nk_rnf(ji,jj) = mbkt(ji,jj)
404	ELSE
405	CALL ctl_stop( 'sbc_rnf_init: runoff depth not positive, and not -999 or -1, rnf value in file fort.999' )
406	WRITE(999,*) 'ji, jj, h_rnf(ji,jj) :', ji, jj, h_rnf(ji,jj)
407	ENDIF
408	END DO
409	END DO
410	DO jj = 1, jpj ! set the associated depth
411	DO ji = 1, jpi
412	h_rnf(ji,jj) = 0._wp
413	DO jk = 1, nk_rnf(ji,jj)
414	h_rnf(ji,jj) = h_rnf(ji,jj) + fse3t(ji,jj,jk)
415	END DO
416	END DO
417	END DO
418	ELSE ! runoffs applied at the surface
419	nk_rnf(:,:) = 1
420	h_rnf (:,:) = fse3t(:,:,1)
421	ENDIF
422	nkrnf_max = MAXVAL( nk_rnf(:,:) )
423	hrnf_max = MAXVAL( h_rnf(:,:) )
424	IF( lk_mpp ) THEN
425	CALL mpp_max( nkrnf_max ) ! max over the global domain
426	CALL mpp_max( hrnf_max ) ! max over the global domain
427	ENDIF
428	IF(lwp) WRITE(numout,*) ' '
429	IF(lwp) WRITE(numout,*) ' max depht of runoff : ', hrnf_max,' max level : ', nkrnf_max
430	IF(lwp) WRITE(numout,*) ' '
431	!
432	CALL dta_dyn( nit000 )
433	!
434	END SUBROUTINE dta_dyn_init
435
436	SUBROUTINE dta_dyn_swp( kt )
437	!!---------------------------------------------------------------------
438	!! * ROUTINE dta_dyn_swp *
439	!!
440	!! ** Purpose : Swap and the data and compute the vertical scale factor at U/V/W point
441	!! and the depht
442	!!
443	!!---------------------------------------------------------------------
444	INTEGER, INTENT(in) :: kt ! time step
445	INTEGER :: ji, jj, jk
446	REAL(wp) :: zcoef
447	!
448	!!---------------------------------------------------------------------
449
450	IF( kt == nit000 ) THEN
451	IF(lwp) WRITE(numout,*)
452	IF(lwp) WRITE(numout,*) 'ssh_swp : Asselin time filter and swap of sea surface height'
453	IF(lwp) WRITE(numout,*) '~~~~~~~ '
454	ENDIF
455
456	sshb(:,:) = sshn(:,:) + atfp * ( sshb(:,:) - 2 * sshn(:,:) + ssha(:,:)) ! before <-- now filtered
457	sshn(:,:) = ssha(:,:)
458
459	fse3t_n(:,:,:) = fse3t_a(:,:,:)
460
461	! Reconstruction of all vertical scale factors at now and before time steps
462	! =============================================================================
463
464	! Horizontal scale factor interpolations
465	! --------------------------------------
466	CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3u_n(:,:,:), 'U' )
467	CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3v_n(:,:,:), 'V' )
468
469	! Vertical scale factor interpolations
470	! ------------------------------------
471	CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3w_n (:,:,:), 'W' )
472
473	fse3t_b(:,:,:) = fse3t_n(:,:,:)
474	fse3u_b(:,:,:) = fse3u_n(:,:,:)
475	fse3v_b(:,:,:) = fse3v_n(:,:,:)
476
477	! t- and w- points depth
478	! ----------------------
479	fsdept_n(:,:,1) = 0.5_wp * fse3w_n(:,:,1)
480	fsdepw_n(:,:,1) = 0.0_wp
481
482	DO jk = 2, jpk
483	DO jj = 1,jpj
484	DO ji = 1,jpi
485	zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk))
486	fsdepw_n(ji,jj,jk) = fsdepw_n(ji,jj,jk-1) + fse3t_n(ji,jj,jk-1)
487	fsdept_n(ji,jj,jk) = zcoef * ( fsdepw_n(ji,jj,jk ) + 0.5 * fse3w_n(ji,jj,jk)) &
488	& + (1-zcoef) * ( fsdept_n(ji,jj,jk-1) + fse3w_n(ji,jj,jk))
489	END DO
490	END DO
491	END DO
492
493	fsdept_b(:,:,:) = fsdept_n(:,:,:)
494	fsdepw_b(:,:,:) = fsdepw_n(:,:,:)
495
496	!
497	END SUBROUTINE dta_dyn_swp
498
499	SUBROUTINE dta_dyn_ssh( kt, phdivtr, psshb, pemp, pssha, pe3ta )
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	!! * Arguments
521	INTEGER, INTENT(in ) :: kt ! time-step
522	REAL(wp), DIMENSION(jpi,jpj,jpk) , INTENT(in ) :: phdivtr ! horizontal divergence transport
523	REAL(wp), DIMENSION(jpi,jpj) , OPTIONAL, INTENT(in ) :: psshb ! now ssh
524	REAL(wp), DIMENSION(jpi,jpj) , OPTIONAL, INTENT(in ) :: pemp ! evaporation minus precipitation
525	REAL(wp), DIMENSION(jpi,jpj) , OPTIONAL, INTENT(inout) :: pssha ! after ssh
526	REAL(wp), DIMENSION(jpi,jpj,jpk), OPTIONAL, INTENT(out) :: pe3ta ! after vertical scale factor
527	!! * Local declarations
528	INTEGER :: jk
529	REAL(wp), DIMENSION(jpi,jpj) :: zhdiv
530	REAL(wp) :: z2dt
531	!!----------------------------------------------------------------------
532
533	!
534	z2dt = 2._wp * rdt
535	!
536	zhdiv(:,:) = 0._wp
537	DO jk = 1, jpkm1
538	zhdiv(:,:) = zhdiv(:,:) + phdivtr(:,:,jk) * tmask(:,:,jk)
539	END DO
540	! ! Sea surface elevation time-stepping
541	pssha(:,:) = ( psshb(:,:) - z2dt * ( r1_rau0 * pemp(:,:) + zhdiv(:,:) ) ) * ssmask(:,:)
542	! !
543	! ! After acale factors at t-points ( z_star coordinate )
544	DO jk = 1, jpkm1
545	pe3ta(:,:,jk) = e3t_0(:,:,jk) * ( 1._wp + pssha(:,:) * tmask(:,:,1) / ( ht_0(:,:) + 1.0 - tmask(:,:,1) ) )
546	END DO
547	!
548	END SUBROUTINE dta_dyn_ssh
549
550
551	SUBROUTINE dta_dyn_hrnf
552	!!----------------------------------------------------------------------
553	!! * ROUTINE sbc_rnf *
554	!!
555	!! ** Purpose : update the horizontal divergence with the runoff inflow
556	!!
557	!! ** Method :
558	!! CAUTION : rnf is positive (inflow) decreasing the
559	!! divergence and expressed in m/s
560	!!
561	!! ** Action : phdivn decreased by the runoff inflow
562	!!----------------------------------------------------------------------
563	!!
564	INTEGER :: ji, jj, jk ! dummy loop indices
565	!!----------------------------------------------------------------------
566	!
567	DO jj = 1, jpj ! update the depth over which runoffs are distributed
568	DO ji = 1, jpi
569	h_rnf(ji,jj) = 0._wp
570	DO jk = 1, nk_rnf(ji,jj) ! recalculates h_rnf to be the depth in metres
571	h_rnf(ji,jj) = h_rnf(ji,jj) + fse3t(ji,jj,jk) ! to the bottom of the relevant grid box
572	END DO
573	END DO
574	END DO
575	!
576	END SUBROUTINE dta_dyn_hrnf
577
578
579
580	SUBROUTINE dta_dyn_slp( kt )
581	!!---------------------------------------------------------------------
582	!! * ROUTINE dta_dyn_slp *
583	!!
584	!! ** Purpose : Computation of slope
585	!!
586	!!---------------------------------------------------------------------
587	USE oce, ONLY: zts => tsa
588	!
589	INTEGER, INTENT(in) :: kt ! time step
590	!
591	INTEGER :: ji, jj ! dummy loop indices
592	REAL(wp) :: ztinta ! ratio applied to after records when doing time interpolation
593	REAL(wp) :: ztintb ! ratio applied to before records when doing time interpolation
594	INTEGER :: iswap
595	REAL(wp), POINTER, DIMENSION(:,:,:) :: zuslp, zvslp, zwslpi, zwslpj
596	!!---------------------------------------------------------------------
597	!
598	CALL wrk_alloc(jpi, jpj, jpk, zuslp, zvslp, zwslpi, zwslpj )
599	!
600	IF( sf_dyn(jf_tem)%ln_tint ) THEN ! Computes slopes (here avt is used as workspace)
601	IF( kt == nit000 ) THEN
602	zts(:,:,:,jp_tem) = sf_dyn(jf_tem)%fdta(:,:,:,1) * tmask(:,:,:) ! temperature
603	zts(:,:,:,jp_sal) = sf_dyn(jf_sal)%fdta(:,:,:,1) * tmask(:,:,:) ! salinity
604	avt(:,:,:) = sf_dyn(jf_avt)%fdta(:,:,:,1) * tmask(:,:,:) ! vertical diffusive coef.
605	CALL compute_slopes( kt, zts, zuslp, zvslp, zwslpi, zwslpj )
606	uslpdta (:,:,:,1) = zuslp (:,:,:)
607	vslpdta (:,:,:,1) = zvslp (:,:,:)
608	wslpidta(:,:,:,1) = zwslpi(:,:,:)
609	wslpjdta(:,:,:,1) = zwslpj(:,:,:)
610	!
611	zts(:,:,:,jp_tem) = sf_dyn(jf_tem)%fdta(:,:,:,2) * tmask(:,:,:) ! temperature
612	zts(:,:,:,jp_sal) = sf_dyn(jf_sal)%fdta(:,:,:,2) * tmask(:,:,:) ! salinity
613	avt(:,:,:) = sf_dyn(jf_avt)%fdta(:,:,:,2) * tmask(:,:,:) ! vertical diffusive coef.
614	CALL compute_slopes( kt, zts, zuslp, zvslp, zwslpi, zwslpj )
615	uslpdta (:,:,:,2) = zuslp (:,:,:)
616	vslpdta (:,:,:,2) = zvslp (:,:,:)
617	wslpidta(:,:,:,2) = zwslpi(:,:,:)
618	wslpjdta(:,:,:,2) = zwslpj(:,:,:)
619	ELSE
620	!
621	iswap = 0
622	IF( sf_dyn(jf_tem)%nrec_a(2) - nprevrec /= 0 ) iswap = 1
623	IF( nsecdyn > sf_dyn(jf_tem)%nrec_b(2) .AND. iswap == 1 ) THEN ! read/update the after data
624	IF(lwp) WRITE(numout,*) ' Compute new slopes at kt = ', kt
625	uslpdta (:,:,:,1) = uslpdta (:,:,:,2) ! swap the data
626	vslpdta (:,:,:,1) = vslpdta (:,:,:,2)
627	wslpidta(:,:,:,1) = wslpidta(:,:,:,2)
628	wslpjdta(:,:,:,1) = wslpjdta(:,:,:,2)
629	!
630	zts(:,:,:,jp_tem) = sf_dyn(jf_tem)%fdta(:,:,:,2) * tmask(:,:,:) ! temperature
631	zts(:,:,:,jp_sal) = sf_dyn(jf_sal)%fdta(:,:,:,2) * tmask(:,:,:) ! salinity
632	avt(:,:,:) = sf_dyn(jf_avt)%fdta(:,:,:,2) * tmask(:,:,:) ! vertical diffusive coef.
633	CALL compute_slopes( kt, zts, zuslp, zvslp, zwslpi, zwslpj )
634	!
635	uslpdta (:,:,:,2) = zuslp (:,:,:)
636	vslpdta (:,:,:,2) = zvslp (:,:,:)
637	wslpidta(:,:,:,2) = zwslpi(:,:,:)
638	wslpjdta(:,:,:,2) = zwslpj(:,:,:)
639	ENDIF
640	ENDIF
641	ENDIF
642	!
643	IF( sf_dyn(jf_tem)%ln_tint ) THEN
644	ztinta = REAL( nsecdyn - sf_dyn(jf_tem)%nrec_b(2), wp ) &
645	& / REAL( sf_dyn(jf_tem)%nrec_a(2) - sf_dyn(jf_tem)%nrec_b(2), wp )
646	ztintb = 1. - ztinta
647	#if defined key_ldfslp && ! defined key_c1d
648	uslp (:,:,:) = ztintb * uslpdta (:,:,:,1) + ztinta * uslpdta (:,:,:,2)
649	vslp (:,:,:) = ztintb * vslpdta (:,:,:,1) + ztinta * vslpdta (:,:,:,2)
650	wslpi(:,:,:) = ztintb * wslpidta(:,:,:,1) + ztinta * wslpidta(:,:,:,2)
651	wslpj(:,:,:) = ztintb * wslpjdta(:,:,:,1) + ztinta * wslpjdta(:,:,:,2)
652	#endif
653	ELSE
654	zts(:,:,:,jp_tem) = sf_dyn(jf_tem)%fnow(:,:,:) * tmask(:,:,:) ! temperature
655	zts(:,:,:,jp_sal) = sf_dyn(jf_sal)%fnow(:,:,:) * tmask(:,:,:) ! salinity
656	avt(:,:,:) = sf_dyn(jf_avt)%fnow(:,:,:) * tmask(:,:,:) ! vertical diffusive coef.
657	CALL compute_slopes( kt, zts, zuslp, zvslp, zwslpi, zwslpj )
658	!
659	#if defined key_ldfslp && ! defined key_c1d
660	uslp (:,:,:) = zuslp (:,:,:)
661	vslp (:,:,:) = zvslp (:,:,:)
662	wslpi(:,:,:) = zwslpi(:,:,:)
663	wslpj(:,:,:) = zwslpj(:,:,:)
664	#endif
665	ENDIF
666	!
667	CALL wrk_dealloc(jpi, jpj, jpk, zuslp, zvslp, zwslpi, zwslpj )
668	!
669	END SUBROUTINE dta_dyn_slp
670
671	SUBROUTINE compute_slopes( kt, pts, puslp, pvslp, pwslpi, pwslpj )
672	!!---------------------------------------------------------------------
673	!! * ROUTINE dta_dyn_slp *
674	!!
675	!! ** Purpose : Computation of slope
676	!!
677	!!---------------------------------------------------------------------
678	INTEGER , INTENT(in ) :: kt ! time step
679	REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in ) :: pts ! temperature/salinity
680	REAL(wp), DIMENSION(jpi,jpj,jpk) , INTENT(out) :: puslp ! zonal isopycnal slopes
681	REAL(wp), DIMENSION(jpi,jpj,jpk) , INTENT(out) :: pvslp ! meridional isopycnal slopes
682	REAL(wp), DIMENSION(jpi,jpj,jpk) , INTENT(out) :: pwslpi ! zonal diapycnal slopes
683	REAL(wp), DIMENSION(jpi,jpj,jpk) , INTENT(out) :: pwslpj ! meridional diapycnal slopes
684	!!---------------------------------------------------------------------
685	#if defined key_ldfslp && ! defined key_c1d
686	CALL eos ( pts, rhd, rhop, gdept_0(:,:,:) )
687	CALL eos_rab( pts, rab_n ) ! now local thermal/haline expension ratio at T-points
688	CALL bn2 ( pts, rab_n, rn2 ) ! now Brunt-Vaisala
689
690	! Partial steps: before Horizontal DErivative
691	IF( ln_zps .AND. .NOT. ln_isfcav) &
692	& CALL zps_hde ( kt, jpts, pts, gtsu, gtsv, & ! Partial steps: before horizontal gradient
693	& rhd, gru , grv ) ! of t, s, rd at the last ocean level
694	IF( ln_zps .AND. ln_isfcav) &
695	& CALL zps_hde_isf( kt, jpts, pts, gtsu, gtsv, & ! Partial steps for top cell (ISF)
696	& rhd, gru , grv , aru , arv , gzu , gzv , ge3ru , ge3rv , &
697	& gtui, gtvi, grui, grvi, arui, arvi, gzui, gzvi, ge3rui, ge3rvi ) ! of t, s, rd at the first ocean level
698
699	rn2b(:,:,:) = rn2(:,:,:) ! need for zdfmxl
700	CALL zdf_mxl( kt ) ! mixed layer depth
701	CALL ldf_slp( kt, rhd, rn2 ) ! slopes
702	puslp (:,:,:) = uslp (:,:,:)
703	pvslp (:,:,:) = vslp (:,:,:)
704	pwslpi(:,:,:) = wslpi(:,:,:)
705	pwslpj(:,:,:) = wslpj(:,:,:)
706	#else
707	puslp (:,:,:) = 0. ! to avoid warning when compiling
708	pvslp (:,:,:) = 0.
709	pwslpi(:,:,:) = 0.
710	pwslpj(:,:,:) = 0.
711	#endif
712	!
713	END SUBROUTINE compute_slopes
714
715	!!======================================================================
716	END MODULE dtadyn

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

Download in other formats: