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

Context Navigation

source: branches/2012/dev_r3309_LOCEAN12_Ediag/NEMOGCM/NEMO/OPA_SRC/TRD/trdmld.F90 @ 3325

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

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

Download in other formats: