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

Context Navigation

source: branches/2015/dev_CMCC_merge_2015/NEMOGCM/NEMO/OPA_SRC/TRD/trdmxl.F90 @ 6051

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

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

Download in other formats: