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.

trdmld.F90 in branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/TRD – NEMO

Context Navigation

source: branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/TRD/trdmld.F90 @ 4409

Last change on this file since 4409 was 4409, checked in by trackstand2, 10 years ago
Changes to allow jpk to be modified to deepest level within a subdomain. jpkorig holds original value.
Property svn:keywords set to `Id`
File size: 52.1 KB

Line
1	MODULE trdmld
2	!!======================================================================
3	!! * MODULE trdmld *
4	!! Ocean diagnostics: mixed layer T-S trends
5	!!=====================================================================
6	!! History : ! 95-04 (J. Vialard) Original code
7	!! ! 97-02 (E. Guilyardi) Adaptation global + base cmo
8	!! ! 99-09 (E. Guilyardi) Re-writing + netCDF output
9	!! 8.5 ! 02-06 (G. Madec) F90: Free form and module
10	!! 9.0 ! 04-08 (C. Talandier) New trends organization
11	!! ! 05-05 (C. Deltel) Diagnose trends of time averaged ML T & S
12	!!----------------------------------------------------------------------
13	#if defined key_trdmld \|\| defined key_esopa
14	!!----------------------------------------------------------------------
15	!! 'key_trdmld' mixed layer trend diagnostics
16	!!----------------------------------------------------------------------
17	!! trd_mld : T and S cumulated trends averaged over the mixed layer
18	!! trd_mld_zint : T and S trends vertical integration
19	!! trd_mld_init : initialization step
20	!!----------------------------------------------------------------------
21	USE oce ! ocean dynamics and tracers variables
22	USE dom_oce ! ocean space and time domain variables
23	USE trdmod_oce ! ocean variables trends
24	USE trdmld_oce ! ocean variables trends
25	USE ldftra_oce ! ocean active tracers lateral physics
26	USE zdf_oce ! ocean vertical physics
27	USE in_out_manager ! I/O manager
28	USE phycst ! Define parameters for the routines
29	USE dianam ! build the name of file (routine)
30	USE ldfslp ! iso-neutral slopes
31	USE zdfmxl ! mixed layer depth
32	USE zdfddm ! ocean vertical physics: double diffusion
33	USE ioipsl ! NetCDF library
34	USE lbclnk ! ocean lateral boundary conditions (or mpp link)
35	USE diadimg ! dimg direct access file format output
36	USE trdmld_rst ! restart for diagnosing the ML trends
37	USE prtctl ! Print control
38	USE restart ! for lrst_oce
39	USE lib_mpp ! MPP library
40
41	IMPLICIT NONE
42	PRIVATE
43
44	PUBLIC trd_mld ! routine called by step.F90
45	PUBLIC trd_mld_init ! routine called by opa.F90
46	PUBLIC trd_mld_zint ! routine called by tracers routines
47
48	CHARACTER (LEN=40) :: clhstnam ! name of the trends NetCDF file
49	INTEGER :: nh_t, nmoymltrd
50	INTEGER :: nidtrd
51	INTEGER, ALLOCATABLE, SAVE, DIMENSION(:) :: ndextrd1
52	INTEGER :: ndimtrd1
53	INTEGER :: ionce, icount
54
55	!! * Control permutation of array indices
56	# include "oce_ftrans.h90"
57	# include "dom_oce_ftrans.h90"
58	# include "trdmld_oce_ftrans.h90"
59	# include "ldftra_oce_ftrans.h90"
60	# include "zdf_oce_ftrans.h90"
61	# include "ldfslp_ftrans.h90"
62	# include "zdfddm_ftrans.h90"
63
64	!! * Substitutions
65	# include "domzgr_substitute.h90"
66	# include "ldftra_substitute.h90"
67	# include "zdfddm_substitute.h90"
68	!!----------------------------------------------------------------------
69	!! NEMO/OPA 3.3 , NEMO Consortium (2010)
70	!! $Id$
71	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
72	!!----------------------------------------------------------------------
73	CONTAINS
74
75	INTEGER FUNCTION trd_mld_alloc()
76	!!----------------------------------------------------------------------
77	!! * ROUTINE trd_mld_alloc *
78	!!----------------------------------------------------------------------
79	ALLOCATE( ndextrd1(jpi*jpj) , STAT=trd_mld_alloc )
80	!
81	IF( lk_mpp ) CALL mpp_sum ( trd_mld_alloc )
82	IF( trd_mld_alloc /= 0 ) CALL ctl_warn('trd_mld_alloc: failed to allocate array ndextrd1')
83	END FUNCTION trd_mld_alloc
84
85
86	SUBROUTINE trd_mld_zint( pttrdmld, pstrdmld, ktrd, ctype )
87	!!----------------------------------------------------------------------
88	!! * ROUTINE trd_mld_zint *
89	!!
90	!! ** Purpose : Compute the vertical average of the 3D fields given as arguments
91	!! to the subroutine. This vertical average is performed from ocean
92	!! surface down to a chosen control surface.
93	!!
94	!! ** Method/usage :
95	!! The control surface can be either a mixed layer depth (time varying)
96	!! or a fixed surface (jk level or bowl).
97	!! Choose control surface with nn_ctls in namelist NAMTRD :
98	!! nn_ctls = 0 : use mixed layer with density criterion
99	!! nn_ctls = 1 : read index from file 'ctlsurf_idx'
100	!! nn_ctls > 1 : use fixed level surface jk = nn_ctls
101	!! Note: in the remainder of the routine, the volume between the
102	!! surface and the control surface is called "mixed-layer"
103	!!----------------------------------------------------------------------
104	USE wrk_nemo, ONLY: wrk_in_use, wrk_not_released
105	USE wrk_nemo, ONLY: zvlmsk => wrk_2d_10 ! 2D workspace
106	!
107	INTEGER , INTENT( in ) :: ktrd ! ocean trend index
108	CHARACTER(len=2) , INTENT( in ) :: ctype ! 2D surface/bottom or 3D interior physics
109
110	!! DCSE_NEMO: This style defeats ftrans
111	! REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( in ) :: pttrdmld ! temperature trend
112	! REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( in ) :: pstrdmld ! salinity trend
113
114	!FTRANS pttrdmld pstrdmld :I :I :z
115	REAL(wp), INTENT( in ) :: pttrdmld(jpi,jpj,jpkorig) ! temperature trend
116	REAL(wp), INTENT( in ) :: pstrdmld(jpi,jpj,jpkorig) ! salinity trend
117	!
118	INTEGER :: ji, jj, jk, isum
119	!!----------------------------------------------------------------------
120
121	IF( wrk_in_use(2, 10) ) THEN
122	CALL ctl_stop('trd_mld_zint : requested workspace arrays unavailable') ; RETURN
123	ENDIF
124
125	! I. Definition of control surface and associated fields
126	! ------------------------------------------------------
127	! ==> only once per time step <==
128
129	IF( icount == 1 ) THEN
130	!
131	tmltrd(:,:,:) = 0.e0 ; smltrd(:,:,:) = 0.e0 ! <<< reset trend arrays to zero
132
133	! ... Set nmld(ji,jj) = index of first T point below control surf. or outside mixed-layer
134	IF( nn_ctls == 0 ) THEN ! * control surface = mixed-layer with density criterion
135	nmld(:,:) = nmln(:,:) ! array nmln computed in zdfmxl.F90
136	ELSE IF( nn_ctls == 1 ) THEN ! * control surface = read index from file
137	nmld(:,:) = nbol(:,:)
138	ELSE IF( nn_ctls >= 2 ) THEN ! * control surface = model level
139	nn_ctls = MIN( nn_ctls, jpktrd - 1 )
140	nmld(:,:) = nn_ctls + 1
141	ENDIF
142
143	! ... Compute ndextrd1 and ndimtrd1 only once
144	IF( ionce == 1 ) THEN
145	!
146	! Check of validity : nmld(ji,jj) <= jpktrd
147	isum = 0
148	zvlmsk(:,:) = 0.e0
149
150	IF( jpktrd < jpk ) THEN
151	DO jj = 1, jpj
152	DO ji = 1, jpi
153	IF( nmld(ji,jj) <= jpktrd ) THEN
154	zvlmsk(ji,jj) = tmask(ji,jj,1)
155	ELSE
156	isum = isum + 1
157	zvlmsk(ji,jj) = 0.
158	ENDIF
159	END DO
160	END DO
161	ENDIF
162
163	! Index of ocean points (2D only)
164	IF( isum > 0 ) THEN
165	WRITE(numout,*)' Number of invalid points nmld > jpktrd', isum
166	CALL wheneq( jpi*jpj, zvlmsk(:,:) , 1, 1., ndextrd1, ndimtrd1 ) ! surface
167	ELSE
168	CALL wheneq( jpi*jpj, tmask(:,:,1), 1, 1., ndextrd1, ndimtrd1 ) ! surface
169	ENDIF
170
171	ionce = 0 ! no more pass here
172	!
173	END IF
174
175	! ... Weights for vertical averaging
176	wkx(:,:,:) = 0.e0
177	#if defined key_z_first
178	DO jj = 1,jpj ! initialize wkx with vertical scale factor in mixed-layer
179	DO ji = 1,jpi
180	DO jk = 1, jpktrd
181	IF( jk < nmld(ji,jj) ) wkx(ji,jj,jk) = fse3t(ji,jj,jk) * tmask(ji,jj,jk)
182	#else
183	DO jk = 1, jpktrd ! initialize wkx with vertical scale factor in mixed-layer
184	DO jj = 1,jpj
185	DO ji = 1,jpi
186	IF( jk - nmld(ji,jj) < 0.e0 ) wkx(ji,jj,jk) = fse3t(ji,jj,jk) * tmask(ji,jj,jk)
187	#endif
188	END DO
189	END DO
190	END DO
191
192	rmld(:,:) = 0.e0 ! compute mixed-layer depth : rmld
193	#if defined key_z_first
194	DO jj = 1, jpj
195	DO ji = 1, jpi
196	DO jk = 1, jpktrd
197	rmld(ji,jj) = rmld(ji,jj) + wkx(ji,jj,jk)
198	END DO
199	END DO
200	END DO
201	#else
202	DO jk = 1, jpktrd
203	rmld(:,:) = rmld(:,:) + wkx(:,:,jk)
204	END DO
205	#endif
206
207	#if defined key_z_first
208	DO jj = 1, jpj
209	DO ji = 1, jpi
210	DO jk = 1, jpktrd ! compute integration weights
211	wkx(ji,jj,jk) = wkx(ji,jj,jk) / MAX( 1., rmld(ji,jj) )
212	END DO
213	END DO
214	END DO
215	#else
216	DO jk = 1, jpktrd ! compute integration weights
217	wkx(:,:,jk) = wkx(:,:,jk) / MAX( 1., rmld(:,:) )
218	END DO
219	#endif
220
221	icount = 0 ! <<< flag = off : control surface & integr. weights
222	! ! computed only once per time step
223	END IF
224
225	! II. Vertical integration of trends in the mixed-layer
226	! -----------------------------------------------------
227
228	SELECT CASE (ctype)
229	CASE ( '3D' ) ! mean T/S trends in the mixed-layer
230	#if defined key_z_first
231	DO jj = 1, jpj
232	DO ji = 1, jpi
233	DO jk = 1, jpktrd
234	tmltrd(ji,jj,ktrd) = tmltrd(ji,jj,ktrd) + pttrdmld(ji,jj,jk) * wkx(ji,jj,jk) ! temperature
235	smltrd(ji,jj,ktrd) = smltrd(ji,jj,ktrd) + pstrdmld(ji,jj,jk) * wkx(ji,jj,jk) ! salinity
236	END DO
237	END DO
238	END DO
239	#else
240	DO jk = 1, jpktrd
241	tmltrd(:,:,ktrd) = tmltrd(:,:,ktrd) + pttrdmld(:,:,jk) * wkx(:,:,jk) ! temperature
242	smltrd(:,:,ktrd) = smltrd(:,:,ktrd) + pstrdmld(:,:,jk) * wkx(:,:,jk) ! salinity
243	END DO
244	#endif
245	CASE ( '2D' ) ! forcing at upper boundary of the mixed-layer
246	tmltrd(:,:,ktrd) = tmltrd(:,:,ktrd) + pttrdmld(:,:,1) * wkx(:,:,1) ! non penetrative
247	smltrd(:,:,ktrd) = smltrd(:,:,ktrd) + pstrdmld(:,:,1) * wkx(:,:,1)
248	END SELECT
249	!
250	IF( wrk_not_released(2, 10) ) CALL ctl_stop('trd_mld_zint: failed to release workspace arrays')
251	!
252	END SUBROUTINE trd_mld_zint
253
254	!! * Reset control of array index permutation
255	!FTRANS CLEAR
256	# include "oce_ftrans.h90"
257	# include "dom_oce_ftrans.h90"
258	# include "trdmld_oce_ftrans.h90"
259	# include "ldftra_oce_ftrans.h90"
260	# include "zdf_oce_ftrans.h90"
261	# include "ldfslp_ftrans.h90"
262	# include "zdfddm_ftrans.h90"
263
264
265	SUBROUTINE trd_mld( kt )
266	!!----------------------------------------------------------------------
267	!! * ROUTINE trd_mld *
268	!!
269	!! ** Purpose : Compute and cumulate the mixed layer trends over an analysis
270	!! period, and write NetCDF (or dimg) outputs.
271	!!
272	!! ** Method/usage :
273	!! The stored trends can be chosen twofold (according to the ln_trdmld_instant
274	!! logical namelist variable) :
275	!! 1) to explain the difference between initial and final
276	!! mixed-layer T & S (where initial and final relate to the
277	!! current analysis window, defined by nn_trd in the namelist)
278	!! 2) to explain the difference between the current and previous
279	!! TIME-AVERAGED mixed-layer T & S (where time-averaging is
280	!! performed over each analysis window).
281	!!
282	!! ** Consistency check :
283	!! If the control surface is fixed ( nn_ctls > 1 ), the residual term (dh/dt
284	!! entrainment) should be zero, at machine accuracy. Note that in the case
285	!! of time-averaged mixed-layer fields, this residual WILL NOT BE ZERO
286	!! over the first two analysis windows (except if restart).
287	!! N.B. For ORCA2_LIM, use e.g. nn_trd=5, rn_ucf=1., nn_ctls=8
288	!! for checking residuals.
289	!! On a NEC-SX5 computer, this typically leads to:
290	!! O(1.e-20) temp. residuals (tml_res) when ln_trdmld_instant=.false.
291	!! O(1.e-21) temp. residuals (tml_res) when ln_trdmld_instant=.true.
292	!!
293	!! ** Action :
294	!! At each time step, mixed-layer averaged trends are stored in the
295	!! tmltrd(:,:,jpmld_xxx) array (see trdmld_oce.F90 for definitions of jpmld_xxx).
296	!! This array is known when trd_mld is called, at the end of the stp subroutine,
297	!! except for the purely vertical K_z diffusion term, which is embedded in the
298	!! lateral diffusion trend.
299	!!
300	!! In I), this K_z term is diagnosed and stored, thus its contribution is removed
301	!! from the lateral diffusion trend.
302	!! In II), the instantaneous mixed-layer T & S are computed, and misc. cumulative
303	!! arrays are updated.
304	!! In III), called only once per analysis window, we compute the total trends,
305	!! along with the residuals and the Asselin correction terms.
306	!! In IV), the appropriate trends are written in the trends NetCDF file.
307	!!
308	!! References :
309	!! - Vialard & al.
310	!! - See NEMO documentation (in preparation)
311	!!----------------------------------------------------------------------
312	USE wrk_nemo, ONLY: wrk_in_use, wrk_not_released
313	USE wrk_nemo, ONLY: ztmltot => wrk_2d_1, zsmltot => wrk_2d_2 ! dT/dt over the anlysis window (including Asselin)
314	USE wrk_nemo, ONLY: ztmlres => wrk_2d_3, zsmlres => wrk_2d_4 ! residual = dh/dt entrainment term
315	USE wrk_nemo, ONLY: ztmlatf => wrk_2d_5, zsmlatf => wrk_2d_6 ! needed for storage only
316	USE wrk_nemo, ONLY: ztmltot2 => wrk_2d_7, ztmlres2 => wrk_2d_8, ztmltrdm2 => wrk_2d_9 ! \ working arrays to diagnose the trends
317	USE wrk_nemo, ONLY: zsmltot2 => wrk_2d_10, zsmlres2 => wrk_2d_11, zsmltrdm2 => wrk_2d_12 ! > associated with the time meaned ML T & S
318	USE wrk_nemo, ONLY: ztmlatf2 => wrk_2d_13, zsmlatf2 => wrk_2d_14
319	USE wrk_nemo, ONLY: wrk_3d_1, wrk_3d_2 ! /
320	!
321	INTEGER, INTENT( in ) :: kt ! ocean time-step index
322	!
323	INTEGER :: ji, jj, jk, jl, ik, it, itmod
324	LOGICAL :: lldebug = .TRUE.
325	REAL(wp) :: zavt, zfn, zfn2
326
327	#if defined key_z_first
328	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: ztmltrd2, zsmltrd2 ! only needed for mean diagnostics
329	#else
330	REAL(wp), POINTER, DIMENSION(:,:,:) :: ztmltrd2, zsmltrd2 ! only needed for mean diagnostics
331	#endif
332
333	#if defined key_dimgout
334	INTEGER :: iyear,imon,iday
335	CHARACTER(LEN=80) :: cltext, clmode
336	#endif
337	!!----------------------------------------------------------------------
338
339	! Check that the workspace arrays are all OK to be used
340	#if defined key_z_first
341	IF( wrk_in_use(2, 1,2,3,4,5,6,7,8,9,10,11,12,13,14) ) THEN
342	CALL ctl_stop('trd_mld : requested workspace arrays unavailable') ; RETURN
343	END IF
344	ALLOCATE( ztmltrd2(jpi,jpj,jpltrd) )
345	ALLOCATE( zsmltrd2(jpi,jpj,jpltrd) )
346	#else
347	IF( wrk_in_use(2, 1,2,3,4,5,6,7,8,9,10,11,12,13,14) .OR. &
348	wrk_in_use(3, 1,2) ) THEN
349	CALL ctl_stop('trd_mld : requested workspace arrays unavailable') ; RETURN
350	ELSE IF(jpltrd > jpk) THEN
351	! ARPDBG, is this reasonable or will this cause trouble in the future?
352	CALL ctl_stop('trd_mld : no. of mixed-layer trends (jpltrd) exceeds no. of model levels so cannot use 3D workspaces.')
353	RETURN
354	END IF
355	! Set-up pointers into sub-arrays of 3d-workspaces
356	ztmltrd2 => wrk_3d_1(1:,:,1:jpltrd)
357	zsmltrd2 => wrk_3d_2(1:,:,1:jpltrd)
358	#endif
359
360	! ======================================================================
361	! I. Diagnose the purely vertical (K_z) diffusion trend
362	! ======================================================================
363
364	! ... These terms can be estimated by flux computation at the lower boundary of the ML
365	! (we compute (-1/h) * K_z * d_z( T ) and (-1/h) * K_z * d_z( S ))
366	IF( ln_traldf_iso ) THEN
367	DO jj = 1,jpj
368	DO ji = 1,jpi
369	ik = nmld(ji,jj)
370	zavt = avt(ji,jj,ik)
371	tmltrd(ji,jj,jpmld_zdf) = - zavt / fse3w(ji,jj,ik) * tmask(ji,jj,ik) &
372	& * ( tn(ji,jj,ik-1) - tn(ji,jj,ik) ) &
373	& / MAX( 1., rmld(ji,jj) ) * tmask(ji,jj,1)
374	zavt = fsavs(ji,jj,ik)
375	smltrd(ji,jj,jpmld_zdf) = - zavt / fse3w(ji,jj,ik) * tmask(ji,jj,ik) &
376	& * ( sn(ji,jj,ik-1) - sn(ji,jj,ik) ) &
377	& / MAX( 1., rmld(ji,jj) ) * tmask(ji,jj,1)
378	END DO
379	END DO
380
381	! ... Remove this K_z trend from the iso-neutral diffusion term (if any)
382	tmltrd(:,:,jpmld_ldf) = tmltrd(:,:,jpmld_ldf) - tmltrd(:,:,jpmld_zdf)
383	smltrd(:,:,jpmld_ldf) = smltrd(:,:,jpmld_ldf) - smltrd(:,:,jpmld_zdf)
384	END IF
385
386	! ... Lateral boundary conditions
387	DO jl = 1, jpltrd
388	CALL lbc_lnk( tmltrd(:,:,jl), 'T', 1. )
389	CALL lbc_lnk( smltrd(:,:,jl), 'T', 1. )
390	END DO
391
392	! ======================================================================
393	! II. Cumulate the trends over the analysis window
394	! ======================================================================
395
396	ztmltrd2(:,:,:) = 0.e0 ; zsmltrd2(:,:,:) = 0.e0 ! <<< reset arrays to zero
397	ztmltot2(:,:) = 0.e0 ; zsmltot2(:,:) = 0.e0
398	ztmlres2(:,:) = 0.e0 ; zsmlres2(:,:) = 0.e0
399	ztmlatf2(:,:) = 0.e0 ; zsmlatf2(:,:) = 0.e0
400
401	! II.1 Set before values of vertically average T and S
402	! ----------------------------------------------------
403	IF( kt > nit000 ) THEN
404	! ... temperature ... ... salinity ...
405	tmlb (:,:) = tml (:,:) ; smlb (:,:) = sml (:,:)
406	tmlatfn(:,:) = tmltrd(:,:,jpmld_atf) ; smlatfn(:,:) = smltrd(:,:,jpmld_atf)
407	END IF
408
409	! II.2 Vertically averaged T and S
410	! --------------------------------
411	tml(:,:) = 0.e0 ; sml(:,:) = 0.e0
412	#if defined key_z_first
413	DO jj = 1, jpj
414	DO ji = 1, jpi
415	DO jk = 1, jpktrd - 1
416	tml(ji,jj) = tml(ji,jj) + wkx(ji,jj,jk) * tn(ji,jj,jk)
417	sml(ji,jj) = sml(ji,jj) + wkx(ji,jj,jk) * sn(ji,jj,jk)
418	END DO
419	END DO
420	END DO
421	#else
422	DO jk = 1, jpktrd - 1
423	tml(:,:) = tml(:,:) + wkx(:,:,jk) * tn(:,:,jk)
424	sml(:,:) = sml(:,:) + wkx(:,:,jk) * sn(:,:,jk)
425	END DO
426	#endif
427
428	! II.3 Initialize mixed-layer "before" arrays for the 1rst analysis window
429	! ------------------------------------------------------------------------
430	IF( kt == 2 ) THEN ! i.e. ( .NOT. ln_rstart ).AND.( kt == nit000 + 1)
431	!
432	! ... temperature ... ... salinity ...
433	tmlbb (:,:) = tmlb (:,:) ; smlbb (:,:) = smlb (:,:)
434	tmlbn (:,:) = tml (:,:) ; smlbn (:,:) = sml (:,:)
435	tmlatfb(:,:) = tmlatfn(:,:) ; smlatfb(:,:) = smlatfn(:,:)
436
437	tmltrd_csum_ub (:,:,:) = 0.e0 ; smltrd_csum_ub (:,:,:) = 0.e0
438	tmltrd_atf_sumb(:,:) = 0.e0 ; smltrd_atf_sumb(:,:) = 0.e0
439
440	rmldbn(:,:) = rmld(:,:)
441
442	IF( ln_ctl ) THEN
443	WRITE(numout,*) ' we reach kt == nit000 + 1 = ', nit000+1
444	CALL prt_ctl(tab2d_1=tmlbb , clinfo1=' tmlbb - : ', mask1=tmask, ovlap=1)
445	CALL prt_ctl(tab2d_1=tmlbn , clinfo1=' tmlbn - : ', mask1=tmask, ovlap=1)
446	CALL prt_ctl(tab2d_1=tmlatfb , clinfo1=' tmlatfb - : ', mask1=tmask, ovlap=1)
447	END IF
448	!
449	END IF
450
451	IF( ( ln_rstart ) .AND. ( kt == nit000 ) .AND. ( ln_ctl ) ) THEN
452	IF( ln_trdmld_instant ) THEN
453	WRITE(numout,*) ' restart from kt == nit000 = ', nit000
454	CALL prt_ctl(tab2d_1=tmlbb , clinfo1=' tmlbb - : ', mask1=tmask, ovlap=1)
455	CALL prt_ctl(tab2d_1=tmlbn , clinfo1=' tmlbn - : ', mask1=tmask, ovlap=1)
456	CALL prt_ctl(tab2d_1=tmlatfb , clinfo1=' tmlatfb - : ', mask1=tmask, ovlap=1)
457	ELSE
458	WRITE(numout,*) ' restart from kt == nit000 = ', nit000
459	CALL prt_ctl(tab2d_1=tmlbn , clinfo1=' tmlbn - : ', mask1=tmask, ovlap=1)
460	CALL prt_ctl(tab2d_1=rmldbn , clinfo1=' rmldbn - : ', mask1=tmask, ovlap=1)
461	CALL prt_ctl(tab2d_1=tml_sumb , clinfo1=' tml_sumb - : ', mask1=tmask, ovlap=1)
462	CALL prt_ctl(tab2d_1=tmltrd_atf_sumb, clinfo1=' tmltrd_atf_sumb - : ', mask1=tmask, ovlap=1)
463	CALL prt_ctl(tab3d_1=tmltrd_csum_ub , clinfo1=' tmltrd_csum_ub - : ', mask1=tmask, ovlap=1, kdim=1)
464	END IF
465	END IF
466
467	! II.4 Cumulated trends over the analysis period
468	! ----------------------------------------------
469	!
470	! [ 1rst analysis window ] [ 2nd analysis window ]
471	!
472	! o---[--o-----o-----o-----o--]-[--o-----o-----o-----o-----o--]---o-----o--> time steps
473	! nn_trd 2*nn_trd etc.
474	! 1 2 3 4 =5 e.g. =10
475	!
476	IF( ( kt >= 2 ).OR.( ln_rstart ) ) THEN
477	!
478	nmoymltrd = nmoymltrd + 1
479
480	! ... Cumulate over BOTH physical contributions AND over time steps
481	DO jl = 1, jpltrd
482	tmltrdm(:,:) = tmltrdm(:,:) + tmltrd(:,:,jl)
483	smltrdm(:,:) = smltrdm(:,:) + smltrd(:,:,jl)
484	END DO
485
486	! ... Special handling of the Asselin trend
487	tmlatfm(:,:) = tmlatfm(:,:) + tmlatfn(:,:)
488	smlatfm(:,:) = smlatfm(:,:) + smlatfn(:,:)
489
490	! ... Trends associated with the time mean of the ML T/S
491	tmltrd_sum (:,:,:) = tmltrd_sum (:,:,:) + tmltrd (:,:,:) ! tem
492	tmltrd_csum_ln(:,:,:) = tmltrd_csum_ln(:,:,:) + tmltrd_sum(:,:,:)
493	tml_sum (:,:) = tml_sum (:,:) + tml (:,:)
494	smltrd_sum (:,:,:) = smltrd_sum (:,:,:) + smltrd (:,:,:) ! sal
495	smltrd_csum_ln(:,:,:) = smltrd_csum_ln(:,:,:) + smltrd_sum(:,:,:)
496	sml_sum (:,:) = sml_sum (:,:) + sml (:,:)
497	rmld_sum (:,:) = rmld_sum (:,:) + rmld (:,:) ! rmld
498	!
499	END IF
500
501	! ======================================================================
502	! III. Prepare fields for output (get here ONCE PER ANALYSIS PERIOD)
503	! ======================================================================
504
505	! Convert to appropriate physical units
506	! N.B. It may be useful to check IOIPSL time averaging with :
507	! tmltrd (:,:,:) = 1. ; smltrd (:,:,:) = 1.
508	tmltrd(:,:,:) = tmltrd(:,:,:) * rn_ucf ! (actually needed for 1:jpltrd-1, but trdmld(:,:,jpltrd)
509	smltrd(:,:,:) = smltrd(:,:,:) * rn_ucf ! is no longer used, and is reset to 0. at next time step)
510
511	! define time axis
512	it = kt
513	itmod = kt - nit000 + 1
514
515	MODULO_NTRD : IF( MOD( itmod, nn_trd ) == 0 ) THEN ! nitend MUST be multiple of nn_trd
516	!
517	ztmltot (:,:) = 0.e0 ; zsmltot (:,:) = 0.e0 ! reset arrays to zero
518	ztmlres (:,:) = 0.e0 ; zsmlres (:,:) = 0.e0
519	ztmltot2(:,:) = 0.e0 ; zsmltot2(:,:) = 0.e0
520	ztmlres2(:,:) = 0.e0 ; zsmlres2(:,:) = 0.e0
521
522	zfn = float(nmoymltrd) ; zfn2 = zfn * zfn
523
524	! III.1 Prepare fields for output ("instantaneous" diagnostics)
525	! -------------------------------------------------------------
526
527	!-- Compute total trends
528	ztmltot(:,:) = ( tml(:,:) - tmlbn(:,:) + tmlb(:,:) - tmlbb(:,:) ) / ( 2.*rdt )
529	zsmltot(:,:) = ( sml(:,:) - smlbn(:,:) + smlb(:,:) - smlbb(:,:) ) / ( 2.*rdt )
530
531	!-- Compute residuals
532	ztmlres(:,:) = ztmltot(:,:) - ( tmltrdm(:,:) - tmlatfn(:,:) + tmlatfb(:,:) )
533	zsmlres(:,:) = zsmltot(:,:) - ( smltrdm(:,:) - smlatfn(:,:) + smlatfb(:,:) )
534
535	!-- Diagnose Asselin trend over the analysis window
536	ztmlatf(:,:) = tmlatfm(:,:) - tmlatfn(:,:) + tmlatfb(:,:)
537	zsmlatf(:,:) = smlatfm(:,:) - smlatfn(:,:) + smlatfb(:,:)
538
539	!-- Lateral boundary conditions
540	! ... temperature ... ... salinity ...
541	CALL lbc_lnk( ztmltot , 'T', 1. ) ; CALL lbc_lnk( zsmltot , 'T', 1. )
542	CALL lbc_lnk( ztmlres , 'T', 1. ) ; CALL lbc_lnk( zsmlres , 'T', 1. )
543	CALL lbc_lnk( ztmlatf , 'T', 1. ) ; CALL lbc_lnk( zsmlatf , 'T', 1. )
544
545	#if defined key_diainstant
546	CALL ctl_stop( 'tml_trd : key_diainstant was never checked within trdmld. Comment this to proceed.')
547	#endif
548	! III.2 Prepare fields for output ("mean" diagnostics)
549	! ----------------------------------------------------
550
551	!-- Update the ML depth time sum (to build the Leap-Frog time mean)
552	rmld_sum(:,:) = rmldbn(:,:) + 2 * ( rmld_sum(:,:) - rmld(:,:) ) + rmld(:,:)
553
554	!-- Compute temperature total trends
555	tml_sum (:,:) = tmlbn(:,:) + 2 * ( tml_sum(:,:) - tml(:,:) ) + tml(:,:)
556	ztmltot2(:,:) = ( tml_sum(:,:) - tml_sumb(:,:) ) / ( 2.*rdt ) ! now in degC/s
557
558	!-- Compute salinity total trends
559	sml_sum (:,:) = smlbn(:,:) + 2 * ( sml_sum(:,:) - sml(:,:) ) + sml(:,:)
560	zsmltot2(:,:) = ( sml_sum(:,:) - sml_sumb(:,:) ) / ( 2.*rdt ) ! now in psu/s
561
562	!-- Compute temperature residuals
563	DO jl = 1, jpltrd
564	ztmltrd2(:,:,jl) = tmltrd_csum_ub(:,:,jl) + tmltrd_csum_ln(:,:,jl)
565	END DO
566
567	ztmltrdm2(:,:) = 0.e0
568	DO jl = 1, jpltrd
569	ztmltrdm2(:,:) = ztmltrdm2(:,:) + ztmltrd2(:,:,jl)
570	END DO
571
572	ztmlres2(:,:) = ztmltot2(:,:) - &
573	( ztmltrdm2(:,:) - tmltrd_sum(:,:,jpmld_atf) + tmltrd_atf_sumb(:,:) )
574
575	!-- Compute salinity residuals
576	DO jl = 1, jpltrd
577	zsmltrd2(:,:,jl) = smltrd_csum_ub(:,:,jl) + smltrd_csum_ln(:,:,jl)
578	END DO
579
580	zsmltrdm2(:,:) = 0.
581	DO jl = 1, jpltrd
582	zsmltrdm2(:,:) = zsmltrdm2(:,:) + zsmltrd2(:,:,jl)
583	END DO
584
585	zsmlres2(:,:) = zsmltot2(:,:) - &
586	( zsmltrdm2(:,:) - smltrd_sum(:,:,jpmld_atf) + smltrd_atf_sumb(:,:) )
587
588	!-- Diagnose Asselin trend over the analysis window
589	ztmlatf2(:,:) = ztmltrd2(:,:,jpmld_atf) - tmltrd_sum(:,:,jpmld_atf) + tmltrd_atf_sumb(:,:)
590	zsmlatf2(:,:) = zsmltrd2(:,:,jpmld_atf) - smltrd_sum(:,:,jpmld_atf) + smltrd_atf_sumb(:,:)
591
592	!-- Lateral boundary conditions
593	! ... temperature ... ... salinity ...
594	CALL lbc_lnk( ztmltot2, 'T', 1. ) ; CALL lbc_lnk( zsmltot2, 'T', 1. )
595	CALL lbc_lnk( ztmlres2, 'T', 1. ) ; CALL lbc_lnk( zsmlres2, 'T', 1. )
596	DO jl = 1, jpltrd
597	CALL lbc_lnk( ztmltrd2(:,:,jl), 'T', 1. ) ! \ these will be output
598	CALL lbc_lnk( zsmltrd2(:,:,jl), 'T', 1. ) ! / in the NetCDF trends file
599	END DO
600
601	! III.3 Time evolution array swap
602	! -------------------------------
603
604	! For T/S instantaneous diagnostics
605	! ... temperature ... ... salinity ...
606	tmlbb (:,:) = tmlb (:,:) ; smlbb (:,:) = smlb (:,:)
607	tmlbn (:,:) = tml (:,:) ; smlbn (:,:) = sml (:,:)
608	tmlatfb(:,:) = tmlatfn(:,:) ; smlatfb(:,:) = smlatfn(:,:)
609
610	! For T mean diagnostics
611	tmltrd_csum_ub (:,:,:) = zfn * tmltrd_sum(:,:,:) - tmltrd_csum_ln(:,:,:)
612	tml_sumb (:,:) = tml_sum(:,:)
613	tmltrd_atf_sumb(:,:) = tmltrd_sum(:,:,jpmld_atf)
614
615	! For S mean diagnostics
616	smltrd_csum_ub (:,:,:) = zfn * smltrd_sum(:,:,:) - smltrd_csum_ln(:,:,:)
617	sml_sumb (:,:) = sml_sum(:,:)
618	smltrd_atf_sumb(:,:) = smltrd_sum(:,:,jpmld_atf)
619
620	! ML depth
621	rmldbn (:,:) = rmld (:,:)
622
623	IF( ln_ctl ) THEN
624	IF( ln_trdmld_instant ) THEN
625	CALL prt_ctl(tab2d_1=tmlbb , clinfo1=' tmlbb - : ', mask1=tmask, ovlap=1)
626	CALL prt_ctl(tab2d_1=tmlbn , clinfo1=' tmlbn - : ', mask1=tmask, ovlap=1)
627	CALL prt_ctl(tab2d_1=tmlatfb , clinfo1=' tmlatfb - : ', mask1=tmask, ovlap=1)
628	ELSE
629	CALL prt_ctl(tab2d_1=tmlbn , clinfo1=' tmlbn - : ', mask1=tmask, ovlap=1)
630	CALL prt_ctl(tab2d_1=rmldbn , clinfo1=' rmldbn - : ', mask1=tmask, ovlap=1)
631	CALL prt_ctl(tab2d_1=tml_sumb , clinfo1=' tml_sumb - : ', mask1=tmask, ovlap=1)
632	CALL prt_ctl(tab2d_1=tmltrd_atf_sumb, clinfo1=' tmltrd_atf_sumb - : ', mask1=tmask, ovlap=1)
633	CALL prt_ctl(tab3d_1=tmltrd_csum_ub , clinfo1=' tmltrd_csum_ub - : ', mask1=tmask, ovlap=1, kdim=1)
634	END IF
635	END IF
636
637	! III.4 Convert to appropriate physical units
638	! -------------------------------------------
639
640	! ... temperature ... ... salinity ...
641	ztmltot (:,:) = ztmltot(:,:) * rn_ucf/zfn ; zsmltot (:,:) = zsmltot(:,:) * rn_ucf/zfn
642	ztmlres (:,:) = ztmlres(:,:) * rn_ucf/zfn ; zsmlres (:,:) = zsmlres(:,:) * rn_ucf/zfn
643	ztmlatf (:,:) = ztmlatf(:,:) * rn_ucf/zfn ; zsmlatf (:,:) = zsmlatf(:,:) * rn_ucf/zfn
644
645	tml_sum (:,:) = tml_sum (:,:) / (2zfn) ; sml_sum (:,:) = sml_sum (:,:) / (2zfn)
646	ztmltot2(:,:) = ztmltot2(:,:) * rn_ucf/zfn2 ; zsmltot2(:,:) = zsmltot2(:,:) * rn_ucf/zfn2
647	ztmltrd2(:,:,:) = ztmltrd2(:,:,:)* rn_ucf/zfn2 ; zsmltrd2(:,:,:) = zsmltrd2(:,:,:)* rn_ucf/zfn2
648	ztmlatf2(:,:) = ztmlatf2(:,:) * rn_ucf/zfn2 ; zsmlatf2(:,:) = zsmlatf2(:,:) * rn_ucf/zfn2
649	ztmlres2(:,:) = ztmlres2(:,:) * rn_ucf/zfn2 ; zsmlres2(:,:) = zsmlres2(:,:) * rn_ucf/zfn2
650
651	rmld_sum(:,:) = rmld_sum(:,:) / (2*zfn) ! similar to tml_sum and sml_sum
652
653	! * Debugging information *
654	IF( lldebug ) THEN
655	!
656	WRITE(numout,*)
657	WRITE(numout,*) 'trd_mld : write trends in the Mixed Layer for debugging process:'
658	WRITE(numout,*) '~~~~~~~ '
659	WRITE(numout,*) ' TRA kt = ', kt, 'nmoymltrd = ', nmoymltrd
660	WRITE(numout,*)
661	WRITE(numout,*) ' >>>>>>>>>>>>>>>>>> TRA TEMPERATURE <<<<<<<<<<<<<<<<<<'
662	WRITE(numout,*) ' TRA ztmlres : ', SUM(ztmlres(:,:))
663	WRITE(numout,*) ' TRA ztmltot : ', SUM(ztmltot(:,:))
664	WRITE(numout,*) ' TRA tmltrdm : ', SUM(tmltrdm(:,:))
665	WRITE(numout,*) ' TRA tmlatfb : ', SUM(tmlatfb(:,:))
666	WRITE(numout,*) ' TRA tmlatfn : ', SUM(tmlatfn(:,:))
667	DO jl = 1, jpltrd
668	WRITE(numout,) ' TRA TREND INDEX jpmld_xxx = jl = ', jl, &
669	& ' tmltrd : ', SUM(tmltrd(:,:,jl))
670	END DO
671	WRITE(numout,*) ' TRA ztmlres (jpi/2,jpj/2) : ', ztmlres (jpi/2,jpj/2)
672	WRITE(numout,*) ' TRA ztmlres2(jpi/2,jpj/2) : ', ztmlres2(jpi/2,jpj/2)
673	WRITE(numout,*)
674	WRITE(numout,*) ' >>>>>>>>>>>>>>>>>> TRA SALINITY <<<<<<<<<<<<<<<<<<'
675	WRITE(numout,*) ' TRA zsmlres : ', SUM(zsmlres(:,:))
676	WRITE(numout,*) ' TRA zsmltot : ', SUM(zsmltot(:,:))
677	WRITE(numout,*) ' TRA smltrdm : ', SUM(smltrdm(:,:))
678	WRITE(numout,*) ' TRA smlatfb : ', SUM(smlatfb(:,:))
679	WRITE(numout,*) ' TRA smlatfn : ', SUM(smlatfn(:,:))
680	DO jl = 1, jpltrd
681	WRITE(numout,) ' TRA TREND INDEX jpmld_xxx = jl = ', jl, &
682	& ' smltrd : ', SUM(smltrd(:,:,jl))
683	END DO
684	WRITE(numout,*) ' TRA zsmlres (jpi/2,jpj/2) : ', zsmlres (jpi/2,jpj/2)
685	WRITE(numout,*) ' TRA zsmlres2(jpi/2,jpj/2) : ', zsmlres2(jpi/2,jpj/2)
686	!
687	END IF
688	!
689	END IF MODULO_NTRD
690
691	! ======================================================================
692	! IV. Write trends in the NetCDF file
693	! ======================================================================
694
695	! IV.1 Code for dimg mpp output
696	! -----------------------------
697
698	#if defined key_dimgout
699
700	IF( MOD( itmod, nn_trd ) == 0 ) THEN
701	iyear = ndastp/10000
702	imon = (ndastp-iyear*10000)/100
703	iday = ndastp - imon100 - iyear10000
704	WRITE(clname,9000) TRIM(cexper),'MLDiags',iyear,imon,iday
705	WRITE(clmode,'(f5.1,a)') nn_trd*rdt/86400.,' days average'
706	cltext = TRIM(cexper)//' mld diags'//TRIM(clmode)
707	CALL dia_wri_dimg (clname, cltext, smltrd, jpltrd, '2')
708	END IF
709
710	9000 FORMAT(a,"_",a,"_y",i4.4,"m",i2.2,"d",i2.2,".dimgproc")
711
712	#else
713
714	! IV.2 Code for IOIPSL/NetCDF output
715	! ----------------------------------
716
717	IF( lwp .AND. MOD( itmod , nn_trd ) == 0 ) THEN
718	WRITE(numout,*) ' '
719	WRITE(numout,*) 'trd_mld : write trends in the NetCDF file :'
720	WRITE(numout,*) '~~~~~~~ '
721	WRITE(numout,*) ' ', TRIM(clhstnam), ' at kt = ', kt
722	WRITE(numout,*) ' N.B. nmoymltrd = ', nmoymltrd
723	WRITE(numout,*) ' '
724	END IF
725
726	!-- Write the trends for T/S instantaneous diagnostics
727	IF( ln_trdmld_instant ) THEN
728
729	CALL histwrite( nidtrd, "mxl_depth", it, rmld(:,:), ndimtrd1, ndextrd1 )
730
731	!................................. ( ML temperature ) ...................................
732
733	!-- Output the fields
734	CALL histwrite( nidtrd, "tml" , it, tml (:,:), ndimtrd1, ndextrd1 )
735	CALL histwrite( nidtrd, "tml_tot" , it, ztmltot(:,:), ndimtrd1, ndextrd1 )
736	CALL histwrite( nidtrd, "tml_res" , it, ztmlres(:,:), ndimtrd1, ndextrd1 )
737
738	DO jl = 1, jpltrd - 1
739	CALL histwrite( nidtrd, trim("tml"//ctrd(jl,2)), &
740	& it, tmltrd (:,:,jl), ndimtrd1, ndextrd1 )
741	END DO
742
743	CALL histwrite( nidtrd, trim("tml"//ctrd(jpmld_atf,2)), &
744	& it, ztmlatf(:,:), ndimtrd1, ndextrd1 )
745
746	!.................................. ( ML salinity ) .....................................
747
748	!-- Output the fields
749	CALL histwrite( nidtrd, "sml" , it, sml (:,:), ndimtrd1, ndextrd1 )
750	CALL histwrite( nidtrd, "sml_tot" , it, zsmltot(:,:), ndimtrd1, ndextrd1 )
751	CALL histwrite( nidtrd, "sml_res" , it, zsmlres(:,:), ndimtrd1, ndextrd1 )
752
753	DO jl = 1, jpltrd - 1
754	CALL histwrite( nidtrd, trim("sml"//ctrd(jl,2)), &
755	& it, smltrd(:,:,jl), ndimtrd1, ndextrd1 )
756	END DO
757
758	CALL histwrite( nidtrd, trim("sml"//ctrd(jpmld_atf,2)), &
759	& it, zsmlatf(:,:), ndimtrd1, ndextrd1 )
760
761	IF( kt == nitend ) CALL histclo( nidtrd )
762
763	!-- Write the trends for T/S mean diagnostics
764	ELSE
765
766	CALL histwrite( nidtrd, "mxl_depth", it, rmld_sum(:,:), ndimtrd1, ndextrd1 )
767
768	!................................. ( ML temperature ) ...................................
769
770	!-- Output the fields
771	CALL histwrite( nidtrd, "tml" , it, tml_sum (:,:), ndimtrd1, ndextrd1 )
772	CALL histwrite( nidtrd, "tml_tot" , it, ztmltot2(:,:), ndimtrd1, ndextrd1 )
773	CALL histwrite( nidtrd, "tml_res" , it, ztmlres2(:,:), ndimtrd1, ndextrd1 )
774
775	DO jl = 1, jpltrd - 1
776	CALL histwrite( nidtrd, trim("tml"//ctrd(jl,2)), &
777	& it, ztmltrd2(:,:,jl), ndimtrd1, ndextrd1 )
778	END DO
779
780	CALL histwrite( nidtrd, trim("tml"//ctrd(jpmld_atf,2)), &
781	& it, ztmlatf2(:,:), ndimtrd1, ndextrd1 )
782
783	!.................................. ( ML salinity ) .....................................
784
785	!-- Output the fields
786	CALL histwrite( nidtrd, "sml" , it, sml_sum (:,:), ndimtrd1, ndextrd1 )
787	CALL histwrite( nidtrd, "sml_tot" , it, zsmltot2(:,:), ndimtrd1, ndextrd1 )
788	CALL histwrite( nidtrd, "sml_res" , it, zsmlres2(:,:), ndimtrd1, ndextrd1 )
789
790	DO jl = 1, jpltrd - 1
791	CALL histwrite( nidtrd, trim("sml"//ctrd(jl,2)), &
792	& it, zsmltrd2(:,:,jl), ndimtrd1, ndextrd1 )
793	END DO
794
795	CALL histwrite( nidtrd, trim("sml"//ctrd(jpmld_atf,2)), &
796	& it, zsmlatf2(:,:), ndimtrd1, ndextrd1 )
797
798	IF( kt == nitend ) CALL histclo( nidtrd )
799
800	END IF
801
802	! Compute the control surface (for next time step) : flag = on
803	icount = 1
804	!
805	#endif
806
807	IF( MOD( itmod, nn_trd ) == 0 ) THEN
808	!
809	! III.5 Reset cumulative arrays to zero
810	! -------------------------------------
811	nmoymltrd = 0
812
813	! ... temperature ... ... salinity ...
814	tmltrdm (:,:) = 0.e0 ; smltrdm (:,:) = 0.e0
815	tmlatfm (:,:) = 0.e0 ; smlatfm (:,:) = 0.e0
816	tml_sum (:,:) = 0.e0 ; sml_sum (:,:) = 0.e0
817	tmltrd_csum_ln (:,:,:) = 0.e0 ; smltrd_csum_ln (:,:,:) = 0.e0
818	tmltrd_sum (:,:,:) = 0.e0 ; smltrd_sum (:,:,:) = 0.e0
819
820	rmld_sum (:,:) = 0.e0
821	!
822	END IF
823
824	! ======================================================================
825	! V. Write restart file
826	! ======================================================================
827
828	IF( lrst_oce ) CALL trd_mld_rst_write( kt )
829
830	#if defined key_z_first
831	IF( wrk_not_released(2, 1,2,3,4,5,6,7,8,9,10,11,12,13,14) ) &
832	CALL ctl_stop('trd_mld : failed to release workspace arrays.')
833	DEALLOCATE( ztmltrd2, zsmltrd2 )
834	#else
835	IF( wrk_not_released(2, 1,2,3,4,5,6,7,8,9,10,11,12,13,14) .OR. &
836	wrk_not_released(3, 1,2) ) &
837	CALL ctl_stop('trd_mld : failed to release workspace arrays.')
838	#endif
839	!
840	END SUBROUTINE trd_mld
841
842
843	SUBROUTINE trd_mld_init
844	!!----------------------------------------------------------------------
845	!! * ROUTINE trd_mld_init *
846	!!
847	!! ** Purpose : computation of vertically integrated T and S budgets
848	!! from ocean surface down to control surface (NetCDF output)
849	!!----------------------------------------------------------------------
850	INTEGER :: jl
851	INTEGER :: inum ! logical unit
852	REAL(wp) :: zjulian, zsto, zout
853	CHARACTER (LEN=40) :: clop
854	CHARACTER (LEN=12) :: clmxl, cltu, clsu
855	!!----------------------------------------------------------------------
856
857	! ======================================================================
858	! I. initialization
859	! ======================================================================
860
861	IF(lwp) THEN
862	WRITE(numout,*)
863	WRITE(numout,*) ' trd_mld_init : Mixed-layer trends'
864	WRITE(numout,*) ' ~~~~~~~~~~~~~'
865	WRITE(numout,*) ' namelist namtrd read in trd_mod_init '
866	WRITE(numout,*)
867	END IF
868
869	! I.1 Check consistency of user defined preferences
870	! -------------------------------------------------
871
872	IF( ( lk_trdmld ) .AND. ( MOD( nitend, nn_trd ) /= 0 ) ) THEN
873	WRITE(numout,cform_err)
874	WRITE(numout,*) ' Your nitend parameter, nitend = ', nitend
875	WRITE(numout,*) ' is no multiple of the trends diagnostics frequency '
876	WRITE(numout,*) ' you defined, nn_trd = ', nn_trd
877	WRITE(numout,*) ' This will not allow you to restart from this simulation. '
878	WRITE(numout,*) ' You should reconsider this choice. '
879	WRITE(numout,*)
880	WRITE(numout,*) ' N.B. the nitend parameter is also constrained to be a '
881	WRITE(numout,*) ' multiple of the sea-ice frequency parameter (typically 5) '
882	nstop = nstop + 1
883	END IF
884
885	IF( ( lk_trdmld ) .AND. ( nn_cla == 1 ) ) THEN
886	WRITE(numout,cform_war)
887	WRITE(numout,*) ' You set n_cla = 1. Note that the Mixed-Layer diagnostics '
888	WRITE(numout,*) ' are not exact along the corresponding straits. '
889	nwarn = nwarn + 1
890	END IF
891
892	! ! allocate trdmld arrays
893	IF( trd_mld_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'trd_mld_init : unable to allocate trdmld arrays' )
894	IF( trdmld_oce_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'trd_mld_init : unable to allocate trdmld_oce arrays' )
895
896	! I.2 Initialize arrays to zero or read a restart file
897	! ----------------------------------------------------
898
899	nmoymltrd = 0
900
901	! ... temperature ... ... salinity ...
902	tml (:,:) = 0.e0 ; sml (:,:) = 0.e0 ! inst.
903	tmltrdm (:,:) = 0.e0 ; smltrdm (:,:) = 0.e0
904	tmlatfm (:,:) = 0.e0 ; smlatfm (:,:) = 0.e0
905	tml_sum (:,:) = 0.e0 ; sml_sum (:,:) = 0.e0 ! mean
906	tmltrd_sum (:,:,:) = 0.e0 ; smltrd_sum (:,:,:) = 0.e0
907	tmltrd_csum_ln (:,:,:) = 0.e0 ; smltrd_csum_ln (:,:,:) = 0.e0
908
909	rmld (:,:) = 0.e0
910	rmld_sum (:,:) = 0.e0
911
912	IF( ln_rstart .AND. ln_trdmld_restart ) THEN
913	CALL trd_mld_rst_read
914	ELSE
915	! ... temperature ... ... salinity ...
916	tmlb (:,:) = 0.e0 ; smlb (:,:) = 0.e0 ! inst.
917	tmlbb (:,:) = 0.e0 ; smlbb (:,:) = 0.e0
918	tmlbn (:,:) = 0.e0 ; smlbn (:,:) = 0.e0
919	tml_sumb (:,:) = 0.e0 ; sml_sumb (:,:) = 0.e0 ! mean
920	tmltrd_csum_ub (:,:,:) = 0.e0 ; smltrd_csum_ub (:,:,:) = 0.e0
921	tmltrd_atf_sumb(:,:) = 0.e0 ; smltrd_atf_sumb(:,:) = 0.e0
922	END IF
923
924	icount = 1 ; ionce = 1 ! open specifier
925
926	! I.3 Read control surface from file ctlsurf_idx
927	! ----------------------------------------------
928
929	IF( nn_ctls == 1 ) THEN
930	CALL ctl_opn( inum, 'ctlsurf_idx', 'OLD', 'UNFORMATTED', 'SEQUENTIAL', -1, numout, lwp )
931	READ ( inum ) nbol
932	CLOSE( inum )
933	END IF
934
935	! ======================================================================
936	! II. netCDF output initialization
937	! ======================================================================
938
939	#if defined key_dimgout
940	???
941	#else
942	! clmxl = legend root for netCDF output
943	IF( nn_ctls == 0 ) THEN ! control surface = mixed-layer with density criterion
944	clmxl = 'Mixed Layer ' ! (array nmln computed in zdfmxl.F90)
945	ELSE IF( nn_ctls == 1 ) THEN ! control surface = read index from file
946	clmxl = ' Bowl '
947	ELSE IF( nn_ctls >= 2 ) THEN ! control surface = model level
948	WRITE(clmxl,'(A10,I2,1X)') 'Levels 1 -', nn_ctls
949	END IF
950
951	! II.1 Define frequency of output and means
952	! -----------------------------------------
953	IF( ln_mskland ) THEN ; clop = "only(x)" ! put 1.e+20 on land (very expensive!!)
954	ELSE ; clop = "x" ! no use of the mask value (require less cpu time)
955	ENDIF
956	# if defined key_diainstant
957	IF( .NOT. ln_trdmld_instant ) THEN
958	CALL ctl_stop( 'trd_mld : this was never checked. Comment this line to proceed...' )
959	END IF
960	zsto = nn_trd * rdt
961	clop = "inst("//TRIM(clop)//")"
962	# else
963	IF( ln_trdmld_instant ) THEN
964	zsto = rdt ! inst. diags : we use IOIPSL time averaging
965	ELSE
966	zsto = nn_trd * rdt ! mean diags : we DO NOT use any IOIPSL time averaging
967	END IF
968	clop = "ave("//TRIM(clop)//")"
969	# endif
970	zout = nn_trd * rdt
971
972	IF(lwp) WRITE (numout,*) ' netCDF initialization'
973
974	! II.2 Compute julian date from starting date of the run
975	! ------------------------------------------------------
976	CALL ymds2ju( nyear, nmonth, nday, rdt, zjulian )
977	zjulian = zjulian - adatrj ! set calendar origin to the beginning of the experiment
978	IF(lwp) WRITE(numout,*)' '
979	IF(lwp) WRITE(numout,*)' Date 0 used :',nit000, &
980	& ' YEAR ', nyear,' MONTH ' , nmonth, &
981	& ' DAY ' , nday, 'Julian day : ', zjulian
982
983
984	! II.3 Define the T grid trend file (nidtrd)
985	! ------------------------------------------
986	!-- Define long and short names for the NetCDF output variables
987	! ==> choose them according to trdmld_oce.F90 <==
988
989	ctrd(jpmld_xad,1) = " Zonal advection" ; ctrd(jpmld_xad,2) = "_xad"
990	ctrd(jpmld_yad,1) = " Meridional advection" ; ctrd(jpmld_yad,2) = "_yad"
991	ctrd(jpmld_zad,1) = " Vertical advection" ; ctrd(jpmld_zad,2) = "_zad"
992	ctrd(jpmld_ldf,1) = " Lateral diffusion" ; ctrd(jpmld_ldf,2) = "_ldf"
993	ctrd(jpmld_for,1) = " Forcing" ; ctrd(jpmld_for,2) = "_for"
994	ctrd(jpmld_zdf,1) = " Vertical diff. (Kz)" ; ctrd(jpmld_zdf,2) = "_zdf"
995	ctrd(jpmld_bbc,1) = " Geothermal flux" ; ctrd(jpmld_bbc,2) = "_bbc"
996	ctrd(jpmld_bbl,1) = " Adv/diff. Bottom boundary layer" ; ctrd(jpmld_bbl,2) = "_bbl"
997	ctrd(jpmld_dmp,1) = " Tracer damping" ; ctrd(jpmld_dmp,2) = "_dmp"
998	ctrd(jpmld_npc,1) = " Non penetrative convec. adjust." ; ctrd(jpmld_npc,2) = "_npc"
999	ctrd(jpmld_atf,1) = " Asselin time filter" ; ctrd(jpmld_atf,2) = "_atf"
1000
1001	!-- Create a NetCDF file and enter the define mode
1002	CALL dia_nam( clhstnam, nn_trd, 'trends' )
1003	IF(lwp) WRITE(numout,*) ' Name of NETCDF file ', clhstnam
1004	CALL histbeg( clhstnam, jpi, glamt, jpj, gphit, &
1005	& 1, jpi, 1, jpj, nit000-1, zjulian, rdt, nh_t, nidtrd, domain_id=nidom, snc4chunks=snc4set )
1006
1007	!-- Define the ML depth variable
1008	CALL histdef(nidtrd, "mxl_depth", clmxl//" Mixed Layer Depth" , "m", &
1009	jpi, jpj, nh_t, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
1010
1011	!-- Define physical units
1012	IF ( rn_ucf == 1. ) THEN ; cltu = "degC/s" ; clsu = "p.s.u./s"
1013	ELSEIF ( rn_ucf == 3600.*24.) THEN ; cltu = "degC/day" ; clsu = "p.s.u./day"
1014	ELSE ; cltu = "unknown?" ; clsu = "unknown?"
1015	END IF
1016
1017
1018	!-- Define miscellaneous T and S mixed-layer variables
1019
1020	IF( jpltrd /= jpmld_atf ) CALL ctl_stop( 'Error : jpltrd /= jpmld_atf' ) ! see below
1021
1022	!................................. ( ML temperature ) ...................................
1023
1024	CALL histdef(nidtrd, "tml" , clmxl//" T Mixed Layer Temperature" , "C", &
1025	jpi, jpj, nh_t, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
1026	CALL histdef(nidtrd, "tml_tot", clmxl//" T Total trend" , cltu, &
1027	jpi, jpj, nh_t, 1 , 1, 1 , -99 , 32, clop, zout, zout )
1028	CALL histdef(nidtrd, "tml_res", clmxl//" T dh/dt Entrainment (Resid.)" , cltu, &
1029	jpi, jpj, nh_t, 1 , 1, 1 , -99 , 32, clop, zout, zout )
1030
1031	DO jl = 1, jpltrd - 1 ! <== only true if jpltrd == jpmld_atf
1032	CALL histdef(nidtrd, trim("tml"//ctrd(jl,2)), clmxl//" T"//ctrd(jl,1), cltu, &
1033	jpi, jpj, nh_t, 1 , 1, 1 , -99 , 32, clop, zsto, zout ) ! IOIPSL: time mean
1034	END DO ! if zsto=rdt above
1035
1036	CALL histdef(nidtrd, trim("tml"//ctrd(jpmld_atf,2)), clmxl//" T"//ctrd(jpmld_atf,1), cltu, &
1037	jpi, jpj, nh_t, 1 , 1, 1 , -99 , 32, clop, zout, zout ) ! IOIPSL: NO time mean
1038
1039	!.................................. ( ML salinity ) .....................................
1040
1041	CALL histdef(nidtrd, "sml" , clmxl//" S Mixed Layer Salinity" , "p.s.u.", &
1042	jpi, jpj, nh_t, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
1043	CALL histdef(nidtrd, "sml_tot", clmxl//" S Total trend" , clsu, &
1044	jpi, jpj, nh_t, 1 , 1, 1 , -99 , 32, clop, zout, zout )
1045	CALL histdef(nidtrd, "sml_res", clmxl//" S dh/dt Entrainment (Resid.)" , clsu, &
1046	jpi, jpj, nh_t, 1 , 1, 1 , -99 , 32, clop, zout, zout )
1047
1048	DO jl = 1, jpltrd - 1 ! <== only true if jpltrd == jpmld_atf
1049	CALL histdef(nidtrd, trim("sml"//ctrd(jl,2)), clmxl//" S"//ctrd(jl,1), clsu, &
1050	jpi, jpj, nh_t, 1 , 1, 1 , -99 , 32, clop, zsto, zout ) ! IOIPSL: time mean
1051	END DO ! if zsto=rdt above
1052
1053	CALL histdef(nidtrd, trim("sml"//ctrd(jpmld_atf,2)), clmxl//" S"//ctrd(jpmld_atf,1), clsu, &
1054	jpi, jpj, nh_t, 1 , 1, 1 , -99 , 32, clop, zout, zout ) ! IOIPSL: NO time mean
1055
1056	!-- Leave IOIPSL/NetCDF define mode
1057	CALL histend( nidtrd, snc4set )
1058
1059	#endif /* key_dimgout */
1060	END SUBROUTINE trd_mld_init
1061
1062	#else
1063	!!----------------------------------------------------------------------
1064	!! Default option : Empty module
1065	!!----------------------------------------------------------------------
1066	CONTAINS
1067	SUBROUTINE trd_mld( kt ) ! Empty routine
1068	INTEGER, INTENT( in) :: kt
1069	WRITE(,) 'trd_mld: You should not have seen this print! error?', kt
1070	END SUBROUTINE trd_mld
1071	SUBROUTINE trd_mld_zint( pttrdmld, pstrdmld, ktrd, ctype )
1072	REAL, DIMENSION(:,:,:), INTENT( in ) :: &
1073	pttrdmld, pstrdmld ! Temperature and Salinity trends
1074	INTEGER, INTENT( in ) :: ktrd ! ocean trend index
1075	CHARACTER(len=2), INTENT( in ) :: &
1076	ctype ! surface/bottom (2D arrays) or
1077	! ! interior (3D arrays) physics
1078	WRITE(,) 'trd_mld_zint: You should not have seen this print! error?', pttrdmld(1,1,1)
1079	WRITE(,) ' " " : You should not have seen this print! error?', pstrdmld(1,1,1)
1080	WRITE(,) ' " " : You should not have seen this print! error?', ctype
1081	WRITE(,) ' " " : You should not have seen this print! error?', ktrd
1082	END SUBROUTINE trd_mld_zint
1083	SUBROUTINE trd_mld_init ! Empty routine
1084	WRITE(,) 'trd_mld_init: You should not have seen this print! error?'
1085	END SUBROUTINE trd_mld_init
1086	#endif
1087
1088	!!======================================================================
1089	END MODULE trdmld

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

Download in other formats: