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/UKMO/test_moci_test_suite_namelist_read/NEMOGCM/NEMO/OPA_SRC/TRD – NEMO

Context Navigation

source: branches/UKMO/test_moci_test_suite_namelist_read/NEMOGCM/NEMO/OPA_SRC/TRD/trdmxl.F90 @ 9366

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

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

Download in other formats: