Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

trdmxl.F90 @ 11738

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

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

New URL for NEMO forge! http://forge.nemo-ocean.eu

Context Navigation

source: branches/UKMO/dev_r5518_GO6_under_ice_relax_dr_hook/NEMOGCM/NEMO/OPA_SRC/TRD/trdmxl.F90 @ 11738

Download in other formats: