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.

diahth.F90 in branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/DIA – NEMO

Context Navigation

source: branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/DIA/diahth.F90 @ 4416

Last change on this file since 4416 was 3211, checked in by spickles2, 12 years ago

Stephen Pickles, 11 Dec 2011

Commit to bring the rest of the DCSE NEMO development branch
in line with the latest development version. This includes
array index re-ordering of all OPA_SRC/.

Property svn:keywords set to Id

File size: 17.5 KB

Rev	Line
[3]	1	MODULE diahth
	2	!!======================================================================
	3	!! * MODULE diahth *
	4	!! Ocean diagnostics: thermocline and 20 degree depth
	5	!!======================================================================
[1485]	6	!! History : OPA ! 1994-09 (J.-P. Boulanger) Original code
	7	!! ! 1996-11 (E. Guilyardi) OPA8
	8	!! ! 1997-08 (G. Madec) optimization
	9	!! ! 1999-07 (E. Guilyardi) hd28 + heat content
	10	!! 8.5 ! 2002-06 (G. Madec) F90: Free form and module
[1577]	11	!! NEMO 3.2 ! 2009-07 (S. Masson) hc300 bugfix + cleaning + add new diag
[1485]	12	!!----------------------------------------------------------------------
[3]	13	#if defined key_diahth \|\| defined key_esopa
	14	!!----------------------------------------------------------------------
	15	!! 'key_diahth' : thermocline depth diag.
	16	!!----------------------------------------------------------------------
[1577]	17	!! dia_hth : Compute varius diagnostics associated with the mixed layer
[3]	18	!!----------------------------------------------------------------------
	19	USE oce ! ocean dynamics and tracers
	20	USE dom_oce ! ocean space and time domain
	21	USE phycst ! physical constants
	22	USE in_out_manager ! I/O manager
[2715]	23	USE lib_mpp ! MPP library
[2528]	24	USE iom ! I/O library
[3]	25
	26	IMPLICIT NONE
	27	PRIVATE
	28
[2715]	29	PUBLIC dia_hth ! routine called by step.F90
	30	PUBLIC dia_hth_alloc ! routine called by nemogcm.F90
[3]	31
[2715]	32	LOGICAL , PUBLIC, PARAMETER :: lk_diahth = .TRUE. !: thermocline-20d depths flag
[1577]	33	! note: following variables should move to local variables once iom_put is always used
[2715]	34	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hth !: depth of the max vertical temperature gradient [m]
	35	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hd20 !: depth of 20 C isotherm [m]
	36	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hd28 !: depth of 28 C isotherm [m]
	37	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: htc3 !: heat content of first 300 m [W]
[3]	38
[3211]	39	!! * Control permutation of array indices
	40	# include "oce_ftrans.h90"
	41	# include "dom_oce_ftrans.h90"
	42
[3]	43	!! * Substitutions
	44	# include "domzgr_substitute.h90"
	45	!!----------------------------------------------------------------------
[2715]	46	!! NEMO/OPA 4.0 , NEMO Consortium (2011)
[1152]	47	!! $Id$
[2528]	48	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
[3]	49	!!----------------------------------------------------------------------
	50	CONTAINS
	51
[2715]	52	FUNCTION dia_hth_alloc()
	53	!!---------------------------------------------------------------------
	54	INTEGER :: dia_hth_alloc
	55	!!---------------------------------------------------------------------
	56	!
	57	ALLOCATE(hth(jpi,jpj), hd20(jpi,jpj), hd28(jpi,jpj), htc3(jpi,jpj), STAT=dia_hth_alloc)
	58	!
	59	IF( lk_mpp ) CALL mpp_sum ( dia_hth_alloc )
	60	IF(dia_hth_alloc /= 0) CALL ctl_warn('dia_hth_alloc: failed to allocate arrays.')
	61	!
	62	END FUNCTION dia_hth_alloc
	63
	64
[3]	65	SUBROUTINE dia_hth( kt )
	66	!!---------------------------------------------------------------------
	67	!! * ROUTINE dia_hth *
	68	!!
[1577]	69	!! ** Purpose : Computes
	70	!! the mixing layer depth (turbocline): avt = 5.e-4
	71	!! the depth of strongest vertical temperature gradient
	72	!! the mixed layer depth with density criteria: rho = rho(10m or surf) + 0.03(or 0.01)
	73	!! the mixed layer depth with temperature criteria: abs( tn - tn(10m) ) = 0.2
	74	!! the top of the thermochine: tn = tn(10m) - ztem2
	75	!! the pycnocline depth with density criteria equivalent to a temperature variation
	76	!! rho = rho10m + (dr/dT)(T,S,10m)*(-0.2 degC)
	77	!! the barrier layer thickness
	78	!! the maximal verical inversion of temperature and its depth max( 0, max of tn - tn(10m) )
	79	!! the depth of the 20 degree isotherm (linear interpolation)
	80	!! the depth of the 28 degree isotherm (linear interpolation)
	81	!! the heat content of first 300 m
[3]	82	!!
	83	!! ** Method :
	84	!!-------------------------------------------------------------------
	85	INTEGER, INTENT( in ) :: kt ! ocean time-step index
[1485]	86	!!
	87	INTEGER :: ji, jj, jk ! dummy loop arguments
[2528]	88	INTEGER :: iid, ilevel ! temporary integers
[2715]	89	INTEGER, ALLOCATABLE, SAVE, DIMENSION(:,:) :: ik20, ik28 ! levels
[1577]	90	REAL(wp) :: zavt5 = 5.e-4_wp ! Kz criterion for the turbocline depth
	91	REAL(wp) :: zrho3 = 0.03_wp ! density criterion for mixed layer depth
	92	REAL(wp) :: zrho1 = 0.01_wp ! density criterion for mixed layer depth
	93	REAL(wp) :: ztem2 = 0.2_wp ! temperature criterion for mixed layer depth
	94	REAL(wp) :: zthick_0, zcoef ! temporary scalars
	95	REAL(wp) :: zztmp, zzdep ! temporary scalars inside do loop
	96	REAL(wp) :: zu, zv, zw, zut, zvt ! temporary workspace
[2715]	97	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: zabs2 ! MLD: abs( tn - tn(10m) ) = ztem2
	98	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: ztm2 ! Top of thermocline: tn = tn(10m) - ztem2
	99	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: zrho10_3 ! MLD: rho = rho10m + zrho3
	100	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: zpycn ! pycnocline: rho = rho10m + (dr/dT)(T,S,10m)*(-0.2 degC)
	101	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: ztinv ! max of temperature inversion
	102	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: zdepinv ! depth of temperature inversion
	103	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: zrho0_3 ! MLD rho = rho(surf) = 0.03
	104	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: zrho0_1 ! MLD rho = rho(surf) = 0.01
	105	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: zmaxdzT ! max of dT/dz
	106	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: zthick ! vertical integration thickness
	107	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: zdelr ! delta rho equivalent to deltaT = 0.2
[3]	108	!!----------------------------------------------------------------------
	109
	110	IF( kt == nit000 ) THEN
[2715]	111	! ! allocate dia_hth array
	112	IF( dia_hth_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'lim_sbc_init : unable to allocate standard arrays' )
	113
	114	IF(.not. ALLOCATED(ik20))THEN
	115	ALLOCATE(ik20(jpi,jpj), ik28(jpi,jpj), &
	116	& zabs2(jpi,jpj), &
	117	& ztm2(jpi,jpj), &
	118	& zrho10_3(jpi,jpj),&
	119	& zpycn(jpi,jpj), &
	120	& ztinv(jpi,jpj), &
	121	& zdepinv(jpi,jpj), &
	122	& zrho0_3(jpi,jpj), &
	123	& zrho0_1(jpi,jpj), &
	124	& zmaxdzT(jpi,jpj), &
	125	& zthick(jpi,jpj), &
	126	& zdelr(jpi,jpj), STAT=ji)
	127	IF( lk_mpp ) CALL mpp_sum(ji)
	128	IF( ji /= 0 ) CALL ctl_stop( 'STOP', 'dia_hth : unable to allocate standard ocean arrays' )
	129	END IF
	130
[3]	131	IF(lwp) WRITE(numout,*)
	132	IF(lwp) WRITE(numout,*) 'dia_hth : diagnostics of the thermocline depth'
	133	IF(lwp) WRITE(numout,*) '~~~~~~~ '
	134	IF(lwp) WRITE(numout,*)
	135	ENDIF
	136
[1577]	137	! initialization
[2528]	138	ztinv (:,:) = 0._wp
	139	zdepinv(:,:) = 0._wp
	140	zmaxdzT(:,:) = 0._wp
[3]	141	DO jj = 1, jpj
	142	DO ji = 1, jpi
[1577]	143	zztmp = bathy(ji,jj)
	144	hth (ji,jj) = zztmp
	145	zabs2 (ji,jj) = zztmp
	146	ztm2 (ji,jj) = zztmp
	147	zrho10_3(ji,jj) = zztmp
	148	zpycn (ji,jj) = zztmp
	149	END DO
[3]	150	END DO
[1577]	151	IF( nla10 > 1 ) THEN
	152	DO jj = 1, jpj
	153	DO ji = 1, jpi
	154	zztmp = bathy(ji,jj)
	155	zrho0_3(ji,jj) = zztmp
	156	zrho0_1(ji,jj) = zztmp
	157	END DO
	158	END DO
	159	ENDIF
	160
	161	! Preliminary computation
	162	! computation of zdelr = (dr/dT)(T,S,10m)*(-0.2 degC)
[2528]	163	DO jj = 1, jpj
	164	DO ji = 1, jpi
[1577]	165	IF( tmask(ji,jj,nla10) == 1. ) THEN
	166	zu = 1779.50 + 11.250tn(ji,jj,nla10) - 3.80sn(ji,jj,nla10) - 0.0745tn(ji,jj,nla10)tn(ji,jj,nla10) &
	167	& - 0.0100tn(ji,jj,nla10)sn(ji,jj,nla10)
	168	zv = 5891.00 + 38.000tn(ji,jj,nla10) + 3.00sn(ji,jj,nla10) - 0.3750tn(ji,jj,nla10)tn(ji,jj,nla10)
	169	zut = 11.25 - 0.149tn(ji,jj,nla10) - 0.01sn(ji,jj,nla10)
	170	zvt = 38.00 - 0.750*tn(ji,jj,nla10)
	171	zw = (zu + 0.698zv) (zu + 0.698*zv)
	172	zdelr(ji,jj) = ztem2 * (1000.(zutzv - zvt*zu)/zw)
	173	ELSE
[2528]	174	zdelr(ji,jj) = 0._wp
[1577]	175	ENDIF
	176	END DO
	177	END DO
[3]	178
[1577]	179	! ------------------------------------------------------------- !
	180	! thermocline depth: strongest vertical gradient of temperature !
	181	! turbocline depth (mixing layer depth): avt = zavt5 !
	182	! MLD: rho = rho(1) + zrho3 !
	183	! MLD: rho = rho(1) + zrho1 !
	184	! ------------------------------------------------------------- !
[3211]	185	#if defined key_z_first
	186	DO jj = 1, jpj
	187	DO ji = 1, jpi
	188	DO jk = jpkm1, 2, -1 ! loop from bottom to 2
	189	#else
[1577]	190	DO jk = jpkm1, 2, -1 ! loop from bottom to 2
	191	DO jj = 1, jpj
	192	DO ji = 1, jpi
[3211]	193	#endif
[2528]	194	!
[1577]	195	zzdep = fsdepw(ji,jj,jk)
	196	zztmp = ( tn(ji,jj,jk-1) - tn(ji,jj,jk) ) / zzdep * tmask(ji,jj,jk) ! vertical gradient of temperature (dT/dz)
	197	zzdep = zzdep * tmask(ji,jj,1)
	198
	199	IF( zztmp > zmaxdzT(ji,jj) ) THEN
	200	zmaxdzT(ji,jj) = zztmp ; hth (ji,jj) = zzdep ! max and depth of dT/dz
	201	ENDIF
	202
	203	IF( nla10 > 1 ) THEN
	204	zztmp = rhop(ji,jj,jk) - rhop(ji,jj,1) ! delta rho(1)
	205	IF( zztmp > zrho3 ) zrho0_3(ji,jj) = zzdep ! > 0.03
	206	IF( zztmp > zrho1 ) zrho0_1(ji,jj) = zzdep ! > 0.01
	207	ENDIF
	208
	209	END DO
	210	END DO
	211	END DO
[3]	212
[1577]	213	CALL iom_put( "mlddzt", hth ) ! depth of the thermocline
	214	IF( nla10 > 1 ) THEN
	215	CALL iom_put( "mldr0_3", zrho0_3 ) ! MLD delta rho(surf) = 0.03
	216	CALL iom_put( "mldr0_1", zrho0_1 ) ! MLD delta rho(surf) = 0.01
	217	ENDIF
	218
	219	! ------------------------------------------------------------- !
	220	! MLD: abs( tn - tn(10m) ) = ztem2 !
	221	! Top of thermocline: tn = tn(10m) - ztem2 !
	222	! MLD: rho = rho10m + zrho3 !
	223	! pycnocline: rho = rho10m + (dr/dT)(T,S,10m)*(-0.2 degC) !
	224	! temperature inversion: max( 0, max of tn - tn(10m) ) !
	225	! depth of temperature inversion !
	226	! ------------------------------------------------------------- !
[3211]	227	#if defined key_z_first
	228	DO jj = 1, jpj
	229	DO ji = 1, jpi
	230	DO jk = jpkm1, nlb10, -1 ! loop from bottom to nlb10
	231	#else
[1577]	232	DO jk = jpkm1, nlb10, -1 ! loop from bottom to nlb10
	233	DO jj = 1, jpj
	234	DO ji = 1, jpi
[3211]	235	#endif
[2528]	236	!
[1577]	237	zzdep = fsdepw(ji,jj,jk) * tmask(ji,jj,1)
[2528]	238	!
[1577]	239	zztmp = tn(ji,jj,nla10) - tn(ji,jj,jk) ! - delta T(10m)
	240	IF( ABS(zztmp) > ztem2 ) zabs2 (ji,jj) = zzdep ! abs > 0.2
	241	IF( zztmp > ztem2 ) ztm2 (ji,jj) = zzdep ! > 0.2
	242	zztmp = -zztmp ! delta T(10m)
	243	IF( zztmp > ztinv(ji,jj) ) THEN ! temperature inversion
	244	ztinv(ji,jj) = zztmp ; zdepinv (ji,jj) = zzdep ! max value and depth
	245	ENDIF
	246
	247	zztmp = rhop(ji,jj,jk) - rhop(ji,jj,nla10) ! delta rho(10m)
	248	IF( zztmp > zrho3 ) zrho10_3(ji,jj) = zzdep ! > 0.03
	249	IF( zztmp > zdelr(ji,jj) ) zpycn (ji,jj) = zzdep ! > equi. delta T(10m) - 0.2
[2528]	250	!
[1577]	251	END DO
[3]	252	END DO
	253	END DO
	254
[2561]	255	CALL iom_put( "mld_dt02", zabs2 ) ! MLD abs(delta t) - 0.2
[1577]	256	CALL iom_put( "topthdep", ztm2 ) ! T(10) - 0.2
	257	CALL iom_put( "mldr10_3", zrho10_3 ) ! MLD delta rho(10m) = 0.03
	258	CALL iom_put( "pycndep" , zpycn ) ! MLD delta rho equi. delta T(10m) = 0.2
	259	CALL iom_put( "BLT" , ztm2 - zpycn ) ! Barrier Layer Thickness
	260	CALL iom_put( "tinv" , ztinv ) ! max. temp. inv. (t10 ref)
	261	CALL iom_put( "depti" , zdepinv ) ! depth of max. temp. inv. (t10 ref)
	262
	263
	264	! ----------------------------------- !
	265	! search deepest level above 20C/28C !
	266	! ----------------------------------- !
	267	ik20(:,:) = 1
	268	ik28(:,:) = 1
[3211]	269	#if defined key_z_first
	270	DO jj = 1, jpj
	271	DO ji = 1, jpi
	272	DO jk = 1, jpkm1 ! beware temperature is not always decreasing with depth =>
	273	! ! loop from top to bottom
	274	#else
[1577]	275	DO jk = 1, jpkm1 ! beware temperature is not always decreasing with depth => loop from top to bottom
	276	DO jj = 1, jpj
	277	DO ji = 1, jpi
[3211]	278	#endif
[1577]	279	zztmp = tn(ji,jj,jk)
	280	IF( zztmp >= 20. ) ik20(ji,jj) = jk
	281	IF( zztmp >= 28. ) ik28(ji,jj) = jk
	282	END DO
	283	END DO
	284	END DO
	285
	286	! --------------------------- !
	287	! Depth of 20C/28C isotherm !
	288	! --------------------------- !
[3]	289	DO jj = 1, jpj
	290	DO ji = 1, jpi
[2528]	291	!
	292	zzdep = fsdepw(ji,jj,mbkt(ji,jj)+1) ! depth of the oean bottom
	293	!
[1577]	294	iid = ik20(ji,jj)
	295	IF( iid /= 1 ) THEN
[2528]	296	zztmp = fsdept(ji,jj,iid ) & ! linear interpolation
[1577]	297	& + ( fsdept(ji,jj,iid+1) - fsdept(ji,jj,iid) ) &
	298	& * ( 20.*tmask(ji,jj,iid+1) - tn(ji,jj,iid) ) &
	299	& / ( tn(ji,jj,iid+1) - tn(ji,jj,iid) + (1.-tmask(ji,jj,1)) )
[2528]	300	hd20(ji,jj) = MIN( zztmp , zzdep) * tmask(ji,jj,1) ! bound by the ocean depth
[1577]	301	ELSE
[2528]	302	hd20(ji,jj) = 0._wp
[1577]	303	ENDIF
[2528]	304	!
[1577]	305	iid = ik28(ji,jj)
	306	IF( iid /= 1 ) THEN
[2528]	307	zztmp = fsdept(ji,jj,iid ) & ! linear interpolation
[1577]	308	& + ( fsdept(ji,jj,iid+1) - fsdept(ji,jj,iid) ) &
	309	& * ( 28.*tmask(ji,jj,iid+1) - tn(ji,jj,iid) ) &
	310	& / ( tn(ji,jj,iid+1) - tn(ji,jj,iid) + (1.-tmask(ji,jj,1)) )
[2528]	311	hd28(ji,jj) = MIN( zztmp , zzdep ) * tmask(ji,jj,1) ! bound by the ocean depth
[1577]	312	ELSE
[2528]	313	hd28(ji,jj) = 0._wp
[1577]	314	ENDIF
	315
[3]	316	END DO
	317	END DO
[1577]	318	CALL iom_put( "20d", hd20 ) ! depth of the 20 isotherm
	319	CALL iom_put( "28d", hd28 ) ! depth of the 28 isotherm
[3]	320
[1485]	321	! ----------------------------- !
	322	! Heat content of first 300 m !
	323	! ----------------------------- !
[3]	324
[1484]	325	! find ilevel with (ilevel+1) the deepest W-level above 300m (we assume we can use e3t_0 to do this search...)
[2528]	326	ilevel = 0
	327	zthick_0 = 0._wp
[1486]	328	DO jk = 1, jpkm1
[1484]	329	zthick_0 = zthick_0 + e3t_0(jk)
	330	IF( zthick_0 < 300. ) ilevel = jk
	331	END DO
	332	! surface boundary condition
[2528]	333	IF( lk_vvl ) THEN ; zthick(:,:) = 0._wp ; htc3(:,:) = 0._wp
[1484]	334	ELSE ; zthick(:,:) = sshn(:,:) ; htc3(:,:) = tn(:,:,jk) * sshn(:,:) * tmask(:,:,jk)
	335	ENDIF
	336	! integration down to ilevel
[3]	337	DO jk = 1, ilevel
[1484]	338	zthick(:,:) = zthick(:,:) + fse3t(:,:,jk)
	339	htc3 (:,:) = htc3 (:,:) + fse3t(:,:,jk) * tn(:,:,jk) * tmask(:,:,jk)
[3]	340	END DO
[1484]	341	! deepest layer
	342	zthick(:,:) = 300. - zthick(:,:) ! remaining thickness to reach 300m
[1551]	343	DO jj = 1, jpj
	344	DO ji = 1, jpi
	345	htc3(ji,jj) = htc3(ji,jj) + tn(ji,jj,ilevel+1) * MIN( fse3t(ji,jj,ilevel+1), zthick(ji,jj) ) * tmask(ji,jj,ilevel+1)
	346	END DO
	347	END DO
[1484]	348	! from temperature to heat contain
	349	zcoef = rau0 * rcp
	350	htc3(:,:) = zcoef * htc3(:,:)
[1577]	351	CALL iom_put( "hc300", htc3 ) ! first 300m heat content
[2528]	352	!
[3]	353	END SUBROUTINE dia_hth
	354
	355	#else
	356	!!----------------------------------------------------------------------
	357	!! Default option : Empty module
	358	!!----------------------------------------------------------------------
[32]	359	LOGICAL , PUBLIC, PARAMETER :: lk_diahth = .FALSE. !: thermocline-20d depths flag
[3]	360	CONTAINS
	361	SUBROUTINE dia_hth( kt ) ! Empty routine
[32]	362	WRITE(,) 'dia_hth: You should not have seen this print! error?', kt
[3]	363	END SUBROUTINE dia_hth
	364	#endif
	365
	366	!!======================================================================
	367	END MODULE diahth

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

Download in other formats: