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.

trdmxl.F90 in branches/2015/dev_merge_2015/NEMOGCM/NEMO/OPA_SRC/TRD – NEMO

Context Navigation

source: branches/2015/dev_merge_2015/NEMOGCM/NEMO/OPA_SRC/TRD/trdmxl.F90 @ 6060

Last change on this file since 6060 was 6060, checked in by timgraham, 8 years ago
Merged dev_r5836_noc2_VVL_BY_DEFAULT into branch
Property svn:keywords set to `Id`
File size: 44.0 KB

Line
1	MODULE trdmxl
2	!!======================================================================
3	!! * MODULE trdmxl *
4	!! Ocean diagnostics: mixed layer T-S trends
5	!!======================================================================
6	!! History : OPA ! 1995-04 (J. Vialard) Original code
7	!! ! 1997-02 (E. Guilyardi) Adaptation global + base cmo
8	!! ! 1999-09 (E. Guilyardi) Re-writing + netCDF output
9	!! NEMO 1.0 ! 2002-06 (G. Madec) F90: Free form and module
10	!! - ! 2004-08 (C. Talandier) New trends organization
11	!! 2.0 ! 2005-05 (C. Deltel) Diagnose trends of time averaged ML T & S
12	!! 3.5 ! 2012-03 (G. Madec) complete reorganisation + change in the time averaging
13	!!----------------------------------------------------------------------
14
15	!!----------------------------------------------------------------------
16	!! trd_mxl : T and S cumulated trends averaged over the mixed layer
17	!! trd_mxl_zint : T and S trends vertical integration
18	!! trd_mxl_init : initialization step
19	!!----------------------------------------------------------------------
20	USE oce ! ocean dynamics and tracers variables
21	USE dom_oce ! ocean space and time domain variables
22	USE trd_oce ! trends: ocean variables
23	USE trdmxl_oce ! ocean variables trends
24	USE ldftra ! lateral diffusion: eddy diffusivity & EIV coeff.
25	USE zdf_oce ! ocean vertical physics
26	USE phycst ! Define parameters for the routines
27	USE dianam ! build the name of file (routine)
28	USE ldfslp ! iso-neutral slopes
29	USE zdfmxl ! mixed layer depth
30	USE zdfddm ! ocean vertical physics: double diffusion
31	USE lbclnk ! ocean lateral boundary conditions (or mpp link)
32	USE diadimg ! dimg direct access file format output
33	USE trdmxl_rst ! restart for diagnosing the ML trends
34	!
35	USE in_out_manager ! I/O manager
36	USE ioipsl ! NetCDF library
37	USE prtctl ! Print control
38	USE restart ! for lrst_oce
39	USE lib_mpp ! MPP library
40	USE wrk_nemo ! Memory allocation
41	USE iom
42
43	IMPLICIT NONE
44	PRIVATE
45
46	PUBLIC trd_mxl ! routine called by step.F90
47	PUBLIC trd_mxl_init ! routine called by opa.F90
48	PUBLIC trd_mxl_zint ! routine called by tracers routines
49
50	INTEGER :: nkstp ! current time step
51
52	!!gm to be moved from trdmxl_oce
53	! REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: hml ! ML depth (sum of e3t over nmln-1 levels) [m]
54	! REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: tml , sml ! now ML averaged T & S
55	! REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: tmlb_nf, smlb_nf ! not filtered before ML averaged T & S
56	!
57	!
58	! REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: hmlb, hmln ! before, now, and after Mixed Layer depths [m]
59	!
60	! REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: tb_mlb, tb_mln ! before (not filtered) tracer averaged over before and now ML
61	!
62	! REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: tn_mln ! now tracer averaged over now ML
63	!!gm end
64
65	CHARACTER (LEN=40) :: clhstnam ! name of the trends NetCDF file
66	INTEGER :: nh_t, nmoymltrd
67	INTEGER :: nidtrd
68	INTEGER, ALLOCATABLE, SAVE, DIMENSION(:) :: ndextrd1
69	INTEGER :: ndimtrd1
70	INTEGER :: ionce, icount
71
72	!! * Substitutions
73	# include "zdfddm_substitute.h90"
74	!!----------------------------------------------------------------------
75	!! NEMO/OPA 3.3 , NEMO Consortium (2010)
76	!! $Id$
77	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
78	!!----------------------------------------------------------------------
79	CONTAINS
80
81	INTEGER FUNCTION trd_mxl_alloc()
82	!!----------------------------------------------------------------------
83	!! * ROUTINE trd_mxl_alloc *
84	!!----------------------------------------------------------------------
85	ALLOCATE( ndextrd1(jpi*jpj) , STAT=trd_mxl_alloc )
86	!
87	IF( lk_mpp ) CALL mpp_sum ( trd_mxl_alloc )
88	IF( trd_mxl_alloc /= 0 ) CALL ctl_warn('trd_mxl_alloc: failed to allocate array ndextrd1')
89	END FUNCTION trd_mxl_alloc
90
91
92	SUBROUTINE trd_tra_mxl( ptrdx, ptrdy, ktrd, kt, p2dt, kmxln )
93	!!----------------------------------------------------------------------
94	!! * ROUTINE trd_tra_mng *
95	!!
96	!! ** Purpose : Dispatch all trends computation, e.g. 3D output, integral
97	!! constraints, barotropic vorticity, kinetic enrgy,
98	!! potential energy, and/or mixed layer budget.
99	!!----------------------------------------------------------------------
100	REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: ptrdx ! Temperature or U trend
101	REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: ptrdy ! Salinity or V trend
102	INTEGER , INTENT(in ) :: ktrd ! tracer trend index
103	INTEGER , INTENT(in ) :: kt ! time step index
104	REAL(wp) , INTENT(in ) :: p2dt ! time step [s]
105	REAL(wp), DIMENSION(:,:) , INTENT(in ) :: kmxln ! number of t-box for the vertical average
106	!
107	INTEGER :: ji, jj, jk ! dummy loop indices
108	!!----------------------------------------------------------------------
109
110	! !==============================!
111	IF ( kt /= nkstp ) THEN != 1st call at kt time step =!
112	! !==============================!
113	nkstp = kt
114
115
116	! !== reset trend arrays to zero ==!
117	tmltrd(:,:,:) = 0._wp ; smltrd(:,:,:) = 0._wp
118
119
120	!
121	wkx(:,:,:) = 0._wp !== now ML weights for vertical averaging ==!
122	DO jk = 1, jpktrd ! initialize wkx with vertical scale factor in mixed-layer
123	DO jj = 1,jpj
124	DO ji = 1,jpi
125	IF( jk - kmxln(ji,jj) < 0 ) wkx(ji,jj,jk) = e3t_n(ji,jj,jk) * tmask(ji,jj,jk)
126	END DO
127	END DO
128	END DO
129	hmxl(:,:) = 0._wp ! NOW mixed-layer depth
130	DO jk = 1, jpktrd
131	hmxl(:,:) = hmxl(:,:) + wkx(:,:,jk)
132	END DO
133	DO jk = 1, jpktrd ! integration weights
134	wkx(:,:,jk) = wkx(:,:,jk) / MAX( 1.e-20_wp, hmxl(:,:) ) * tmask(:,:,1)
135	END DO
136
137
138	!
139	! !== Vertically averaged T and S ==!
140	tml(:,:) = 0._wp ; sml(:,:) = 0._wp
141	DO jk = 1, jpktrd
142	tml(:,:) = tml(:,:) + wkx(:,:,jk) * tsn(:,:,jk,jp_tem)
143	sml(:,:) = sml(:,:) + wkx(:,:,jk) * tsn(:,:,jk,jp_sal)
144	END DO
145	!
146	ENDIF
147
148
149
150	! mean now trends over the now ML
151	tmltrd(:,:,ktrd) = tmltrd(:,:,ktrd) + ptrdx(:,:,jk) * wkx(:,:,jk) ! temperature
152	smltrd(:,:,ktrd) = smltrd(:,:,ktrd) + ptrdy(:,:,jk) * wkx(:,:,jk) ! salinity
153
154
155
156	!!gm to be put juste before the output !
157	! ! Lateral boundary conditions
158	! CALL lbc_lnk( tmltrd(:,:,jl), 'T', 1. )
159	! CALL lbc_lnk( smltrd(:,:,jl), 'T', 1. )
160	!!gm end
161
162
163
164	SELECT CASE( ktrd )
165	CASE( jptra_npc ) ! non-penetrative convection: regrouped with zdf
166	!!gm : to be completed !
167	! IF( ....
168	!!gm end
169	CASE( jptra_zdfp ) ! iso-neutral diffusion: "pure" vertical diffusion
170	! ! regroup iso-neutral diffusion in one term
171	tmltrd(:,:,jpmxl_ldf) = tmltrd(:,:,jpmxl_ldf) + ( tmltrd(:,:,jpmxl_zdf) - tmltrd(:,:,jpmxl_zdfp) )
172	smltrd(:,:,jpmxl_ldf) = smltrd(:,:,jpmxl_ldf) + ( smltrd(:,:,jpmxl_zdf) - smltrd(:,:,jpmxl_zdfp) )
173	! ! put in zdf the dia-neutral diffusion
174	tmltrd(:,:,jpmxl_zdf) = tmltrd(:,:,jpmxl_zdfp)
175	smltrd(:,:,jpmxl_zdf) = smltrd(:,:,jpmxl_zdfp)
176	IF( ln_zdfnpc ) THEN
177	tmltrd(:,:,jpmxl_zdf) = tmltrd(:,:,jpmxl_zdf) + tmltrd(:,:,jpmxl_npc)
178	smltrd(:,:,jpmxl_zdf) = smltrd(:,:,jpmxl_zdf) + smltrd(:,:,jpmxl_npc)
179	ENDIF
180	!
181	CASE( jptra_atf ) ! last trends of the current time step: perform the time averaging & output
182	!
183	! after ML : zhmla NB will be swaped to provide hmln and hmlb
184	!
185	! entrainement ent_1 : tb_mln - tb_mlb ==>> use previous timestep ztn_mla = tb_mln
186	! " " " tn_mln = tb_mlb (unfiltered tb!)
187	! NB: tn_mln itself comes from the 2 time step before (ta_mla)
188	!
189	! atf trend : ztbf_mln - tb_mln ==>> use previous timestep tn_mla = tb_mln
190	! need to compute tbf_mln, using the current tb
191	! which is the before fitered tracer
192	!
193	! entrainement ent_2 : zta_mla - zta_mln ==>> need to compute zta_mla and zta_mln
194	!
195	! time averaging : mean: CALL trd_mean( kt, ptrd, ptrdm )
196	! and out put the starting mean value and the total trends
197	! (i.e. difference between starting and ending values)
198	! hat : CALL trd_hat ( kt, ptrd, ptrdm )
199	! and output the starting hat value and the total hat trends
200	!
201	! swaps : hmlb <== hmln <== zhmla
202	! tb_mlb <== tn_mln <== zta_mla
203	! tb_mln <== ztn_mla ==>> now T over after h, need to be computed here
204	! to be used at next time step (unfiltered before)
205	!
206	END SELECT
207	!
208	END SUBROUTINE trd_tra_mxl
209
210
211	SUBROUTINE trd_mean( kt, ptrd, ptrdm )
212	!!----------------------------------------------------------------------
213	!! * ROUTINE trd_mean *
214	!!
215	!! ** Purpose :
216	!!----------------------------------------------------------------------
217	REAL(wp), DIMENSION(:,:,:), INTENT(in ) :: ptrd ! trend at kt
218	REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: ptrdm ! cumulative trends at kt
219	INTEGER , INTENT(in ) :: kt ! time step index
220	!!----------------------------------------------------------------------
221	!
222	IF ( kt == nn_it000 ) ptrdm(:,:,:) = 0._wp
223	!
224	ptrdm(:,:,:) = ptrdm(:,:,:) + ptrd(:,:,:)
225	!
226	IF ( MOD( kt - nn_it000 + 1, nn_trd ) == 0 ) THEN
227	!
228	! call iom put???? avec en argument le tableau de nom des trends?
229	!
230	ENDIF
231	!
232	END SUBROUTINE trd_mean
233
234
235	SUBROUTINE trd_mxl_zint( pttrdmxl, pstrdmxl, ktrd, ctype )
236	!!----------------------------------------------------------------------
237	!! * ROUTINE trd_mxl_zint *
238	!!
239	!! ** Purpose : Compute the vertical average of the 3D fields given as arguments
240	!! to the subroutine. This vertical average is performed from ocean
241	!! surface down to a chosen control surface.
242	!!
243	!! ** Method/usage :
244	!! The control surface can be either a mixed layer depth (time varying)
245	!! or a fixed surface (jk level or bowl).
246	!! Choose control surface with nn_ctls in namelist NAMTRD :
247	!! nn_ctls = 0 : use mixed layer with density criterion
248	!! nn_ctls = 1 : read index from file 'ctlsurf_idx'
249	!! nn_ctls > 1 : use fixed level surface jk = nn_ctls
250	!! Note: in the remainder of the routine, the volume between the
251	!! surface and the control surface is called "mixed-layer"
252	!!----------------------------------------------------------------------
253	INTEGER , INTENT( in ) :: ktrd ! ocean trend index
254	CHARACTER(len=2) , INTENT( in ) :: ctype ! 2D surface/bottom or 3D interior physics
255	REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( in ) :: pttrdmxl ! temperature trend
256	REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( in ) :: pstrdmxl ! salinity trend
257	!
258	INTEGER :: ji, jj, jk, isum
259	REAL(wp), POINTER, DIMENSION(:,:) :: zvlmsk
260	!!----------------------------------------------------------------------
261
262	CALL wrk_alloc( jpi, jpj, zvlmsk )
263
264	! I. Definition of control surface and associated fields
265	! ------------------------------------------------------
266	! ==> only once per time step <==
267
268	IF( icount == 1 ) THEN
269	!
270
271	!!gm BUG?
272	!!gm CAUTION: double check the definition of nmln it is the nb of w-level, not t-level I guess
273
274
275	! ... Set nmxl(ji,jj) = index of first T point below control surf. or outside mixed-layer
276	IF( nn_ctls == 0 ) THEN ! * control surface = mixed-layer with density criterion
277	nmxl(:,:) = nmln(:,:) ! array nmln computed in zdfmxl.F90
278	ELSEIF( nn_ctls == 1 ) THEN ! * control surface = read index from file
279	nmxl(:,:) = nbol(:,:)
280	ELSEIF( nn_ctls >= 2 ) THEN ! * control surface = model level
281	nn_ctls = MIN( nn_ctls, jpktrd - 1 )
282	nmxl(:,:) = nn_ctls + 1
283	ENDIF
284
285	END IF
286	!
287	CALL wrk_dealloc( jpi, jpj, zvlmsk )
288	!
289	END SUBROUTINE trd_mxl_zint
290
291
292	SUBROUTINE trd_mxl( kt, p2dt )
293	!!----------------------------------------------------------------------
294	!! * ROUTINE trd_mxl *
295	!!
296	!! ** Purpose : Compute and cumulate the mixed layer trends over an analysis
297	!! period, and write NetCDF (or dimg) outputs.
298	!!
299	!! ** Method/usage :
300	!! The stored trends can be chosen twofold (according to the ln_trdmxl_instant
301	!! logical namelist variable) :
302	!! 1) to explain the difference between initial and final
303	!! mixed-layer T & S (where initial and final relate to the
304	!! current analysis window, defined by nn_trd in the namelist)
305	!! 2) to explain the difference between the current and previous
306	!! TIME-AVERAGED mixed-layer T & S (where time-averaging is
307	!! performed over each analysis window).
308	!!
309	!! ** Consistency check :
310	!! If the control surface is fixed ( nn_ctls > 1 ), the residual term (dh/dt
311	!! entrainment) should be zero, at machine accuracy. Note that in the case
312	!! of time-averaged mixed-layer fields, this residual WILL NOT BE ZERO
313	!! over the first two analysis windows (except if restart).
314	!! N.B. For ORCA2_LIM, use e.g. nn_trd=5, rn_ucf=1., nn_ctls=8
315	!! for checking residuals.
316	!! On a NEC-SX5 computer, this typically leads to:
317	!! O(1.e-20) temp. residuals (tml_res) when ln_trdmxl_instant=.false.
318	!! O(1.e-21) temp. residuals (tml_res) when ln_trdmxl_instant=.true.
319	!!
320	!! ** Action :
321	!! At each time step, mixed-layer averaged trends are stored in the
322	!! tmltrd(:,:,jpmxl_xxx) array (see trdmxl_oce.F90 for definitions of jpmxl_xxx).
323	!! This array is known when trd_mxl is called, at the end of the stp subroutine,
324	!! except for the purely vertical K_z diffusion term, which is embedded in the
325	!! lateral diffusion trend.
326	!!
327	!! In I), this K_z term is diagnosed and stored, thus its contribution is removed
328	!! from the lateral diffusion trend.
329	!! In II), the instantaneous mixed-layer T & S are computed, and misc. cumulative
330	!! arrays are updated.
331	!! In III), called only once per analysis window, we compute the total trends,
332	!! along with the residuals and the Asselin correction terms.
333	!! In IV), the appropriate trends are written in the trends NetCDF file.
334	!!
335	!! References : Vialard et al.,2001, JPO.
336	!!----------------------------------------------------------------------
337	INTEGER , INTENT(in ) :: kt ! ocean time-step index
338	REAL(wp), INTENT(in ) :: p2dt ! time step [s]
339	!
340	INTEGER :: ji, jj, jk, jl, ik, it, itmod
341	LOGICAL :: lldebug = .TRUE.
342	REAL(wp) :: zavt, zfn, zfn2
343	! ! z(ts)mltot : dT/dt over the anlysis window (including Asselin)
344	! ! z(ts)mlres : residual = dh/dt entrainment term
345	REAL(wp), POINTER, DIMENSION(:,: ) :: ztmltot , zsmltot , ztmlres , zsmlres , ztmlatf , zsmlatf
346	REAL(wp), POINTER, DIMENSION(:,: ) :: ztmltot2, zsmltot2, ztmlres2, zsmlres2, ztmlatf2, zsmlatf2, ztmltrdm2, zsmltrdm2
347	REAL(wp), POINTER, DIMENSION(:,:,:) :: ztmltrd2, zsmltrd2 ! only needed for mean diagnostics
348	#if defined key_dimgout
349	INTEGER :: iyear,imon,iday
350	CHARACTER(LEN=80) :: cltext, clmode
351	#endif
352	!!----------------------------------------------------------------------
353
354	CALL wrk_alloc( jpi, jpj, ztmltot , zsmltot , ztmlres , zsmlres , ztmlatf , zsmlatf )
355	CALL wrk_alloc( jpi, jpj, ztmltot2, zsmltot2, ztmlres2, zsmlres2, ztmlatf2, zsmlatf2, ztmltrdm2, zsmltrdm2 )
356	CALL wrk_alloc( jpi, jpj, jpltrd, ztmltrd2, zsmltrd2 )
357
358	! ======================================================================
359	! II. Cumulate the trends over the analysis window
360	! ======================================================================
361
362	ztmltrd2(:,:,:) = 0.e0 ; zsmltrd2(:,:,:) = 0.e0 ! <<< reset arrays to zero
363	ztmltot2(:,:) = 0.e0 ; zsmltot2(:,:) = 0.e0
364	ztmlres2(:,:) = 0.e0 ; zsmlres2(:,:) = 0.e0
365	ztmlatf2(:,:) = 0.e0 ; zsmlatf2(:,:) = 0.e0
366
367	! II.1 Set before values of vertically average T and S
368	! ----------------------------------------------------
369	IF( kt > nit000 ) THEN
370	! ... temperature ... ... salinity ...
371	tmlb (:,:) = tml (:,:) ; smlb (:,:) = sml (:,:)
372	tmlatfn(:,:) = tmltrd(:,:,jpmxl_atf) ; smlatfn(:,:) = smltrd(:,:,jpmxl_atf)
373	END IF
374
375
376	! II.3 Initialize mixed-layer "before" arrays for the 1rst analysis window
377	! ------------------------------------------------------------------------
378	IF( kt == 2 ) THEN ! i.e. ( .NOT. ln_rstart ).AND.( kt == nit000 + 1)
379	!
380	! ... temperature ... ... salinity ...
381	tmlbb (:,:) = tmlb (:,:) ; smlbb (:,:) = smlb (:,:)
382	tmlbn (:,:) = tml (:,:) ; smlbn (:,:) = sml (:,:)
383	tmlatfb(:,:) = tmlatfn(:,:) ; smlatfb(:,:) = smlatfn(:,:)
384
385	tmltrd_csum_ub (:,:,:) = 0.e0 ; smltrd_csum_ub (:,:,:) = 0.e0
386	tmltrd_atf_sumb(:,:) = 0.e0 ; smltrd_atf_sumb(:,:) = 0.e0
387
388	hmxlbn(:,:) = hmxl(:,:)
389
390	IF( ln_ctl ) THEN
391	WRITE(numout,*) ' we reach kt == nit000 + 1 = ', nit000+1
392	CALL prt_ctl(tab2d_1=tmlbb , clinfo1=' tmlbb - : ', mask1=tmask, ovlap=1)
393	CALL prt_ctl(tab2d_1=tmlbn , clinfo1=' tmlbn - : ', mask1=tmask, ovlap=1)
394	CALL prt_ctl(tab2d_1=tmlatfb , clinfo1=' tmlatfb - : ', mask1=tmask, ovlap=1)
395	END IF
396	!
397	END IF
398
399	IF( ( ln_rstart ) .AND. ( kt == nit000 ) .AND. ( ln_ctl ) ) THEN
400	IF( ln_trdmxl_instant ) THEN
401	WRITE(numout,*) ' restart from kt == nit000 = ', nit000
402	CALL prt_ctl(tab2d_1=tmlbb , clinfo1=' tmlbb - : ', mask1=tmask, ovlap=1)
403	CALL prt_ctl(tab2d_1=tmlbn , clinfo1=' tmlbn - : ', mask1=tmask, ovlap=1)
404	CALL prt_ctl(tab2d_1=tmlatfb , clinfo1=' tmlatfb - : ', mask1=tmask, ovlap=1)
405	ELSE
406	WRITE(numout,*) ' restart from kt == nit000 = ', nit000
407	CALL prt_ctl(tab2d_1=tmlbn , clinfo1=' tmlbn - : ', mask1=tmask, ovlap=1)
408	CALL prt_ctl(tab2d_1=hmxlbn , clinfo1=' hmxlbn - : ', mask1=tmask, ovlap=1)
409	CALL prt_ctl(tab2d_1=tml_sumb , clinfo1=' tml_sumb - : ', mask1=tmask, ovlap=1)
410	CALL prt_ctl(tab2d_1=tmltrd_atf_sumb, clinfo1=' tmltrd_atf_sumb - : ', mask1=tmask, ovlap=1)
411	CALL prt_ctl(tab3d_1=tmltrd_csum_ub , clinfo1=' tmltrd_csum_ub - : ', mask1=tmask, ovlap=1, kdim=1)
412	END IF
413	END IF
414
415	! II.4 Cumulated trends over the analysis period
416	! ----------------------------------------------
417	!
418	! [ 1rst analysis window ] [ 2nd analysis window ]
419	!
420	! o---[--o-----o-----o-----o--]-[--o-----o-----o-----o-----o--]---o-----o--> time steps
421	! nn_trd 2*nn_trd etc.
422	! 1 2 3 4 =5 e.g. =10
423	!
424	IF( ( kt >= 2 ).OR.( ln_rstart ) ) THEN
425	!
426	nmoymltrd = nmoymltrd + 1
427
428	! ... Cumulate over BOTH physical contributions AND over time steps
429	DO jl = 1, jpltrd
430	tmltrdm(:,:) = tmltrdm(:,:) + tmltrd(:,:,jl)
431	smltrdm(:,:) = smltrdm(:,:) + smltrd(:,:,jl)
432	END DO
433
434	! ... Special handling of the Asselin trend
435	tmlatfm(:,:) = tmlatfm(:,:) + tmlatfn(:,:)
436	smlatfm(:,:) = smlatfm(:,:) + smlatfn(:,:)
437
438	! ... Trends associated with the time mean of the ML T/S
439	tmltrd_sum (:,:,:) = tmltrd_sum (:,:,:) + tmltrd (:,:,:) ! tem
440	tmltrd_csum_ln(:,:,:) = tmltrd_csum_ln(:,:,:) + tmltrd_sum(:,:,:)
441	tml_sum (:,:) = tml_sum (:,:) + tml (:,:)
442	smltrd_sum (:,:,:) = smltrd_sum (:,:,:) + smltrd (:,:,:) ! sal
443	smltrd_csum_ln(:,:,:) = smltrd_csum_ln(:,:,:) + smltrd_sum(:,:,:)
444	sml_sum (:,:) = sml_sum (:,:) + sml (:,:)
445	hmxl_sum (:,:) = hmxl_sum (:,:) + hmxl (:,:) ! rmxl
446	!
447	END IF
448
449	! ======================================================================
450	! III. Prepare fields for output (get here ONCE PER ANALYSIS PERIOD)
451	! ======================================================================
452
453	! Convert to appropriate physical units
454	! N.B. It may be useful to check IOIPSL time averaging with :
455	! tmltrd (:,:,:) = 1. ; smltrd (:,:,:) = 1.
456	tmltrd(:,:,:) = tmltrd(:,:,:) * rn_ucf ! (actually needed for 1:jpltrd-1, but trdmxl(:,:,jpltrd)
457	smltrd(:,:,:) = smltrd(:,:,:) * rn_ucf ! is no longer used, and is reset to 0. at next time step)
458
459	! define time axis
460	it = kt
461	itmod = kt - nit000 + 1
462
463	MODULO_NTRD : IF( MOD( itmod, nn_trd ) == 0 ) THEN ! nitend MUST be multiple of nn_trd
464	!
465	ztmltot (:,:) = 0.e0 ; zsmltot (:,:) = 0.e0 ! reset arrays to zero
466	ztmlres (:,:) = 0.e0 ; zsmlres (:,:) = 0.e0
467	ztmltot2(:,:) = 0.e0 ; zsmltot2(:,:) = 0.e0
468	ztmlres2(:,:) = 0.e0 ; zsmlres2(:,:) = 0.e0
469
470	zfn = REAL( nmoymltrd, wp ) ; zfn2 = zfn * zfn
471
472	! III.1 Prepare fields for output ("instantaneous" diagnostics)
473	! -------------------------------------------------------------
474
475	!-- Compute total trends
476	ztmltot(:,:) = ( tml(:,:) - tmlbn(:,:) + tmlb(:,:) - tmlbb(:,:) ) / p2dt
477	zsmltot(:,:) = ( sml(:,:) - smlbn(:,:) + smlb(:,:) - smlbb(:,:) ) / p2dt
478
479	!-- Compute residuals
480	ztmlres(:,:) = ztmltot(:,:) - ( tmltrdm(:,:) - tmlatfn(:,:) + tmlatfb(:,:) )
481	zsmlres(:,:) = zsmltot(:,:) - ( smltrdm(:,:) - smlatfn(:,:) + smlatfb(:,:) )
482
483	!-- Diagnose Asselin trend over the analysis window
484	ztmlatf(:,:) = tmlatfm(:,:) - tmlatfn(:,:) + tmlatfb(:,:)
485	zsmlatf(:,:) = smlatfm(:,:) - smlatfn(:,:) + smlatfb(:,:)
486
487	!-- Lateral boundary conditions
488	! ... temperature ... ... salinity ...
489	CALL lbc_lnk( ztmltot , 'T', 1. ) ; CALL lbc_lnk( zsmltot , 'T', 1. )
490	CALL lbc_lnk( ztmlres , 'T', 1. ) ; CALL lbc_lnk( zsmlres , 'T', 1. )
491	CALL lbc_lnk( ztmlatf , 'T', 1. ) ; CALL lbc_lnk( zsmlatf , 'T', 1. )
492
493
494	! III.2 Prepare fields for output ("mean" diagnostics)
495	! ----------------------------------------------------
496
497	!-- Update the ML depth time sum (to build the Leap-Frog time mean)
498	hmxl_sum(:,:) = hmxlbn(:,:) + 2 * ( hmxl_sum(:,:) - hmxl(:,:) ) + hmxl(:,:)
499
500	!-- Compute temperature total trends
501	tml_sum (:,:) = tmlbn(:,:) + 2 * ( tml_sum(:,:) - tml(:,:) ) + tml(:,:)
502	ztmltot2(:,:) = ( tml_sum(:,:) - tml_sumb(:,:) ) / p2dt ! now in degC/s
503
504	!-- Compute salinity total trends
505	sml_sum (:,:) = smlbn(:,:) + 2 * ( sml_sum(:,:) - sml(:,:) ) + sml(:,:)
506	zsmltot2(:,:) = ( sml_sum(:,:) - sml_sumb(:,:) ) / p2dt ! now in psu/s
507
508	!-- Compute temperature residuals
509	DO jl = 1, jpltrd
510	ztmltrd2(:,:,jl) = tmltrd_csum_ub(:,:,jl) + tmltrd_csum_ln(:,:,jl)
511	END DO
512
513	ztmltrdm2(:,:) = 0.e0
514	DO jl = 1, jpltrd
515	ztmltrdm2(:,:) = ztmltrdm2(:,:) + ztmltrd2(:,:,jl)
516	END DO
517
518	ztmlres2(:,:) = ztmltot2(:,:) - &
519	( ztmltrdm2(:,:) - tmltrd_sum(:,:,jpmxl_atf) + tmltrd_atf_sumb(:,:) )
520
521	!-- Compute salinity residuals
522	DO jl = 1, jpltrd
523	zsmltrd2(:,:,jl) = smltrd_csum_ub(:,:,jl) + smltrd_csum_ln(:,:,jl)
524	END DO
525
526	zsmltrdm2(:,:) = 0.
527	DO jl = 1, jpltrd
528	zsmltrdm2(:,:) = zsmltrdm2(:,:) + zsmltrd2(:,:,jl)
529	END DO
530
531	zsmlres2(:,:) = zsmltot2(:,:) - &
532	( zsmltrdm2(:,:) - smltrd_sum(:,:,jpmxl_atf) + smltrd_atf_sumb(:,:) )
533
534	!-- Diagnose Asselin trend over the analysis window
535	ztmlatf2(:,:) = ztmltrd2(:,:,jpmxl_atf) - tmltrd_sum(:,:,jpmxl_atf) + tmltrd_atf_sumb(:,:)
536	zsmlatf2(:,:) = zsmltrd2(:,:,jpmxl_atf) - smltrd_sum(:,:,jpmxl_atf) + smltrd_atf_sumb(:,:)
537
538	!-- Lateral boundary conditions
539	! ... temperature ... ... salinity ...
540	CALL lbc_lnk( ztmltot2, 'T', 1. ) ; CALL lbc_lnk( zsmltot2, 'T', 1. )
541	CALL lbc_lnk( ztmlres2, 'T', 1. ) ; CALL lbc_lnk( zsmlres2, 'T', 1. )
542	DO jl = 1, jpltrd
543	CALL lbc_lnk( ztmltrd2(:,:,jl), 'T', 1. ) ! \ these will be output
544	CALL lbc_lnk( zsmltrd2(:,:,jl), 'T', 1. ) ! / in the NetCDF trends file
545	END DO
546
547	! III.3 Time evolution array swap
548	! -------------------------------
549
550	! For T/S instantaneous diagnostics
551	! ... temperature ... ... salinity ...
552	tmlbb (:,:) = tmlb (:,:) ; smlbb (:,:) = smlb (:,:)
553	tmlbn (:,:) = tml (:,:) ; smlbn (:,:) = sml (:,:)
554	tmlatfb(:,:) = tmlatfn(:,:) ; smlatfb(:,:) = smlatfn(:,:)
555
556	! For T mean diagnostics
557	tmltrd_csum_ub (:,:,:) = zfn * tmltrd_sum(:,:,:) - tmltrd_csum_ln(:,:,:)
558	tml_sumb (:,:) = tml_sum(:,:)
559	tmltrd_atf_sumb(:,:) = tmltrd_sum(:,:,jpmxl_atf)
560
561	! For S mean diagnostics
562	smltrd_csum_ub (:,:,:) = zfn * smltrd_sum(:,:,:) - smltrd_csum_ln(:,:,:)
563	sml_sumb (:,:) = sml_sum(:,:)
564	smltrd_atf_sumb(:,:) = smltrd_sum(:,:,jpmxl_atf)
565
566	! ML depth
567	hmxlbn (:,:) = hmxl (:,:)
568
569	IF( ln_ctl ) THEN
570	IF( ln_trdmxl_instant ) THEN
571	CALL prt_ctl(tab2d_1=tmlbb , clinfo1=' tmlbb - : ', mask1=tmask, ovlap=1)
572	CALL prt_ctl(tab2d_1=tmlbn , clinfo1=' tmlbn - : ', mask1=tmask, ovlap=1)
573	CALL prt_ctl(tab2d_1=tmlatfb , clinfo1=' tmlatfb - : ', mask1=tmask, ovlap=1)
574	ELSE
575	CALL prt_ctl(tab2d_1=tmlbn , clinfo1=' tmlbn - : ', mask1=tmask, ovlap=1)
576	CALL prt_ctl(tab2d_1=hmxlbn , clinfo1=' hmxlbn - : ', mask1=tmask, ovlap=1)
577	CALL prt_ctl(tab2d_1=tml_sumb , clinfo1=' tml_sumb - : ', mask1=tmask, ovlap=1)
578	CALL prt_ctl(tab2d_1=tmltrd_atf_sumb, clinfo1=' tmltrd_atf_sumb - : ', mask1=tmask, ovlap=1)
579	CALL prt_ctl(tab3d_1=tmltrd_csum_ub , clinfo1=' tmltrd_csum_ub - : ', mask1=tmask, ovlap=1, kdim=1)
580	END IF
581	END IF
582
583	! III.4 Convert to appropriate physical units
584	! -------------------------------------------
585
586	! ... temperature ... ... salinity ...
587	ztmltot (:,:) = ztmltot(:,:) * rn_ucf/zfn ; zsmltot (:,:) = zsmltot(:,:) * rn_ucf/zfn
588	ztmlres (:,:) = ztmlres(:,:) * rn_ucf/zfn ; zsmlres (:,:) = zsmlres(:,:) * rn_ucf/zfn
589	ztmlatf (:,:) = ztmlatf(:,:) * rn_ucf/zfn ; zsmlatf (:,:) = zsmlatf(:,:) * rn_ucf/zfn
590
591	tml_sum (:,:) = tml_sum (:,:) / (2zfn) ; sml_sum (:,:) = sml_sum (:,:) / (2zfn)
592	ztmltot2(:,:) = ztmltot2(:,:) * rn_ucf/zfn2 ; zsmltot2(:,:) = zsmltot2(:,:) * rn_ucf/zfn2
593	ztmltrd2(:,:,:) = ztmltrd2(:,:,:)* rn_ucf/zfn2 ; zsmltrd2(:,:,:) = zsmltrd2(:,:,:)* rn_ucf/zfn2
594	ztmlatf2(:,:) = ztmlatf2(:,:) * rn_ucf/zfn2 ; zsmlatf2(:,:) = zsmlatf2(:,:) * rn_ucf/zfn2
595	ztmlres2(:,:) = ztmlres2(:,:) * rn_ucf/zfn2 ; zsmlres2(:,:) = zsmlres2(:,:) * rn_ucf/zfn2
596
597	hmxl_sum(:,:) = hmxl_sum(:,:) / (2*zfn) ! similar to tml_sum and sml_sum
598
599	! * Debugging information *
600	IF( lldebug ) THEN
601	!
602	WRITE(numout,*)
603	WRITE(numout,*) 'trd_mxl : write trends in the Mixed Layer for debugging process:'
604	WRITE(numout,*) '~~~~~~~ '
605	WRITE(numout,*) ' TRA kt = ', kt, 'nmoymltrd = ', nmoymltrd
606	WRITE(numout,*)
607	WRITE(numout,*) ' >>>>>>>>>>>>>>>>>> TRA TEMPERATURE <<<<<<<<<<<<<<<<<<'
608	WRITE(numout,*) ' TRA ztmlres : ', SUM(ztmlres(:,:))
609	WRITE(numout,*) ' TRA ztmltot : ', SUM(ztmltot(:,:))
610	WRITE(numout,*) ' TRA tmltrdm : ', SUM(tmltrdm(:,:))
611	WRITE(numout,*) ' TRA tmlatfb : ', SUM(tmlatfb(:,:))
612	WRITE(numout,*) ' TRA tmlatfn : ', SUM(tmlatfn(:,:))
613	DO jl = 1, jpltrd
614	WRITE(numout,) ' TRA TREND INDEX jpmxl_xxx = jl = ', jl, &
615	& ' tmltrd : ', SUM(tmltrd(:,:,jl))
616	END DO
617	WRITE(numout,*) ' TRA ztmlres (jpi/2,jpj/2) : ', ztmlres (jpi/2,jpj/2)
618	WRITE(numout,*) ' TRA ztmlres2(jpi/2,jpj/2) : ', ztmlres2(jpi/2,jpj/2)
619	WRITE(numout,*)
620	WRITE(numout,*) ' >>>>>>>>>>>>>>>>>> TRA SALINITY <<<<<<<<<<<<<<<<<<'
621	WRITE(numout,*) ' TRA zsmlres : ', SUM(zsmlres(:,:))
622	WRITE(numout,*) ' TRA zsmltot : ', SUM(zsmltot(:,:))
623	WRITE(numout,*) ' TRA smltrdm : ', SUM(smltrdm(:,:))
624	WRITE(numout,*) ' TRA smlatfb : ', SUM(smlatfb(:,:))
625	WRITE(numout,*) ' TRA smlatfn : ', SUM(smlatfn(:,:))
626	DO jl = 1, jpltrd
627	WRITE(numout,) ' TRA TREND INDEX jpmxl_xxx = jl = ', jl, &
628	& ' smltrd : ', SUM(smltrd(:,:,jl))
629	END DO
630	WRITE(numout,*) ' TRA zsmlres (jpi/2,jpj/2) : ', zsmlres (jpi/2,jpj/2)
631	WRITE(numout,*) ' TRA zsmlres2(jpi/2,jpj/2) : ', zsmlres2(jpi/2,jpj/2)
632	!
633	END IF
634	!
635	END IF MODULO_NTRD
636
637	! ======================================================================
638	! IV. Write trends in the NetCDF file
639	! ======================================================================
640
641	!-- Write the trends for T/S instantaneous diagnostics
642
643	IF( ln_trdmxl_instant ) THEN
644
645	CALL iom_put( "mxl_depth", hmxl(:,:) )
646
647	!................................. ( ML temperature ) ...................................
648
649	!-- Output the fields
650	CALL iom_put( "tml" , tml (:,:) )
651	CALL iom_put( "tml_tot" , ztmltot(:,:) )
652	CALL iom_put( "tml_res" , ztmlres(:,:) )
653
654	DO jl = 1, jpltrd - 1
655	CALL iom_put( trim("tml"//ctrd(jl,2)), tmltrd (:,:,jl) )
656	END DO
657
658	CALL iom_put( trim("tml"//ctrd(jpmxl_atf,2)), ztmlatf(:,:) )
659
660	!.................................. ( ML salinity ) .....................................
661
662	!-- Output the fields
663	CALL iom_put( "sml" , sml (:,:) )
664	CALL iom_put( "sml_tot", zsmltot(:,:) )
665	CALL iom_put( "sml_res", zsmlres(:,:) )
666
667	DO jl = 1, jpltrd - 1
668	CALL iom_put( trim("sml"//ctrd(jl,2)), smltrd(:,:,jl) )
669	END DO
670
671	CALL iom_put( trim("sml"//ctrd(jpmxl_atf,2)), zsmlatf(:,:) )
672
673
674
675	ELSE !-- Write the trends for T/S mean diagnostics
676
677	CALL iom_put( "mxl_depth", hmxl_sum(:,:) )
678
679	!................................. ( ML temperature ) ...................................
680
681	!-- Output the fields
682	CALL iom_put( "tml" , tml_sum (:,:) )
683	CALL iom_put( "tml_tot" , ztmltot2(:,:) )
684	CALL iom_put( "tml_res" , ztmlres2(:,:) )
685
686	DO jl = 1, jpltrd - 1
687	CALL iom_put( trim("tml"//ctrd(jl,2)), ztmltrd2(:,:,jl) )
688	END DO
689
690	CALL iom_put( trim("tml"//ctrd(jpmxl_atf,2)), ztmlatf2(:,:) )
691
692	!.................................. ( ML salinity ) .....................................
693
694	!-- Output the fields
695	CALL iom_put( "sml" , sml_sum (:,:) )
696	CALL iom_put( "sml_tot", zsmltot2(:,:) )
697	CALL iom_put( "sml_res", zsmlres2(:,:) )
698
699	DO jl = 1, jpltrd - 1
700	CALL iom_put( trim("sml"//ctrd(jl,2)), zsmltrd2(:,:,jl) )
701	END DO
702
703	CALL iom_put( trim("sml"//ctrd(jpmxl_atf,2)), zsmlatf2(:,:) )
704	!
705	END IF
706	!
707
708	IF( MOD( itmod, nn_trd ) == 0 ) THEN
709	!
710	! III.5 Reset cumulative arrays to zero
711	! -------------------------------------
712	nmoymltrd = 0
713
714	! ... temperature ... ... salinity ...
715	tmltrdm (:,:) = 0.e0 ; smltrdm (:,:) = 0.e0
716	tmlatfm (:,:) = 0.e0 ; smlatfm (:,:) = 0.e0
717	tml_sum (:,:) = 0.e0 ; sml_sum (:,:) = 0.e0
718	tmltrd_csum_ln (:,:,:) = 0.e0 ; smltrd_csum_ln (:,:,:) = 0.e0
719	tmltrd_sum (:,:,:) = 0.e0 ; smltrd_sum (:,:,:) = 0.e0
720
721	hmxl_sum (:,:) = 0.e0
722	!
723	END IF
724
725	! ======================================================================
726	! V. Write restart file
727	! ======================================================================
728
729	IF( lrst_oce ) CALL trd_mxl_rst_write( kt )
730
731	CALL wrk_dealloc( jpi, jpj, ztmltot , zsmltot , ztmlres , zsmlres , ztmlatf , zsmlatf )
732	CALL wrk_dealloc( jpi, jpj, ztmltot2, zsmltot2, ztmlres2, zsmlres2, ztmlatf2, zsmlatf2, ztmltrdm2, zsmltrdm2 )
733	CALL wrk_dealloc( jpi, jpj, jpltrd, ztmltrd2, zsmltrd2 )
734	!
735	END SUBROUTINE trd_mxl
736
737
738	SUBROUTINE trd_mxl_init
739	!!----------------------------------------------------------------------
740	!! * ROUTINE trd_mxl_init *
741	!!
742	!! ** Purpose : computation of vertically integrated T and S budgets
743	!! from ocean surface down to control surface (NetCDF output)
744	!!----------------------------------------------------------------------
745	INTEGER :: jl ! dummy loop indices
746	INTEGER :: inum ! logical unit
747	INTEGER :: ios ! local integer
748	REAL(wp) :: zjulian, zsto, zout
749	CHARACTER (LEN=40) :: clop
750	CHARACTER (LEN=12) :: clmxl, cltu, clsu
751	!!
752	NAMELIST/namtrd_mxl/ nn_trd , cn_trdrst_in , ln_trdmxl_restart, &
753	& nn_ctls, cn_trdrst_out, ln_trdmxl_instant, rn_ucf, rn_rho_c
754	!!----------------------------------------------------------------------
755	!
756	REWIND( numnam_ref ) ! Namelist namtrd_mxl in reference namelist : mixed layer trends diagnostic
757	READ ( numnam_ref, namtrd_mxl, IOSTAT = ios, ERR = 901 )
758	901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namtrd_mxl in reference namelist', lwp )
759
760	REWIND( numnam_cfg ) ! Namelist namtrd_mxl in configuration namelist : mixed layer trends diagnostic
761	READ ( numnam_cfg, namtrd_mxl, IOSTAT = ios, ERR = 902 )
762	902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namtrd_mxl in configuration namelist', lwp )
763	IF(lwm) WRITE( numond, namtrd_mxl )
764	!
765	IF(lwp) THEN ! control print
766	WRITE(numout,*)
767	WRITE(numout,*) ' trd_mxl_init : Mixed-layer trends'
768	WRITE(numout,*) ' ~~~~~~~~~~'
769	WRITE(numout,*) ' Namelist namtrd : set trends parameters'
770	WRITE(numout,*) ' frequency of trends diagnostics (glo) nn_trd = ', nn_trd
771	WRITE(numout,*) ' density criteria used to defined the MLD rn_rho_c = ', rn_rho_c
772	WRITE(numout,*) ' control surface type (mld) nn_ctls = ', nn_ctls
773	WRITE(numout,*) ' restart for ML diagnostics ln_trdmxl_restart = ', ln_trdmxl_restart
774	WRITE(numout,*) ' instantaneous or mean ML T/S ln_trdmxl_instant = ', ln_trdmxl_instant
775	WRITE(numout,*) ' unit conversion factor rn_ucf = ', rn_ucf
776	WRITE(numout,*) ' criteria to compute the MLD rn_rho_c = ', rn_rho_c
777	ENDIF
778
779
780
781	! I.1 Check consistency of user defined preferences
782	! -------------------------------------------------
783
784	IF ( rn_rho_c /= rho_c ) CALL ctl_warn( 'Unless you have good reason to do so, you should use the value ', &
785	& 'defined in zdfmxl.F90 module to calculate the mixed layer depth' )
786
787	IF( MOD( nitend, nn_trd ) /= 0 ) THEN
788	WRITE(numout,cform_err)
789	WRITE(numout,*) ' Your nitend parameter, nitend = ', nitend
790	WRITE(numout,*) ' is no multiple of the trends diagnostics frequency '
791	WRITE(numout,*) ' you defined, nn_trd = ', nn_trd
792	WRITE(numout,*) ' This will not allow you to restart from this simulation. '
793	WRITE(numout,*) ' You should reconsider this choice. '
794	WRITE(numout,*)
795	WRITE(numout,*) ' N.B. the nitend parameter is also constrained to be a '
796	WRITE(numout,*) ' multiple of the nn_fsbc parameter '
797	nstop = nstop + 1
798	END IF
799
800	! ! allocate trdmxl arrays
801	IF( trd_mxl_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'trd_mxl_init : unable to allocate trdmxl arrays' )
802	IF( trdmxl_oce_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'trd_mxl_init : unable to allocate trdmxl_oce arrays' )
803
804
805
806	nkstp = nit000 - 1 ! current time step indicator initialization
807
808
809
810
811	! I.2 Initialize arrays to zero or read a restart file
812	! ----------------------------------------------------
813
814	nmoymltrd = 0
815
816	! ... temperature ... ... salinity ...
817	tml (:,:) = 0.e0 ; sml (:,:) = 0.e0 ! inst.
818	tmltrdm (:,:) = 0.e0 ; smltrdm (:,:) = 0.e0
819	tmlatfm (:,:) = 0.e0 ; smlatfm (:,:) = 0.e0
820	tml_sum (:,:) = 0.e0 ; sml_sum (:,:) = 0.e0 ! mean
821	tmltrd_sum (:,:,:) = 0.e0 ; smltrd_sum (:,:,:) = 0.e0
822	tmltrd_csum_ln (:,:,:) = 0.e0 ; smltrd_csum_ln (:,:,:) = 0.e0
823
824	hmxl (:,:) = 0.e0
825	hmxl_sum (:,:) = 0.e0
826
827	IF( ln_rstart .AND. ln_trdmxl_restart ) THEN
828	CALL trd_mxl_rst_read
829	ELSE
830	! ... temperature ... ... salinity ...
831	tmlb (:,:) = 0.e0 ; smlb (:,:) = 0.e0 ! inst.
832	tmlbb (:,:) = 0.e0 ; smlbb (:,:) = 0.e0
833	tmlbn (:,:) = 0.e0 ; smlbn (:,:) = 0.e0
834	tml_sumb (:,:) = 0.e0 ; sml_sumb (:,:) = 0.e0 ! mean
835	tmltrd_csum_ub (:,:,:) = 0.e0 ; smltrd_csum_ub (:,:,:) = 0.e0
836	tmltrd_atf_sumb(:,:) = 0.e0 ; smltrd_atf_sumb(:,:) = 0.e0
837	END IF
838
839	icount = 1 ; ionce = 1 ! open specifier
840
841	! I.3 Read control surface from file ctlsurf_idx
842	! ----------------------------------------------
843
844	IF( nn_ctls == 1 ) THEN
845	CALL ctl_opn( inum, 'ctlsurf_idx', 'OLD', 'UNFORMATTED', 'SEQUENTIAL', -1, numout, lwp )
846	READ ( inum, * ) nbol
847	CLOSE( inum )
848	END IF
849
850	! ======================================================================
851	! II. netCDF output initialization
852	! ======================================================================
853
854	! clmxl = legend root for netCDF output
855	IF( nn_ctls == 0 ) THEN ! control surface = mixed-layer with density criterion
856	clmxl = 'Mixed Layer ' ! (array nmln computed in zdfmxl.F90)
857	ELSE IF( nn_ctls == 1 ) THEN ! control surface = read index from file
858	clmxl = ' Bowl '
859	ELSE IF( nn_ctls >= 2 ) THEN ! control surface = model level
860	WRITE(clmxl,'(A10,I2,1X)') 'Levels 1 -', nn_ctls
861	END IF
862
863
864
865	! II.3 Define the T grid trend file (nidtrd)
866	! ------------------------------------------
867	!-- Define long and short names for the NetCDF output variables
868	! ==> choose them according to trdmxl_oce.F90 <==
869
870	ctrd(jpmxl_xad,1) = " Zonal advection" ; ctrd(jpmxl_xad,2) = "_xad"
871	ctrd(jpmxl_yad,1) = " Meridional advection" ; ctrd(jpmxl_yad,2) = "_yad"
872	ctrd(jpmxl_zad,1) = " Vertical advection" ; ctrd(jpmxl_zad,2) = "_zad"
873	ctrd(jpmxl_ldf,1) = " Lateral diffusion" ; ctrd(jpmxl_ldf,2) = "_ldf"
874	ctrd(jpmxl_for,1) = " Forcing" ; ctrd(jpmxl_for,2) = "_for"
875	ctrd(jpmxl_zdf,1) = " Vertical diff. (Kz)" ; ctrd(jpmxl_zdf,2) = "_zdf"
876	ctrd(jpmxl_bbc,1) = " Geothermal flux" ; ctrd(jpmxl_bbc,2) = "_bbc"
877	ctrd(jpmxl_bbl,1) = " Adv/diff. Bottom boundary layer" ; ctrd(jpmxl_bbl,2) = "_bbl"
878	ctrd(jpmxl_dmp,1) = " Tracer damping" ; ctrd(jpmxl_dmp,2) = "_dmp"
879	ctrd(jpmxl_npc,1) = " Non penetrative convec. adjust." ; ctrd(jpmxl_npc,2) = "_npc"
880	ctrd(jpmxl_atf,1) = " Asselin time filter" ; ctrd(jpmxl_atf,2) = "_atf"
881
882
883	!-- Define physical units
884	IF ( rn_ucf == 1. ) THEN ; cltu = "degC/s" ; clsu = "p.s.u./s"
885	ELSEIF ( rn_ucf == 3600.*24.) THEN ; cltu = "degC/day" ; clsu = "p.s.u./day"
886	ELSE ; cltu = "unknown?" ; clsu = "unknown?"
887	END IF
888	!
889	END SUBROUTINE trd_mxl_init
890
891	!!======================================================================
892	END MODULE trdmxl

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

Download in other formats: