Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

dtadyn.F90 @ 7910

Last change on this file since 7910 was 7910, checked in by timgraham, 7 years ago
All wrk_alloc removed
Property svn:keywords set to `Id`
File size: 35.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	!!----------------------------------------------------------------------
16
17	!!----------------------------------------------------------------------
18	!! dta_dyn_init : initialization, namelist read, and SAVEs control
19	!! dta_dyn : Interpolation of the fields
20	!!----------------------------------------------------------------------
21	USE oce ! ocean dynamics and tracers variables
22	USE c1d ! 1D configuration: lk_c1d
23	USE dom_oce ! ocean domain: variables
24	USE domvvl ! variable volume
25	USE zdf_oce ! ocean vertical physics: variables
26	USE sbc_oce ! surface module: variables
27	USE trc_oce ! share ocean/biogeo variables
28	USE phycst ! physical constants
29	USE trabbl ! active tracer: bottom boundary layer
30	USE ldfslp ! lateral diffusion: iso-neutral slopes
31	USE sbcrnf ! river runoffs
32	USE ldftra ! ocean tracer lateral physics
33	USE zdfmxl ! vertical physics: mixed layer depth
34	USE eosbn2 ! equation of state - Brunt Vaisala frequency
35	USE lbclnk ! ocean lateral boundary conditions (or mpp link)
36	USE zpshde ! z-coord. with partial steps: horizontal derivatives
37	USE in_out_manager ! I/O manager
38	USE iom ! I/O library
39	USE lib_mpp ! distributed memory computing library
40	USE prtctl ! print control
41	USE fldread ! read input fields
42	USE timing ! Timing
43	USE trc, ONLY : ln_rsttr, numrtr, numrtw, lrst_trc
44
45	IMPLICIT NONE
46	PRIVATE
47
48	PUBLIC dta_dyn_init ! called by opa.F90
49	PUBLIC dta_dyn ! called by step.F90
50	PUBLIC dta_dyn_swp ! called by step.F90
51
52	CHARACTER(len=100) :: cn_dir !: Root directory for location of ssr files
53	LOGICAL :: ln_dynrnf !: read runoff data in file (T) or set to zero (F)
54	LOGICAL :: ln_dynrnf_depth !: read runoff data in file (T) or set to zero (F)
55	REAL(wp) :: fwbcorr
56
57
58	INTEGER , PARAMETER :: jpfld = 20 ! maximum number of fields to read
59	INTEGER , SAVE :: jf_tem ! index of temperature
60	INTEGER , SAVE :: jf_sal ! index of salinity
61	INTEGER , SAVE :: jf_uwd ! index of u-transport
62	INTEGER , SAVE :: jf_vwd ! index of v-transport
63	INTEGER , SAVE :: jf_wwd ! index of v-transport
64	INTEGER , SAVE :: jf_avt ! index of Kz
65	INTEGER , SAVE :: jf_mld ! index of mixed layer deptht
66	INTEGER , SAVE :: jf_emp ! index of water flux
67	INTEGER , SAVE :: jf_empb ! index of water flux
68	INTEGER , SAVE :: jf_qsr ! index of solar radiation
69	INTEGER , SAVE :: jf_wnd ! index of wind speed
70	INTEGER , SAVE :: jf_ice ! index of sea ice cover
71	INTEGER , SAVE :: jf_rnf ! index of river runoff
72	INTEGER , SAVE :: jf_fmf ! index of downward salt flux
73	INTEGER , SAVE :: jf_ubl ! index of u-bbl coef
74	INTEGER , SAVE :: jf_vbl ! index of v-bbl coef
75	INTEGER , SAVE :: jf_div ! index of e3t
76
77
78	TYPE(FLD), ALLOCATABLE, SAVE, DIMENSION(:) :: sf_dyn ! structure of input fields (file informations, fields read)
79	! !
80	REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: uslpdta ! zonal isopycnal slopes
81	REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: vslpdta ! meridional isopycnal slopes
82	REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: wslpidta ! zonal diapycnal slopes
83	REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: wslpjdta ! meridional diapycnal slopes
84
85	INTEGER, SAVE :: nprevrec, nsecdyn
86
87
88	!!----------------------------------------------------------------------
89	!! NEMO/OFF 3.3 , NEMO Consortium (2010)
90	!! $Id$
91	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
92	!!----------------------------------------------------------------------
93	CONTAINS
94
95	SUBROUTINE dta_dyn( kt )
96	!!----------------------------------------------------------------------
97	!! * ROUTINE dta_dyn *
98	!!
99	!! ** Purpose : Prepares dynamics and physics fields from a NEMO run
100	!! for an off-line simulation of passive tracers
101	!!
102	!! ** Method : calculates the position of data
103	!! - computes slopes if needed
104	!! - interpolates data if needed
105	!!----------------------------------------------------------------------
106	!
107	USE oce, ONLY: zhdivtr => ua
108	INTEGER, INTENT(in) :: kt ! ocean time-step index
109	INTEGER :: ji, jj, jk
110	REAL(wp), DIMENSION(jpi,jpj) :: zemp
111	!
112	!!----------------------------------------------------------------------
113
114	!
115	IF( nn_timing == 1 ) CALL timing_start( 'dta_dyn')
116	!
117	nsecdyn = nsec_year + nsec1jan000 ! number of seconds between Jan. 1st 00h of nit000 year and the middle of time step
118	!
119	IF( kt == nit000 ) THEN ; nprevrec = 0
120	ELSE ; nprevrec = sf_dyn(jf_tem)%nrec_a(2)
121	ENDIF
122	CALL fld_read( kt, 1, sf_dyn ) != read data at kt time step ==!
123	!
124	IF( l_ldfslp .AND. .NOT.lk_c1d ) CALL dta_dyn_slp( kt ) ! Computation of slopes
125	!
126	tsn(:,:,:,jp_tem) = sf_dyn(jf_tem)%fnow(:,:,:) * tmask(:,:,:) ! temperature
127	tsn(:,:,:,jp_sal) = sf_dyn(jf_sal)%fnow(:,:,:) * tmask(:,:,:) ! salinity
128	wndm(:,:) = sf_dyn(jf_wnd)%fnow(:,:,1) * tmask(:,:,1) ! wind speed - needed for gas exchange
129	fmmflx(:,:) = sf_dyn(jf_fmf)%fnow(:,:,1) * tmask(:,:,1) ! downward salt flux (v3.5+)
130	fr_i(:,:) = sf_dyn(jf_ice)%fnow(:,:,1) * tmask(:,:,1) ! Sea-ice fraction
131	qsr (:,:) = sf_dyn(jf_qsr)%fnow(:,:,1) * tmask(:,:,1) ! solar radiation
132	emp (:,:) = sf_dyn(jf_emp)%fnow(:,:,1) * tmask(:,:,1) ! E-P
133	IF( ln_dynrnf ) THEN
134	rnf (:,:) = sf_dyn(jf_rnf)%fnow(:,:,1) * tmask(:,:,1) ! E-P
135	IF( ln_dynrnf_depth .AND. .NOT. ln_linssh ) CALL dta_dyn_hrnf
136	ENDIF
137	!
138	un(:,:,:) = sf_dyn(jf_uwd)%fnow(:,:,:) * umask(:,:,:) ! effective u-transport
139	vn(:,:,:) = sf_dyn(jf_vwd)%fnow(:,:,:) * vmask(:,:,:) ! effective v-transport
140	wn(:,:,:) = sf_dyn(jf_wwd)%fnow(:,:,:) * tmask(:,:,:) ! effective v-transport
141	!
142	IF( .NOT.ln_linssh ) THEN
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	! Write in the tracer restart file
148	! *******************************
149	IF( lrst_trc ) THEN
150	IF(lwp) WRITE(numout,*)
151	IF(lwp) WRITE(numout,*) 'dta_dyn_ssh : ssh field written in tracer restart file ', &
152	& '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	!
169	#if defined key_trabbl && ! defined key_c1d
170	ahu_bbl(:,:) = sf_dyn(jf_ubl)%fnow(:,:,1) * umask(:,:,1) ! bbl diffusive coef
171	ahv_bbl(:,:) = sf_dyn(jf_vbl)%fnow(:,:,1) * vmask(:,:,1)
172	#endif
173	!
174	!
175	CALL eos( tsn, rhd, rhop, gdept_0(:,:,:) ) ! In any case, we need rhop
176	!
177	IF(ln_ctl) THEN ! print control
178	CALL prt_ctl(tab3d_1=tsn(:,:,:,jp_tem), clinfo1=' tn - : ', mask1=tmask, ovlap=1, kdim=jpk )
179	CALL prt_ctl(tab3d_1=tsn(:,:,:,jp_sal), clinfo1=' sn - : ', mask1=tmask, ovlap=1, kdim=jpk )
180	CALL prt_ctl(tab3d_1=un , clinfo1=' un - : ', mask1=umask, ovlap=1, kdim=jpk )
181	CALL prt_ctl(tab3d_1=vn , clinfo1=' vn - : ', mask1=vmask, ovlap=1, kdim=jpk )
182	CALL prt_ctl(tab3d_1=wn , clinfo1=' wn - : ', mask1=tmask, ovlap=1, kdim=jpk )
183	CALL prt_ctl(tab3d_1=avt , clinfo1=' kz - : ', mask1=tmask, ovlap=1, kdim=jpk )
184	CALL prt_ctl(tab3d_1=uslp , clinfo1=' slp - u : ', tab3d_2=vslp, clinfo2=' v : ', kdim=jpk)
185	CALL prt_ctl(tab3d_1=wslpi , clinfo1=' slp - wi: ', tab3d_2=wslpj, clinfo2=' wj: ', kdim=jpk)
186	! CALL prt_ctl(tab2d_1=fr_i , clinfo1=' fr_i - : ', mask1=tmask, ovlap=1 )
187	! CALL prt_ctl(tab2d_1=hmld , clinfo1=' hmld - : ', mask1=tmask, ovlap=1 )
188	! CALL prt_ctl(tab2d_1=fmmflx , clinfo1=' fmmflx - : ', mask1=tmask, ovlap=1 )
189	! CALL prt_ctl(tab2d_1=emp , clinfo1=' emp - : ', mask1=tmask, ovlap=1 )
190	! CALL prt_ctl(tab2d_1=wndm , clinfo1=' wspd - : ', mask1=tmask, ovlap=1 )
191	! CALL prt_ctl(tab2d_1=qsr , clinfo1=' qsr - : ', mask1=tmask, ovlap=1 )
192	ENDIF
193	!
194	IF( nn_timing == 1 ) CALL timing_stop( 'dta_dyn')
195	!
196	END SUBROUTINE dta_dyn
197
198
199	SUBROUTINE dta_dyn_init
200	!!----------------------------------------------------------------------
201	!! * ROUTINE dta_dyn_init *
202	!!
203	!! ** Purpose : Initialisation of the dynamical data
204	!! ** Method : - read the data namdta_dyn namelist
205	!!
206	!! ** Action : - read parameters
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	!!----------------------------------------------------------------------
227	!
228	NAMELIST/namdta_dyn/cn_dir, ln_dynrnf, ln_dynrnf_depth, fwbcorr, &
229	& sn_uwd, sn_vwd, sn_wwd, sn_emp, &
230	& sn_avt, sn_tem, sn_sal, sn_mld , sn_qsr , &
231	& sn_wnd, sn_ice, sn_fmf, &
232	& sn_ubl, sn_vbl, sn_rnf, &
233	& sn_empb, sn_div
234	!
235	REWIND( numnam_ref ) ! Namelist namdta_dyn in reference namelist : Offline: init. of dynamical data
236	READ ( numnam_ref, namdta_dyn, IOSTAT = ios, ERR = 901)
237	901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdta_dyn in reference namelist', lwp )
238
239	REWIND( numnam_cfg ) ! Namelist namdta_dyn in configuration namelist : Offline: init. of dynamical data
240	READ ( numnam_cfg, namdta_dyn, IOSTAT = ios, ERR = 902 )
241	902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdta_dyn in configuration namelist', lwp )
242	IF(lwm) WRITE ( numond, namdta_dyn )
243	! ! store namelist information in an array
244	! ! Control print
245	IF(lwp) THEN
246	WRITE(numout,*)
247	WRITE(numout,*) 'dta_dyn : offline dynamics '
248	WRITE(numout,*) '~~~~~~~ '
249	WRITE(numout,*) ' Namelist namdta_dyn'
250	WRITE(numout,*) ' runoffs option enabled (T) or not (F) ln_dynrnf = ', ln_dynrnf
251	WRITE(numout,*) ' runoffs is spread in vertical ln_dynrnf_depth = ', ln_dynrnf_depth
252	WRITE(numout,*) ' annual global mean of empmr for ssh correction fwbcorr = ', fwbcorr
253	WRITE(numout,*)
254	ENDIF
255	!
256
257	jf_uwd = 1 ; jf_vwd = 2 ; jf_wwd = 3 ; jf_emp = 4 ; jf_avt = 5
258	jf_tem = 6 ; jf_sal = 7 ; jf_mld = 8 ; jf_qsr = 9
259	jf_wnd = 10 ; jf_ice = 11 ; jf_fmf = 12 ; jfld = jf_fmf
260
261	!
262	slf_d(jf_uwd) = sn_uwd ; slf_d(jf_vwd) = sn_vwd ; slf_d(jf_wwd) = sn_wwd
263	slf_d(jf_emp) = sn_emp ; slf_d(jf_avt) = sn_avt
264	slf_d(jf_tem) = sn_tem ; slf_d(jf_sal) = sn_sal ; slf_d(jf_mld) = sn_mld
265	slf_d(jf_qsr) = sn_qsr ; slf_d(jf_wnd) = sn_wnd ; slf_d(jf_ice) = sn_ice
266	slf_d(jf_fmf) = sn_fmf
267
268	!
269	IF( .NOT.ln_linssh ) THEN
270	jf_div = jfld + 1 ; jf_empb = jfld + 2 ; jfld = jf_empb
271	slf_d(jf_div) = sn_div ; slf_d(jf_empb) = sn_empb
272	ENDIF
273	!
274	IF( lk_trabbl ) THEN
275	jf_ubl = jfld + 1 ; jf_vbl = jfld + 2 ; jfld = jf_vbl
276	slf_d(jf_ubl) = sn_ubl ; slf_d(jf_vbl) = sn_vbl
277	ENDIF
278	!
279	IF( ln_dynrnf ) THEN
280	jf_rnf = jfld + 1 ; jfld = jf_rnf
281	slf_d(jf_rnf) = sn_rnf
282	ELSE
283	rnf(:,:) = 0._wp
284	ENDIF
285
286
287	ALLOCATE( sf_dyn(jfld), STAT=ierr ) ! set sf structure
288	IF( ierr > 0 ) THEN
289	CALL ctl_stop( 'dta_dyn: unable to allocate sf structure' ) ; RETURN
290	ENDIF
291	! ! fill sf with slf_i and control print
292	CALL fld_fill( sf_dyn, slf_d, cn_dir, 'dta_dyn_init', 'Data in file', 'namdta_dyn' )
293	!
294	! Open file for each variable to get his number of dimension
295	DO ifpr = 1, jfld
296	CALL fld_clopn( sf_dyn(ifpr), nyear, nmonth, nday )
297	idv = iom_varid( sf_dyn(ifpr)%num , slf_d(ifpr)%clvar ) ! id of the variable sdjf%clvar
298	idimv = iom_file ( sf_dyn(ifpr)%num )%ndims(idv) ! number of dimension for variable sdjf%clvar
299	IF( sf_dyn(ifpr)%num /= 0 ) CALL iom_close( sf_dyn(ifpr)%num ) ! close file if already open
300	ierr1=0
301	IF( idimv == 3 ) THEN ! 2D variable
302	ALLOCATE( sf_dyn(ifpr)%fnow(jpi,jpj,1) , STAT=ierr0 )
303	IF( slf_d(ifpr)%ln_tint ) ALLOCATE( sf_dyn(ifpr)%fdta(jpi,jpj,1,2) , STAT=ierr1 )
304	ELSE ! 3D variable
305	ALLOCATE( sf_dyn(ifpr)%fnow(jpi,jpj,jpk) , STAT=ierr0 )
306	IF( slf_d(ifpr)%ln_tint ) ALLOCATE( sf_dyn(ifpr)%fdta(jpi,jpj,jpk,2), STAT=ierr1 )
307	ENDIF
308	IF( ierr0 + ierr1 > 0 ) THEN
309	CALL ctl_stop( 'dta_dyn_init : unable to allocate sf_dyn array structure' ) ; RETURN
310	ENDIF
311	END DO
312	!
313	IF( l_ldfslp .AND. .NOT.lk_c1d ) THEN ! slopes
314	IF( sf_dyn(jf_tem)%ln_tint ) THEN ! time interpolation
315	ALLOCATE( uslpdta (jpi,jpj,jpk,2), vslpdta (jpi,jpj,jpk,2), &
316	& wslpidta(jpi,jpj,jpk,2), wslpjdta(jpi,jpj,jpk,2), STAT=ierr2 )
317	!
318	IF( ierr2 > 0 ) THEN
319	CALL ctl_stop( 'dta_dyn_init : unable to allocate slope arrays' ) ; RETURN
320	ENDIF
321	ENDIF
322	ENDIF
323	!
324	IF( .NOT.ln_linssh ) THEN
325	IF( .NOT. sf_dyn(jf_uwd)%ln_clim .AND. ln_rsttr .AND. & ! Restart: read in restart file
326	iom_varid( numrtr, 'sshn', ldstop = .FALSE. ) > 0 ) THEN
327	IF(lwp) WRITE(numout,*) ' sshn forcing fields read in the restart file for initialisation'
328	CALL iom_get( numrtr, jpdom_autoglo, 'sshn', sshn(:,:) )
329	CALL iom_get( numrtr, jpdom_autoglo, 'sshb', sshb(:,:) )
330	ELSE
331	IF(lwp) WRITE(numout,*) ' sshn forcing fields read in the restart file for initialisation'
332	CALL iom_open( 'restart', inum )
333	CALL iom_get( inum, jpdom_autoglo, 'sshn', sshn(:,:) )
334	CALL iom_get( inum, jpdom_autoglo, 'sshb', sshb(:,:) )
335	CALL iom_close( inum ) ! close file
336	ENDIF
337	!
338	DO jk = 1, jpkm1
339	e3t_n(:,:,jk) = e3t_0(:,:,jk) * ( 1._wp + sshn(:,:) * tmask(:,:,1) / ( ht_0(:,:) + 1.0 - tmask(:,:,1) ) )
340	ENDDO
341	e3t_a(:,:,jpk) = e3t_0(:,:,jpk)
342
343	! Horizontal scale factor interpolations
344	! --------------------------------------
345	CALL dom_vvl_interpol( e3t_n(:,:,:), e3u_n(:,:,:), 'U' )
346	CALL dom_vvl_interpol( e3t_n(:,:,:), e3v_n(:,:,:), 'V' )
347
348	! Vertical scale factor interpolations
349	! ------------------------------------
350	CALL dom_vvl_interpol( e3t_n(:,:,:), e3w_n(:,:,:), 'W' )
351
352	e3t_b(:,:,:) = e3t_n(:,:,:)
353	e3u_b(:,:,:) = e3u_n(:,:,:)
354	e3v_b(:,:,:) = e3v_n(:,:,:)
355
356	! t- and w- points depth
357	! ----------------------
358	gdept_n(:,:,1) = 0.5_wp * e3w_n(:,:,1)
359	gdepw_n(:,:,1) = 0.0_wp
360
361	DO jk = 2, jpk
362	DO jj = 1,jpj
363	DO ji = 1,jpi
364	! zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk)) ! 0 everywhere
365	! tmask = wmask, ie everywhere expect at jk = mikt
366	! 1 for jk =
367	! mikt
368	zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk))
369	gdepw_n(ji,jj,jk) = gdepw_n(ji,jj,jk-1) + e3t_n(ji,jj,jk-1)
370	gdept_n(ji,jj,jk) = zcoef * ( gdepw_n(ji,jj,jk ) + 0.5 * e3w_n(ji,jj,jk)) &
371	& + (1-zcoef) * ( gdept_n(ji,jj,jk-1) + e3w_n(ji,jj,jk))
372	END DO
373	END DO
374	END DO
375
376	gdept_b(:,:,:) = gdept_n(:,:,:)
377	gdepw_b(:,:,:) = gdepw_n(:,:,:)
378	!
379	ENDIF
380	!
381	IF( ln_dynrnf .AND. ln_dynrnf_depth ) THEN ! read depht over which runoffs are distributed
382	IF(lwp) WRITE(numout,*)
383	IF(lwp) WRITE(numout,*) ' read in the file depht over which runoffs are distributed'
384	CALL iom_open ( "runoffs", inum ) ! open file
385	CALL iom_get ( inum, jpdom_data, 'rodepth', h_rnf ) ! read the river mouth array
386	CALL iom_close( inum ) ! close file
387	!
388	nk_rnf(:,:) = 0 ! set the number of level over which river runoffs are applied
389	DO jj = 1, jpj
390	DO ji = 1, jpi
391	IF( h_rnf(ji,jj) > 0._wp ) THEN
392	jk = 2
393	DO WHILE ( jk /= mbkt(ji,jj) .AND. gdept_0(ji,jj,jk) < h_rnf(ji,jj) ) ; jk = jk + 1
394	END DO
395	nk_rnf(ji,jj) = jk
396	ELSEIF( h_rnf(ji,jj) == -1._wp ) THEN ; nk_rnf(ji,jj) = 1
397	ELSEIF( h_rnf(ji,jj) == -999._wp ) THEN ; nk_rnf(ji,jj) = mbkt(ji,jj)
398	ELSE
399	CALL ctl_stop( 'sbc_rnf_init: runoff depth not positive, and not -999 or -1, rnf value in file fort.999' )
400	WRITE(999,*) 'ji, jj, h_rnf(ji,jj) :', ji, jj, h_rnf(ji,jj)
401	ENDIF
402	END DO
403	END DO
404	DO jj = 1, jpj ! set the associated depth
405	DO ji = 1, jpi
406	h_rnf(ji,jj) = 0._wp
407	DO jk = 1, nk_rnf(ji,jj)
408	h_rnf(ji,jj) = h_rnf(ji,jj) + e3t_n(ji,jj,jk)
409	END DO
410	END DO
411	END DO
412	ELSE ! runoffs applied at the surface
413	nk_rnf(:,:) = 1
414	h_rnf (:,:) = e3t_n(:,:,1)
415	ENDIF
416	nkrnf_max = MAXVAL( nk_rnf(:,:) )
417	hrnf_max = MAXVAL( h_rnf(:,:) )
418	IF( lk_mpp ) THEN
419	CALL mpp_max( nkrnf_max ) ! max over the global domain
420	CALL mpp_max( hrnf_max ) ! max over the global domain
421	ENDIF
422	IF(lwp) WRITE(numout,*) ' '
423	IF(lwp) WRITE(numout,*) ' max depht of runoff : ', hrnf_max,' max level : ', nkrnf_max
424	IF(lwp) WRITE(numout,*) ' '
425	!
426	CALL dta_dyn( nit000 )
427	!
428	END SUBROUTINE dta_dyn_init
429
430	SUBROUTINE dta_dyn_swp( kt )
431	!!---------------------------------------------------------------------
432	!! * ROUTINE dta_dyn_swp *
433	!!
434	!! ** Purpose : Swap and the data and compute the vertical scale factor at U/V/W point
435	!! and the depht
436	!!
437	!!---------------------------------------------------------------------
438	INTEGER, INTENT(in) :: kt ! time step
439	INTEGER :: ji, jj, jk
440	REAL(wp) :: zcoef
441	!
442	!!---------------------------------------------------------------------
443
444	IF( kt == nit000 ) THEN
445	IF(lwp) WRITE(numout,*)
446	IF(lwp) WRITE(numout,*) 'ssh_swp : Asselin time filter and swap of sea surface height'
447	IF(lwp) WRITE(numout,*) '~~~~~~~ '
448	ENDIF
449
450	sshb(:,:) = sshn(:,:) + atfp * ( sshb(:,:) - 2 * sshn(:,:) + ssha(:,:)) ! before <-- now filtered
451	sshn(:,:) = ssha(:,:)
452
453	e3t_n(:,:,:) = e3t_a(:,:,:)
454
455	! Reconstruction of all vertical scale factors at now and before time steps
456	! =============================================================================
457
458	! Horizontal scale factor interpolations
459	! --------------------------------------
460	CALL dom_vvl_interpol( e3t_n(:,:,:), e3u_n(:,:,:), 'U' )
461	CALL dom_vvl_interpol( e3t_n(:,:,:), e3v_n(:,:,:), 'V' )
462
463	! Vertical scale factor interpolations
464	! ------------------------------------
465	CALL dom_vvl_interpol( e3t_n(:,:,:), e3w_n (:,:,:), 'W' )
466
467	e3t_b(:,:,:) = e3t_n(:,:,:)
468	e3u_b(:,:,:) = e3u_n(:,:,:)
469	e3v_b(:,:,:) = e3v_n(:,:,:)
470
471	! t- and w- points depth
472	! ----------------------
473	gdept_n(:,:,1) = 0.5_wp * e3w_n(:,:,1)
474	gdepw_n(:,:,1) = 0.0_wp
475
476	DO jk = 2, jpk
477	DO jj = 1,jpj
478	DO ji = 1,jpi
479	zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk))
480	gdepw_n(ji,jj,jk) = gdepw_n(ji,jj,jk-1) + e3t_n(ji,jj,jk-1)
481	gdept_n(ji,jj,jk) = zcoef * ( gdepw_n(ji,jj,jk ) + 0.5 * e3w_n(ji,jj,jk)) &
482	& + (1-zcoef) * ( gdept_n(ji,jj,jk-1) + e3w_n(ji,jj,jk))
483	END DO
484	END DO
485	END DO
486
487	gdept_b(:,:,:) = gdept_n(:,:,:)
488	gdepw_b(:,:,:) = gdepw_n(:,:,:)
489
490	!
491	END SUBROUTINE dta_dyn_swp
492
493	SUBROUTINE dta_dyn_ssh( kt, phdivtr, psshb, pemp, pssha, pe3ta )
494	!!----------------------------------------------------------------------
495	!! * ROUTINE dta_dyn_wzv *
496	!!
497	!! ** Purpose : compute the after ssh (ssha) and the now vertical velocity
498	!!
499	!! ** Method : Using the incompressibility hypothesis,
500	!! - the ssh increment is computed by integrating the horizontal divergence
501	!! and multiply by the time step.
502	!!
503	!! - compute the after scale factor : repartition of ssh INCREMENT proportionnaly
504	!! to the level thickness ( z-star case )
505	!!
506	!! - the vertical velocity is computed by integrating the horizontal divergence
507	!! from the bottom to the surface minus the scale factor evolution.
508	!! The boundary conditions are w=0 at the bottom (no flux)
509	!!
510	!! ** action : ssha / e3t_a / wn
511	!!
512	!! Reference : Leclair, M., and G. Madec, 2009, Ocean Modelling.
513	!!----------------------------------------------------------------------
514	!! * Arguments
515	INTEGER, INTENT(in ) :: kt ! time-step
516	REAL(wp), DIMENSION(jpi,jpj,jpk) , INTENT(in ) :: phdivtr ! horizontal divergence transport
517	REAL(wp), DIMENSION(jpi,jpj) , OPTIONAL, INTENT(in ) :: psshb ! now ssh
518	REAL(wp), DIMENSION(jpi,jpj) , OPTIONAL, INTENT(in ) :: pemp ! evaporation minus precipitation
519	REAL(wp), DIMENSION(jpi,jpj) , OPTIONAL, INTENT(inout) :: pssha ! after ssh
520	REAL(wp), DIMENSION(jpi,jpj,jpk), OPTIONAL, INTENT(out) :: pe3ta ! after vertical scale factor
521	!! * Local declarations
522	INTEGER :: jk
523	REAL(wp), DIMENSION(jpi,jpj) :: zhdiv
524	REAL(wp) :: z2dt
525	!!----------------------------------------------------------------------
526
527	!
528	z2dt = 2._wp * rdt
529	!
530	zhdiv(:,:) = 0._wp
531	DO jk = 1, jpkm1
532	zhdiv(:,:) = zhdiv(:,:) + phdivtr(:,:,jk) * tmask(:,:,jk)
533	END DO
534	! ! Sea surface elevation time-stepping
535	pssha(:,:) = ( psshb(:,:) - z2dt * ( r1_rau0 * pemp(:,:) + zhdiv(:,:) ) ) * ssmask(:,:)
536	! !
537	! ! After acale factors at t-points ( z_star coordinate )
538	DO jk = 1, jpkm1
539	pe3ta(:,:,jk) = e3t_0(:,:,jk) * ( 1._wp + pssha(:,:) * tmask(:,:,1) / ( ht_0(:,:) + 1.0 - tmask(:,:,1) ) )
540	END DO
541	!
542	END SUBROUTINE dta_dyn_ssh
543
544
545	SUBROUTINE dta_dyn_hrnf
546	!!----------------------------------------------------------------------
547	!! * ROUTINE sbc_rnf *
548	!!
549	!! ** Purpose : update the horizontal divergence with the runoff inflow
550	!!
551	!! ** Method :
552	!! CAUTION : rnf is positive (inflow) decreasing the
553	!! divergence and expressed in m/s
554	!!
555	!! ** Action : phdivn decreased by the runoff inflow
556	!!----------------------------------------------------------------------
557	!!
558	INTEGER :: ji, jj, jk ! dummy loop indices
559	!!----------------------------------------------------------------------
560	!
561	DO jj = 1, jpj ! update the depth over which runoffs are distributed
562	DO ji = 1, jpi
563	h_rnf(ji,jj) = 0._wp
564	DO jk = 1, nk_rnf(ji,jj) ! recalculates h_rnf to be the depth in metres
565	h_rnf(ji,jj) = h_rnf(ji,jj) + e3t_n(ji,jj,jk) ! to the bottom of the relevant grid box
566	END DO
567	END DO
568	END DO
569	!
570	END SUBROUTINE dta_dyn_hrnf
571
572
573
574	SUBROUTINE dta_dyn_slp( kt )
575	!!---------------------------------------------------------------------
576	!! * ROUTINE dta_dyn_slp *
577	!!
578	!! ** Purpose : Computation of slope
579	!!
580	!!---------------------------------------------------------------------
581	USE oce, ONLY: zts => tsa
582	!
583	INTEGER, INTENT(in) :: kt ! time step
584	!
585	INTEGER :: ji, jj ! dummy loop indices
586	REAL(wp) :: ztinta ! ratio applied to after records when doing time interpolation
587	REAL(wp) :: ztintb ! ratio applied to before records when doing time interpolation
588	INTEGER :: iswap
589	REAL(wp), DIMENSION(jpi,jpj,jpk) :: zuslp, zvslp, zwslpi, zwslpj
590	!!---------------------------------------------------------------------
591	!
592	!
593	IF( sf_dyn(jf_tem)%ln_tint ) THEN ! Computes slopes (here avt is used as workspace)
594	IF( kt == nit000 ) THEN
595	IF(lwp) WRITE(numout,*) ' Compute new slopes at kt = ', kt
596	zts(:,:,:,jp_tem) = sf_dyn(jf_tem)%fdta(:,:,:,1) * tmask(:,:,:) ! temperature
597	zts(:,:,:,jp_sal) = sf_dyn(jf_sal)%fdta(:,:,:,1) * tmask(:,:,:) ! salinity
598	avt(:,:,:) = sf_dyn(jf_avt)%fdta(:,:,:,1) * tmask(:,:,:) ! vertical diffusive coef.
599	CALL compute_slopes( kt, zts, zuslp, zvslp, zwslpi, zwslpj )
600	uslpdta (:,:,:,1) = zuslp (:,:,:)
601	vslpdta (:,:,:,1) = zvslp (:,:,:)
602	wslpidta(:,:,:,1) = zwslpi(:,:,:)
603	wslpjdta(:,:,:,1) = zwslpj(:,:,:)
604	!
605	zts(:,:,:,jp_tem) = sf_dyn(jf_tem)%fdta(:,:,:,2) * tmask(:,:,:) ! temperature
606	zts(:,:,:,jp_sal) = sf_dyn(jf_sal)%fdta(:,:,:,2) * tmask(:,:,:) ! salinity
607	avt(:,:,:) = sf_dyn(jf_avt)%fdta(:,:,:,2) * tmask(:,:,:) ! vertical diffusive coef.
608	CALL compute_slopes( kt, zts, zuslp, zvslp, zwslpi, zwslpj )
609	uslpdta (:,:,:,2) = zuslp (:,:,:)
610	vslpdta (:,:,:,2) = zvslp (:,:,:)
611	wslpidta(:,:,:,2) = zwslpi(:,:,:)
612	wslpjdta(:,:,:,2) = zwslpj(:,:,:)
613	ELSE
614	!
615	iswap = 0
616	IF( sf_dyn(jf_tem)%nrec_a(2) - nprevrec /= 0 ) iswap = 1
617	IF( nsecdyn > sf_dyn(jf_tem)%nrec_b(2) .AND. iswap == 1 ) THEN ! read/update the after data
618	IF(lwp) WRITE(numout,*) ' Compute new slopes at kt = ', kt
619	uslpdta (:,:,:,1) = uslpdta (:,:,:,2) ! swap the data
620	vslpdta (:,:,:,1) = vslpdta (:,:,:,2)
621	wslpidta(:,:,:,1) = wslpidta(:,:,:,2)
622	wslpjdta(:,:,:,1) = wslpjdta(:,:,:,2)
623	!
624	zts(:,:,:,jp_tem) = sf_dyn(jf_tem)%fdta(:,:,:,2) * tmask(:,:,:) ! temperature
625	zts(:,:,:,jp_sal) = sf_dyn(jf_sal)%fdta(:,:,:,2) * tmask(:,:,:) ! salinity
626	avt(:,:,:) = sf_dyn(jf_avt)%fdta(:,:,:,2) * tmask(:,:,:) ! vertical diffusive coef.
627	CALL compute_slopes( kt, zts, zuslp, zvslp, zwslpi, zwslpj )
628	!
629	uslpdta (:,:,:,2) = zuslp (:,:,:)
630	vslpdta (:,:,:,2) = zvslp (:,:,:)
631	wslpidta(:,:,:,2) = zwslpi(:,:,:)
632	wslpjdta(:,:,:,2) = zwslpj(:,:,:)
633	ENDIF
634	ENDIF
635	ENDIF
636	!
637	IF( sf_dyn(jf_tem)%ln_tint ) THEN
638	ztinta = REAL( nsecdyn - sf_dyn(jf_tem)%nrec_b(2), wp ) &
639	& / REAL( sf_dyn(jf_tem)%nrec_a(2) - sf_dyn(jf_tem)%nrec_b(2), wp )
640	ztintb = 1. - ztinta
641	IF( l_ldfslp .AND. .NOT.lk_c1d ) THEN ! Computes slopes (here avt is used as workspace)
642	uslp (:,:,:) = ztintb * uslpdta (:,:,:,1) + ztinta * uslpdta (:,:,:,2)
643	vslp (:,:,:) = ztintb * vslpdta (:,:,:,1) + ztinta * vslpdta (:,:,:,2)
644	wslpi(:,:,:) = ztintb * wslpidta(:,:,:,1) + ztinta * wslpidta(:,:,:,2)
645	wslpj(:,:,:) = ztintb * wslpjdta(:,:,:,1) + ztinta * wslpjdta(:,:,:,2)
646	ENDIF
647	ELSE
648	zts(:,:,:,jp_tem) = sf_dyn(jf_tem)%fnow(:,:,:) * tmask(:,:,:) ! temperature
649	zts(:,:,:,jp_sal) = sf_dyn(jf_sal)%fnow(:,:,:) * tmask(:,:,:) ! salinity
650	avt(:,:,:) = sf_dyn(jf_avt)%fnow(:,:,:) * tmask(:,:,:) ! vertical diffusive coef.
651	CALL compute_slopes( kt, zts, zuslp, zvslp, zwslpi, zwslpj )
652	!
653	IF( l_ldfslp .AND. .NOT.lk_c1d ) THEN ! Computes slopes (here avt is used as workspace)
654	uslp (:,:,:) = zuslp (:,:,:)
655	vslp (:,:,:) = zvslp (:,:,:)
656	wslpi(:,:,:) = zwslpi(:,:,:)
657	wslpj(:,:,:) = zwslpj(:,:,:)
658	ENDIF
659	ENDIF
660	!
661	!
662	END SUBROUTINE dta_dyn_slp
663
664	SUBROUTINE compute_slopes( kt, pts, puslp, pvslp, pwslpi, pwslpj )
665	!!---------------------------------------------------------------------
666	!! * ROUTINE dta_dyn_slp *
667	!!
668	!! ** Purpose : Computation of slope
669	!!
670	!!---------------------------------------------------------------------
671	INTEGER , INTENT(in ) :: kt ! time step
672	REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in ) :: pts ! temperature/salinity
673	REAL(wp), DIMENSION(jpi,jpj,jpk) , INTENT(out) :: puslp ! zonal isopycnal slopes
674	REAL(wp), DIMENSION(jpi,jpj,jpk) , INTENT(out) :: pvslp ! meridional isopycnal slopes
675	REAL(wp), DIMENSION(jpi,jpj,jpk) , INTENT(out) :: pwslpi ! zonal diapycnal slopes
676	REAL(wp), DIMENSION(jpi,jpj,jpk) , INTENT(out) :: pwslpj ! meridional diapycnal slopes
677	!!---------------------------------------------------------------------
678	IF( l_ldfslp .AND. .NOT.lk_c1d ) THEN ! Computes slopes (here avt is used as workspace)
679	CALL eos ( pts, rhd, rhop, gdept_0(:,:,:) )
680	CALL eos_rab( pts, rab_n ) ! now local thermal/haline expension ratio at T-points
681	CALL bn2 ( pts, rab_n, rn2 ) ! now Brunt-Vaisala
682
683	! Partial steps: before Horizontal DErivative
684	IF( ln_zps .AND. .NOT. ln_isfcav) &
685	& CALL zps_hde ( kt, jpts, pts, gtsu, gtsv, & ! Partial steps: before horizontal gradient
686	& rhd, gru , grv ) ! of t, s, rd at the last ocean level
687	IF( ln_zps .AND. ln_isfcav) &
688	& CALL zps_hde_isf( kt, jpts, pts, gtsu, gtsv, gtui, gtvi, & ! Partial steps for top cell (ISF)
689	& rhd, gru , grv , grui, grvi ) ! of t, s, rd at the first ocean level
690
691	rn2b(:,:,:) = rn2(:,:,:) ! need for zdfmxl
692	CALL zdf_mxl( kt ) ! mixed layer depth
693	CALL ldf_slp( kt, rhd, rn2 ) ! slopes
694	puslp (:,:,:) = uslp (:,:,:)
695	pvslp (:,:,:) = vslp (:,:,:)
696	pwslpi(:,:,:) = wslpi(:,:,:)
697	pwslpj(:,:,:) = wslpj(:,:,:)
698	ELSE
699	puslp (:,:,:) = 0. ! to avoid warning when compiling
700	pvslp (:,:,:) = 0.
701	pwslpi(:,:,:) = 0.
702	pwslpj(:,:,:) = 0.
703	ENDIF
704	!
705	END SUBROUTINE compute_slopes
706	!!======================================================================
707	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: branches/2017/dev_r7881_no_wrk_alloc/NEMOGCM/NEMO/OFF_SRC/dtadyn.F90 @ 7910

Download in other formats: