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

Context Navigation

source: branches/2013/dev_r3858_CNRS3_Ediag/NEMOGCM/NEMO/OPA_SRC/TRD/trdmxl.F90 @ 3878

Last change on this file since 3878 was 3876, checked in by gm, 11 years ago
dev_r3858_CNRS3_Ediag: #927 phasing with 2011/dev_r3309_LOCEAN12_Ediag branche + mxl diag update
File size: 44.0 KB

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