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/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/TRD – NEMO

Context Navigation

source: branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/TRD/trdmxl.F90 @ 9019

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

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

Download in other formats: