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

Context Navigation

source: trunk/NEMOGCM/NEMO/OPA_SRC/TRD/trdmld.F90 @ 3606

Last change on this file since 3606 was 3294, checked in by rblod, 12 years ago
Merge of 3.4beta into the trunk
Property svn:keywords set to `Id`
File size: 48.6 KB

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

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

Download in other formats: