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

Context Navigation

source: branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/OPA_SRC/TRD/trdmxl.F90 @ 5770

Last change on this file since 5770 was 5770, checked in by gm, 9 years ago
#1593: LDF-ADV, step II.2: phasing the improvements/simplifications of advective tracer trend (see wiki page)
Property svn:keywords set to `Id`
File size: 44.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 ! lateral diffusion: eddy diffusivity & EIV coeff.
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
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 (or dimg) 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	#if defined key_dimgout
351	INTEGER :: iyear,imon,iday
352	CHARACTER(LEN=80) :: cltext, clmode
353	#endif
354	!!----------------------------------------------------------------------
355
356	CALL wrk_alloc( jpi, jpj, ztmltot , zsmltot , ztmlres , zsmlres , ztmlatf , zsmlatf )
357	CALL wrk_alloc( jpi, jpj, ztmltot2, zsmltot2, ztmlres2, zsmlres2, ztmlatf2, zsmlatf2, ztmltrdm2, zsmltrdm2 )
358	CALL wrk_alloc( jpi, jpj, jpltrd, ztmltrd2, zsmltrd2 )
359
360	! ======================================================================
361	! II. Cumulate the trends over the analysis window
362	! ======================================================================
363
364	ztmltrd2(:,:,:) = 0.e0 ; zsmltrd2(:,:,:) = 0.e0 ! <<< reset arrays to zero
365	ztmltot2(:,:) = 0.e0 ; zsmltot2(:,:) = 0.e0
366	ztmlres2(:,:) = 0.e0 ; zsmlres2(:,:) = 0.e0
367	ztmlatf2(:,:) = 0.e0 ; zsmlatf2(:,:) = 0.e0
368
369	! II.1 Set before values of vertically average T and S
370	! ----------------------------------------------------
371	IF( kt > nit000 ) THEN
372	! ... temperature ... ... salinity ...
373	tmlb (:,:) = tml (:,:) ; smlb (:,:) = sml (:,:)
374	tmlatfn(:,:) = tmltrd(:,:,jpmxl_atf) ; smlatfn(:,:) = smltrd(:,:,jpmxl_atf)
375	END IF
376
377
378	! II.3 Initialize mixed-layer "before" arrays for the 1rst analysis window
379	! ------------------------------------------------------------------------
380	IF( kt == 2 ) THEN ! i.e. ( .NOT. ln_rstart ).AND.( kt == nit000 + 1)
381	!
382	! ... temperature ... ... salinity ...
383	tmlbb (:,:) = tmlb (:,:) ; smlbb (:,:) = smlb (:,:)
384	tmlbn (:,:) = tml (:,:) ; smlbn (:,:) = sml (:,:)
385	tmlatfb(:,:) = tmlatfn(:,:) ; smlatfb(:,:) = smlatfn(:,:)
386
387	tmltrd_csum_ub (:,:,:) = 0.e0 ; smltrd_csum_ub (:,:,:) = 0.e0
388	tmltrd_atf_sumb(:,:) = 0.e0 ; smltrd_atf_sumb(:,:) = 0.e0
389
390	hmxlbn(:,:) = hmxl(:,:)
391
392	IF( ln_ctl ) THEN
393	WRITE(numout,*) ' we reach kt == nit000 + 1 = ', nit000+1
394	CALL prt_ctl(tab2d_1=tmlbb , clinfo1=' tmlbb - : ', mask1=tmask, ovlap=1)
395	CALL prt_ctl(tab2d_1=tmlbn , clinfo1=' tmlbn - : ', mask1=tmask, ovlap=1)
396	CALL prt_ctl(tab2d_1=tmlatfb , clinfo1=' tmlatfb - : ', mask1=tmask, ovlap=1)
397	END IF
398	!
399	END IF
400
401	IF( ( ln_rstart ) .AND. ( kt == nit000 ) .AND. ( ln_ctl ) ) THEN
402	IF( ln_trdmxl_instant ) THEN
403	WRITE(numout,*) ' restart from kt == nit000 = ', nit000
404	CALL prt_ctl(tab2d_1=tmlbb , clinfo1=' tmlbb - : ', mask1=tmask, ovlap=1)
405	CALL prt_ctl(tab2d_1=tmlbn , clinfo1=' tmlbn - : ', mask1=tmask, ovlap=1)
406	CALL prt_ctl(tab2d_1=tmlatfb , clinfo1=' tmlatfb - : ', mask1=tmask, ovlap=1)
407	ELSE
408	WRITE(numout,*) ' restart from kt == nit000 = ', nit000
409	CALL prt_ctl(tab2d_1=tmlbn , clinfo1=' tmlbn - : ', mask1=tmask, ovlap=1)
410	CALL prt_ctl(tab2d_1=hmxlbn , clinfo1=' hmxlbn - : ', mask1=tmask, ovlap=1)
411	CALL prt_ctl(tab2d_1=tml_sumb , clinfo1=' tml_sumb - : ', mask1=tmask, ovlap=1)
412	CALL prt_ctl(tab2d_1=tmltrd_atf_sumb, clinfo1=' tmltrd_atf_sumb - : ', mask1=tmask, ovlap=1)
413	CALL prt_ctl(tab3d_1=tmltrd_csum_ub , clinfo1=' tmltrd_csum_ub - : ', mask1=tmask, ovlap=1, kdim=1)
414	END IF
415	END IF
416
417	! II.4 Cumulated trends over the analysis period
418	! ----------------------------------------------
419	!
420	! [ 1rst analysis window ] [ 2nd analysis window ]
421	!
422	! o---[--o-----o-----o-----o--]-[--o-----o-----o-----o-----o--]---o-----o--> time steps
423	! nn_trd 2*nn_trd etc.
424	! 1 2 3 4 =5 e.g. =10
425	!
426	IF( ( kt >= 2 ).OR.( ln_rstart ) ) THEN
427	!
428	nmoymltrd = nmoymltrd + 1
429
430	! ... Cumulate over BOTH physical contributions AND over time steps
431	DO jl = 1, jpltrd
432	tmltrdm(:,:) = tmltrdm(:,:) + tmltrd(:,:,jl)
433	smltrdm(:,:) = smltrdm(:,:) + smltrd(:,:,jl)
434	END DO
435
436	! ... Special handling of the Asselin trend
437	tmlatfm(:,:) = tmlatfm(:,:) + tmlatfn(:,:)
438	smlatfm(:,:) = smlatfm(:,:) + smlatfn(:,:)
439
440	! ... Trends associated with the time mean of the ML T/S
441	tmltrd_sum (:,:,:) = tmltrd_sum (:,:,:) + tmltrd (:,:,:) ! tem
442	tmltrd_csum_ln(:,:,:) = tmltrd_csum_ln(:,:,:) + tmltrd_sum(:,:,:)
443	tml_sum (:,:) = tml_sum (:,:) + tml (:,:)
444	smltrd_sum (:,:,:) = smltrd_sum (:,:,:) + smltrd (:,:,:) ! sal
445	smltrd_csum_ln(:,:,:) = smltrd_csum_ln(:,:,:) + smltrd_sum(:,:,:)
446	sml_sum (:,:) = sml_sum (:,:) + sml (:,:)
447	hmxl_sum (:,:) = hmxl_sum (:,:) + hmxl (:,:) ! rmxl
448	!
449	END IF
450
451	! ======================================================================
452	! III. Prepare fields for output (get here ONCE PER ANALYSIS PERIOD)
453	! ======================================================================
454
455	! Convert to appropriate physical units
456	! N.B. It may be useful to check IOIPSL time averaging with :
457	! tmltrd (:,:,:) = 1. ; smltrd (:,:,:) = 1.
458	tmltrd(:,:,:) = tmltrd(:,:,:) * rn_ucf ! (actually needed for 1:jpltrd-1, but trdmxl(:,:,jpltrd)
459	smltrd(:,:,:) = smltrd(:,:,:) * rn_ucf ! is no longer used, and is reset to 0. at next time step)
460
461	! define time axis
462	it = kt
463	itmod = kt - nit000 + 1
464
465	MODULO_NTRD : IF( MOD( itmod, nn_trd ) == 0 ) THEN ! nitend MUST be multiple of nn_trd
466	!
467	ztmltot (:,:) = 0.e0 ; zsmltot (:,:) = 0.e0 ! reset arrays to zero
468	ztmlres (:,:) = 0.e0 ; zsmlres (:,:) = 0.e0
469	ztmltot2(:,:) = 0.e0 ; zsmltot2(:,:) = 0.e0
470	ztmlres2(:,:) = 0.e0 ; zsmlres2(:,:) = 0.e0
471
472	zfn = REAL( nmoymltrd, wp ) ; zfn2 = zfn * zfn
473
474	! III.1 Prepare fields for output ("instantaneous" diagnostics)
475	! -------------------------------------------------------------
476
477	!-- Compute total trends
478	ztmltot(:,:) = ( tml(:,:) - tmlbn(:,:) + tmlb(:,:) - tmlbb(:,:) ) / p2dt
479	zsmltot(:,:) = ( sml(:,:) - smlbn(:,:) + smlb(:,:) - smlbb(:,:) ) / p2dt
480
481	!-- Compute residuals
482	ztmlres(:,:) = ztmltot(:,:) - ( tmltrdm(:,:) - tmlatfn(:,:) + tmlatfb(:,:) )
483	zsmlres(:,:) = zsmltot(:,:) - ( smltrdm(:,:) - smlatfn(:,:) + smlatfb(:,:) )
484
485	!-- Diagnose Asselin trend over the analysis window
486	ztmlatf(:,:) = tmlatfm(:,:) - tmlatfn(:,:) + tmlatfb(:,:)
487	zsmlatf(:,:) = smlatfm(:,:) - smlatfn(:,:) + smlatfb(:,:)
488
489	!-- Lateral boundary conditions
490	! ... temperature ... ... salinity ...
491	CALL lbc_lnk( ztmltot , 'T', 1. ) ; CALL lbc_lnk( zsmltot , 'T', 1. )
492	CALL lbc_lnk( ztmlres , 'T', 1. ) ; CALL lbc_lnk( zsmlres , 'T', 1. )
493	CALL lbc_lnk( ztmlatf , 'T', 1. ) ; CALL lbc_lnk( zsmlatf , 'T', 1. )
494
495
496	! III.2 Prepare fields for output ("mean" diagnostics)
497	! ----------------------------------------------------
498
499	!-- Update the ML depth time sum (to build the Leap-Frog time mean)
500	hmxl_sum(:,:) = hmxlbn(:,:) + 2 * ( hmxl_sum(:,:) - hmxl(:,:) ) + hmxl(:,:)
501
502	!-- Compute temperature total trends
503	tml_sum (:,:) = tmlbn(:,:) + 2 * ( tml_sum(:,:) - tml(:,:) ) + tml(:,:)
504	ztmltot2(:,:) = ( tml_sum(:,:) - tml_sumb(:,:) ) / p2dt ! now in degC/s
505
506	!-- Compute salinity total trends
507	sml_sum (:,:) = smlbn(:,:) + 2 * ( sml_sum(:,:) - sml(:,:) ) + sml(:,:)
508	zsmltot2(:,:) = ( sml_sum(:,:) - sml_sumb(:,:) ) / p2dt ! now in psu/s
509
510	!-- Compute temperature residuals
511	DO jl = 1, jpltrd
512	ztmltrd2(:,:,jl) = tmltrd_csum_ub(:,:,jl) + tmltrd_csum_ln(:,:,jl)
513	END DO
514
515	ztmltrdm2(:,:) = 0.e0
516	DO jl = 1, jpltrd
517	ztmltrdm2(:,:) = ztmltrdm2(:,:) + ztmltrd2(:,:,jl)
518	END DO
519
520	ztmlres2(:,:) = ztmltot2(:,:) - &
521	( ztmltrdm2(:,:) - tmltrd_sum(:,:,jpmxl_atf) + tmltrd_atf_sumb(:,:) )
522
523	!-- Compute salinity residuals
524	DO jl = 1, jpltrd
525	zsmltrd2(:,:,jl) = smltrd_csum_ub(:,:,jl) + smltrd_csum_ln(:,:,jl)
526	END DO
527
528	zsmltrdm2(:,:) = 0.
529	DO jl = 1, jpltrd
530	zsmltrdm2(:,:) = zsmltrdm2(:,:) + zsmltrd2(:,:,jl)
531	END DO
532
533	zsmlres2(:,:) = zsmltot2(:,:) - &
534	( zsmltrdm2(:,:) - smltrd_sum(:,:,jpmxl_atf) + smltrd_atf_sumb(:,:) )
535
536	!-- Diagnose Asselin trend over the analysis window
537	ztmlatf2(:,:) = ztmltrd2(:,:,jpmxl_atf) - tmltrd_sum(:,:,jpmxl_atf) + tmltrd_atf_sumb(:,:)
538	zsmlatf2(:,:) = zsmltrd2(:,:,jpmxl_atf) - smltrd_sum(:,:,jpmxl_atf) + smltrd_atf_sumb(:,:)
539
540	!-- Lateral boundary conditions
541	! ... temperature ... ... salinity ...
542	CALL lbc_lnk( ztmltot2, 'T', 1. ) ; CALL lbc_lnk( zsmltot2, 'T', 1. )
543	CALL lbc_lnk( ztmlres2, 'T', 1. ) ; CALL lbc_lnk( zsmlres2, 'T', 1. )
544	DO jl = 1, jpltrd
545	CALL lbc_lnk( ztmltrd2(:,:,jl), 'T', 1. ) ! \ these will be output
546	CALL lbc_lnk( zsmltrd2(:,:,jl), 'T', 1. ) ! / in the NetCDF trends file
547	END DO
548
549	! III.3 Time evolution array swap
550	! -------------------------------
551
552	! For T/S instantaneous diagnostics
553	! ... temperature ... ... salinity ...
554	tmlbb (:,:) = tmlb (:,:) ; smlbb (:,:) = smlb (:,:)
555	tmlbn (:,:) = tml (:,:) ; smlbn (:,:) = sml (:,:)
556	tmlatfb(:,:) = tmlatfn(:,:) ; smlatfb(:,:) = smlatfn(:,:)
557
558	! For T mean diagnostics
559	tmltrd_csum_ub (:,:,:) = zfn * tmltrd_sum(:,:,:) - tmltrd_csum_ln(:,:,:)
560	tml_sumb (:,:) = tml_sum(:,:)
561	tmltrd_atf_sumb(:,:) = tmltrd_sum(:,:,jpmxl_atf)
562
563	! For S mean diagnostics
564	smltrd_csum_ub (:,:,:) = zfn * smltrd_sum(:,:,:) - smltrd_csum_ln(:,:,:)
565	sml_sumb (:,:) = sml_sum(:,:)
566	smltrd_atf_sumb(:,:) = smltrd_sum(:,:,jpmxl_atf)
567
568	! ML depth
569	hmxlbn (:,:) = hmxl (:,:)
570
571	IF( ln_ctl ) THEN
572	IF( ln_trdmxl_instant ) THEN
573	CALL prt_ctl(tab2d_1=tmlbb , clinfo1=' tmlbb - : ', mask1=tmask, ovlap=1)
574	CALL prt_ctl(tab2d_1=tmlbn , clinfo1=' tmlbn - : ', mask1=tmask, ovlap=1)
575	CALL prt_ctl(tab2d_1=tmlatfb , clinfo1=' tmlatfb - : ', mask1=tmask, ovlap=1)
576	ELSE
577	CALL prt_ctl(tab2d_1=tmlbn , clinfo1=' tmlbn - : ', mask1=tmask, ovlap=1)
578	CALL prt_ctl(tab2d_1=hmxlbn , clinfo1=' hmxlbn - : ', mask1=tmask, ovlap=1)
579	CALL prt_ctl(tab2d_1=tml_sumb , clinfo1=' tml_sumb - : ', mask1=tmask, ovlap=1)
580	CALL prt_ctl(tab2d_1=tmltrd_atf_sumb, clinfo1=' tmltrd_atf_sumb - : ', mask1=tmask, ovlap=1)
581	CALL prt_ctl(tab3d_1=tmltrd_csum_ub , clinfo1=' tmltrd_csum_ub - : ', mask1=tmask, ovlap=1, kdim=1)
582	END IF
583	END IF
584
585	! III.4 Convert to appropriate physical units
586	! -------------------------------------------
587
588	! ... temperature ... ... salinity ...
589	ztmltot (:,:) = ztmltot(:,:) * rn_ucf/zfn ; zsmltot (:,:) = zsmltot(:,:) * rn_ucf/zfn
590	ztmlres (:,:) = ztmlres(:,:) * rn_ucf/zfn ; zsmlres (:,:) = zsmlres(:,:) * rn_ucf/zfn
591	ztmlatf (:,:) = ztmlatf(:,:) * rn_ucf/zfn ; zsmlatf (:,:) = zsmlatf(:,:) * rn_ucf/zfn
592
593	tml_sum (:,:) = tml_sum (:,:) / (2zfn) ; sml_sum (:,:) = sml_sum (:,:) / (2zfn)
594	ztmltot2(:,:) = ztmltot2(:,:) * rn_ucf/zfn2 ; zsmltot2(:,:) = zsmltot2(:,:) * rn_ucf/zfn2
595	ztmltrd2(:,:,:) = ztmltrd2(:,:,:)* rn_ucf/zfn2 ; zsmltrd2(:,:,:) = zsmltrd2(:,:,:)* rn_ucf/zfn2
596	ztmlatf2(:,:) = ztmlatf2(:,:) * rn_ucf/zfn2 ; zsmlatf2(:,:) = zsmlatf2(:,:) * rn_ucf/zfn2
597	ztmlres2(:,:) = ztmlres2(:,:) * rn_ucf/zfn2 ; zsmlres2(:,:) = zsmlres2(:,:) * rn_ucf/zfn2
598
599	hmxl_sum(:,:) = hmxl_sum(:,:) / (2*zfn) ! similar to tml_sum and sml_sum
600
601	! * Debugging information *
602	IF( lldebug ) THEN
603	!
604	WRITE(numout,*)
605	WRITE(numout,*) 'trd_mxl : write trends in the Mixed Layer for debugging process:'
606	WRITE(numout,*) '~~~~~~~ '
607	WRITE(numout,*) ' TRA kt = ', kt, 'nmoymltrd = ', nmoymltrd
608	WRITE(numout,*)
609	WRITE(numout,*) ' >>>>>>>>>>>>>>>>>> TRA TEMPERATURE <<<<<<<<<<<<<<<<<<'
610	WRITE(numout,*) ' TRA ztmlres : ', SUM(ztmlres(:,:))
611	WRITE(numout,*) ' TRA ztmltot : ', SUM(ztmltot(:,:))
612	WRITE(numout,*) ' TRA tmltrdm : ', SUM(tmltrdm(:,:))
613	WRITE(numout,*) ' TRA tmlatfb : ', SUM(tmlatfb(:,:))
614	WRITE(numout,*) ' TRA tmlatfn : ', SUM(tmlatfn(:,:))
615	DO jl = 1, jpltrd
616	WRITE(numout,) ' TRA TREND INDEX jpmxl_xxx = jl = ', jl, &
617	& ' tmltrd : ', SUM(tmltrd(:,:,jl))
618	END DO
619	WRITE(numout,*) ' TRA ztmlres (jpi/2,jpj/2) : ', ztmlres (jpi/2,jpj/2)
620	WRITE(numout,*) ' TRA ztmlres2(jpi/2,jpj/2) : ', ztmlres2(jpi/2,jpj/2)
621	WRITE(numout,*)
622	WRITE(numout,*) ' >>>>>>>>>>>>>>>>>> TRA SALINITY <<<<<<<<<<<<<<<<<<'
623	WRITE(numout,*) ' TRA zsmlres : ', SUM(zsmlres(:,:))
624	WRITE(numout,*) ' TRA zsmltot : ', SUM(zsmltot(:,:))
625	WRITE(numout,*) ' TRA smltrdm : ', SUM(smltrdm(:,:))
626	WRITE(numout,*) ' TRA smlatfb : ', SUM(smlatfb(:,:))
627	WRITE(numout,*) ' TRA smlatfn : ', SUM(smlatfn(:,:))
628	DO jl = 1, jpltrd
629	WRITE(numout,) ' TRA TREND INDEX jpmxl_xxx = jl = ', jl, &
630	& ' smltrd : ', SUM(smltrd(:,:,jl))
631	END DO
632	WRITE(numout,*) ' TRA zsmlres (jpi/2,jpj/2) : ', zsmlres (jpi/2,jpj/2)
633	WRITE(numout,*) ' TRA zsmlres2(jpi/2,jpj/2) : ', zsmlres2(jpi/2,jpj/2)
634	!
635	END IF
636	!
637	END IF MODULO_NTRD
638
639	! ======================================================================
640	! IV. Write trends in the NetCDF file
641	! ======================================================================
642
643	!-- Write the trends for T/S instantaneous diagnostics
644
645	IF( ln_trdmxl_instant ) THEN
646
647	CALL iom_put( "mxl_depth", hmxl(:,:) )
648
649	!................................. ( ML temperature ) ...................................
650
651	!-- Output the fields
652	CALL iom_put( "tml" , tml (:,:) )
653	CALL iom_put( "tml_tot" , ztmltot(:,:) )
654	CALL iom_put( "tml_res" , ztmlres(:,:) )
655
656	DO jl = 1, jpltrd - 1
657	CALL iom_put( trim("tml"//ctrd(jl,2)), tmltrd (:,:,jl) )
658	END DO
659
660	CALL iom_put( trim("tml"//ctrd(jpmxl_atf,2)), ztmlatf(:,:) )
661
662	!.................................. ( ML salinity ) .....................................
663
664	!-- Output the fields
665	CALL iom_put( "sml" , sml (:,:) )
666	CALL iom_put( "sml_tot", zsmltot(:,:) )
667	CALL iom_put( "sml_res", zsmlres(:,:) )
668
669	DO jl = 1, jpltrd - 1
670	CALL iom_put( trim("sml"//ctrd(jl,2)), smltrd(:,:,jl) )
671	END DO
672
673	CALL iom_put( trim("sml"//ctrd(jpmxl_atf,2)), zsmlatf(:,:) )
674
675
676
677	ELSE !-- Write the trends for T/S mean diagnostics
678
679	CALL iom_put( "mxl_depth", hmxl_sum(:,:) )
680
681	!................................. ( ML temperature ) ...................................
682
683	!-- Output the fields
684	CALL iom_put( "tml" , tml_sum (:,:) )
685	CALL iom_put( "tml_tot" , ztmltot2(:,:) )
686	CALL iom_put( "tml_res" , ztmlres2(:,:) )
687
688	DO jl = 1, jpltrd - 1
689	CALL iom_put( trim("tml"//ctrd(jl,2)), ztmltrd2(:,:,jl) )
690	END DO
691
692	CALL iom_put( trim("tml"//ctrd(jpmxl_atf,2)), ztmlatf2(:,:) )
693
694	!.................................. ( ML salinity ) .....................................
695
696	!-- Output the fields
697	CALL iom_put( "sml" , sml_sum (:,:) )
698	CALL iom_put( "sml_tot", zsmltot2(:,:) )
699	CALL iom_put( "sml_res", zsmlres2(:,:) )
700
701	DO jl = 1, jpltrd - 1
702	CALL iom_put( trim("sml"//ctrd(jl,2)), zsmltrd2(:,:,jl) )
703	END DO
704
705	CALL iom_put( trim("sml"//ctrd(jpmxl_atf,2)), zsmlatf2(:,:) )
706	!
707	END IF
708	!
709
710	IF( MOD( itmod, nn_trd ) == 0 ) THEN
711	!
712	! III.5 Reset cumulative arrays to zero
713	! -------------------------------------
714	nmoymltrd = 0
715
716	! ... temperature ... ... salinity ...
717	tmltrdm (:,:) = 0.e0 ; smltrdm (:,:) = 0.e0
718	tmlatfm (:,:) = 0.e0 ; smlatfm (:,:) = 0.e0
719	tml_sum (:,:) = 0.e0 ; sml_sum (:,:) = 0.e0
720	tmltrd_csum_ln (:,:,:) = 0.e0 ; smltrd_csum_ln (:,:,:) = 0.e0
721	tmltrd_sum (:,:,:) = 0.e0 ; smltrd_sum (:,:,:) = 0.e0
722
723	hmxl_sum (:,:) = 0.e0
724	!
725	END IF
726
727	! ======================================================================
728	! V. Write restart file
729	! ======================================================================
730
731	IF( lrst_oce ) CALL trd_mxl_rst_write( kt )
732
733	CALL wrk_dealloc( jpi, jpj, ztmltot , zsmltot , ztmlres , zsmlres , ztmlatf , zsmlatf )
734	CALL wrk_dealloc( jpi, jpj, ztmltot2, zsmltot2, ztmlres2, zsmlres2, ztmlatf2, zsmlatf2, ztmltrdm2, zsmltrdm2 )
735	CALL wrk_dealloc( jpi, jpj, jpltrd, ztmltrd2, zsmltrd2 )
736	!
737	END SUBROUTINE trd_mxl
738
739
740	SUBROUTINE trd_mxl_init
741	!!----------------------------------------------------------------------
742	!! * ROUTINE trd_mxl_init *
743	!!
744	!! ** Purpose : computation of vertically integrated T and S budgets
745	!! from ocean surface down to control surface (NetCDF output)
746	!!----------------------------------------------------------------------
747	INTEGER :: jl ! dummy loop indices
748	INTEGER :: inum ! logical unit
749	INTEGER :: ios ! local integer
750	REAL(wp) :: zjulian, zsto, zout
751	CHARACTER (LEN=40) :: clop
752	CHARACTER (LEN=12) :: clmxl, cltu, clsu
753	!!
754	NAMELIST/namtrd_mxl/ nn_trd , cn_trdrst_in , ln_trdmxl_restart, &
755	& nn_ctls, cn_trdrst_out, ln_trdmxl_instant, rn_ucf, rn_rho_c
756	!!----------------------------------------------------------------------
757	!
758	REWIND( numnam_ref ) ! Namelist namtrd_mxl in reference namelist : mixed layer trends diagnostic
759	READ ( numnam_ref, namtrd_mxl, IOSTAT = ios, ERR = 901 )
760	901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namtrd_mxl in reference namelist', lwp )
761
762	REWIND( numnam_cfg ) ! Namelist namtrd_mxl in configuration namelist : mixed layer trends diagnostic
763	READ ( numnam_cfg, namtrd_mxl, IOSTAT = ios, ERR = 902 )
764	902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namtrd_mxl in configuration namelist', lwp )
765	IF(lwm) WRITE( numond, namtrd_mxl )
766	!
767	IF(lwp) THEN ! control print
768	WRITE(numout,*)
769	WRITE(numout,*) ' trd_mxl_init : Mixed-layer trends'
770	WRITE(numout,*) ' ~~~~~~~~~~'
771	WRITE(numout,*) ' Namelist namtrd : set trends parameters'
772	WRITE(numout,*) ' frequency of trends diagnostics (glo) nn_trd = ', nn_trd
773	WRITE(numout,*) ' density criteria used to defined the MLD rn_rho_c = ', rn_rho_c
774	WRITE(numout,*) ' control surface type (mld) nn_ctls = ', nn_ctls
775	WRITE(numout,*) ' restart for ML diagnostics ln_trdmxl_restart = ', ln_trdmxl_restart
776	WRITE(numout,*) ' instantaneous or mean ML T/S ln_trdmxl_instant = ', ln_trdmxl_instant
777	WRITE(numout,*) ' unit conversion factor rn_ucf = ', rn_ucf
778	WRITE(numout,*) ' criteria to compute the MLD rn_rho_c = ', rn_rho_c
779	ENDIF
780
781
782
783	! I.1 Check consistency of user defined preferences
784	! -------------------------------------------------
785
786	IF ( rn_rho_c /= rho_c ) CALL ctl_warn( 'Unless you have good reason to do so, you should use the value ', &
787	& 'defined in zdfmxl.F90 module to calculate the mixed layer depth' )
788
789	IF( MOD( nitend, nn_trd ) /= 0 ) THEN
790	WRITE(numout,cform_err)
791	WRITE(numout,*) ' Your nitend parameter, nitend = ', nitend
792	WRITE(numout,*) ' is no multiple of the trends diagnostics frequency '
793	WRITE(numout,*) ' you defined, nn_trd = ', nn_trd
794	WRITE(numout,*) ' This will not allow you to restart from this simulation. '
795	WRITE(numout,*) ' You should reconsider this choice. '
796	WRITE(numout,*)
797	WRITE(numout,*) ' N.B. the nitend parameter is also constrained to be a '
798	WRITE(numout,*) ' multiple of the nn_fsbc parameter '
799	nstop = nstop + 1
800	END IF
801
802	! ! allocate trdmxl arrays
803	IF( trd_mxl_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'trd_mxl_init : unable to allocate trdmxl arrays' )
804	IF( trdmxl_oce_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'trd_mxl_init : unable to allocate trdmxl_oce arrays' )
805
806
807
808	nkstp = nit000 - 1 ! current time step indicator initialization
809
810
811
812
813	! I.2 Initialize arrays to zero or read a restart file
814	! ----------------------------------------------------
815
816	nmoymltrd = 0
817
818	! ... temperature ... ... salinity ...
819	tml (:,:) = 0.e0 ; sml (:,:) = 0.e0 ! inst.
820	tmltrdm (:,:) = 0.e0 ; smltrdm (:,:) = 0.e0
821	tmlatfm (:,:) = 0.e0 ; smlatfm (:,:) = 0.e0
822	tml_sum (:,:) = 0.e0 ; sml_sum (:,:) = 0.e0 ! mean
823	tmltrd_sum (:,:,:) = 0.e0 ; smltrd_sum (:,:,:) = 0.e0
824	tmltrd_csum_ln (:,:,:) = 0.e0 ; smltrd_csum_ln (:,:,:) = 0.e0
825
826	hmxl (:,:) = 0.e0
827	hmxl_sum (:,:) = 0.e0
828
829	IF( ln_rstart .AND. ln_trdmxl_restart ) THEN
830	CALL trd_mxl_rst_read
831	ELSE
832	! ... temperature ... ... salinity ...
833	tmlb (:,:) = 0.e0 ; smlb (:,:) = 0.e0 ! inst.
834	tmlbb (:,:) = 0.e0 ; smlbb (:,:) = 0.e0
835	tmlbn (:,:) = 0.e0 ; smlbn (:,:) = 0.e0
836	tml_sumb (:,:) = 0.e0 ; sml_sumb (:,:) = 0.e0 ! mean
837	tmltrd_csum_ub (:,:,:) = 0.e0 ; smltrd_csum_ub (:,:,:) = 0.e0
838	tmltrd_atf_sumb(:,:) = 0.e0 ; smltrd_atf_sumb(:,:) = 0.e0
839	END IF
840
841	icount = 1 ; ionce = 1 ! open specifier
842
843	! I.3 Read control surface from file ctlsurf_idx
844	! ----------------------------------------------
845
846	IF( nn_ctls == 1 ) THEN
847	CALL ctl_opn( inum, 'ctlsurf_idx', 'OLD', 'UNFORMATTED', 'SEQUENTIAL', -1, numout, lwp )
848	READ ( inum, * ) nbol
849	CLOSE( inum )
850	END IF
851
852	! ======================================================================
853	! II. netCDF output initialization
854	! ======================================================================
855
856	! clmxl = legend root for netCDF output
857	IF( nn_ctls == 0 ) THEN ! control surface = mixed-layer with density criterion
858	clmxl = 'Mixed Layer ' ! (array nmln computed in zdfmxl.F90)
859	ELSE IF( nn_ctls == 1 ) THEN ! control surface = read index from file
860	clmxl = ' Bowl '
861	ELSE IF( nn_ctls >= 2 ) THEN ! control surface = model level
862	WRITE(clmxl,'(A10,I2,1X)') 'Levels 1 -', nn_ctls
863	END IF
864
865
866
867	! II.3 Define the T grid trend file (nidtrd)
868	! ------------------------------------------
869	!-- Define long and short names for the NetCDF output variables
870	! ==> choose them according to trdmxl_oce.F90 <==
871
872	ctrd(jpmxl_xad,1) = " Zonal advection" ; ctrd(jpmxl_xad,2) = "_xad"
873	ctrd(jpmxl_yad,1) = " Meridional advection" ; ctrd(jpmxl_yad,2) = "_yad"
874	ctrd(jpmxl_zad,1) = " Vertical advection" ; ctrd(jpmxl_zad,2) = "_zad"
875	ctrd(jpmxl_ldf,1) = " Lateral diffusion" ; ctrd(jpmxl_ldf,2) = "_ldf"
876	ctrd(jpmxl_for,1) = " Forcing" ; ctrd(jpmxl_for,2) = "_for"
877	ctrd(jpmxl_zdf,1) = " Vertical diff. (Kz)" ; ctrd(jpmxl_zdf,2) = "_zdf"
878	ctrd(jpmxl_bbc,1) = " Geothermal flux" ; ctrd(jpmxl_bbc,2) = "_bbc"
879	ctrd(jpmxl_bbl,1) = " Adv/diff. Bottom boundary layer" ; ctrd(jpmxl_bbl,2) = "_bbl"
880	ctrd(jpmxl_dmp,1) = " Tracer damping" ; ctrd(jpmxl_dmp,2) = "_dmp"
881	ctrd(jpmxl_npc,1) = " Non penetrative convec. adjust." ; ctrd(jpmxl_npc,2) = "_npc"
882	ctrd(jpmxl_atf,1) = " Asselin time filter" ; ctrd(jpmxl_atf,2) = "_atf"
883
884
885	!-- Define physical units
886	IF ( rn_ucf == 1. ) THEN ; cltu = "degC/s" ; clsu = "p.s.u./s"
887	ELSEIF ( rn_ucf == 3600.*24.) THEN ; cltu = "degC/day" ; clsu = "p.s.u./day"
888	ELSE ; cltu = "unknown?" ; clsu = "unknown?"
889	END IF
890	!
891	END SUBROUTINE trd_mxl_init
892
893	!!======================================================================
894	END MODULE trdmxl

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

Download in other formats: