Context Navigation

pbl_parameters.F @ 652

Last change on this file since 652 was 648, checked in by emillour, 13 years ago

Mars GCM:

some syntax corrections in thermcall_main_mars, vdif_cd, pbl_parameters which cause problems when compiling with some strict compilers (g95, gfortran)
added an initialisation of 'varian' in initracer for cases when using conrath dust; because that value can be is used elsewhere (e.g. surfacearea)

JYC+EM

File size: 17.1 KB

Rev	Line
[529]	1	SUBROUTINE pbl_parameters(ngrid,nlay,ps,pplay,pz0,
[566]	2	& pg,zzlay,zzlev,pu,pv,wstar_in,hfmax,zmax,pts,ph,z_out,n_out,
[605]	3	& T_out,u_out,ustar,tstar,L_mo,vhf,vvv)
[496]	4	IMPLICIT NONE
	5	!=======================================================================
	6	!
	7	! Anlysis of the PBL from input temperature, wind field and thermals outputs.
	8	!
	9	! -------
	10	!
	11	! Author: Arnaud Colaitis 09/01/12
	12	! -------
	13	!
	14	! Arguments:
	15	! ----------
	16	!
	17	! inputs:
	18	! ------
	19	! ngrid size of the horizontal grid
	20	! nlay size of the vertical grid
	21	! pz0(ngrid) surface roughness length
	22	! pg gravity (m s -2)
[529]	23	! zzlay(ngrid,nlay) height of mid-layers
	24	! zzlev(ngrid,nlay+1) height of layers interfaces
[496]	25	! pu(ngrid,nlay) u component of the wind
	26	! pv(ngrid,nlay) v component of the wind
[529]	27	! wstar_in(ngrid) free convection velocity in thermals
[496]	28	! hfmax(ngrid) maximum vertical turbulent heat flux in thermals
	29	! zmax(ngrid) height reached by the thermals (pbl height)
	30	! pts(ngrid) surface temperature
	31	! ph(ngrid,nlay) potential temperature T*(p/ps)^kappa
[566]	32	! z_out(n_out) heights of interpolation
	33	! n_out number of points for interpolation
[496]	34	!
	35	! outputs:
	36	! ------
	37	!
[566]	38	! Teta_out(ngrid,n_out) interpolated teta
	39	! u_out(ngrid,n_out) interpolated u
[496]	40	! ustar(ngrid) friction velocity
	41	! tstar(ngrid) friction temperature
	42	! wstar(ngrid) free convection velocity
	43	! L_mo(ngrid) monin_obukhov length
	44	!
	45	!
	46	!=======================================================================
	47	!
	48	!-----------------------------------------------------------------------
	49	! Declarations:
	50	! -------------
	51
	52	#include "comcstfi.h"
	53	#include "callkeys.h"
	54
	55	! Arguments:
	56	! ----------
	57
[566]	58	INTEGER, INTENT(IN) :: ngrid,nlay,n_out
[529]	59	REAL, INTENT(IN) :: pz0(ngrid),ps(ngrid),pplay(ngrid,nlay)
	60	REAL, INTENT(IN) :: pg,zzlay(ngrid,nlay),zzlev(ngrid,nlay)
[496]	61	REAL, INTENT(IN) :: pu(ngrid,nlay),pv(ngrid,nlay)
[529]	62	REAL, INTENT(IN) :: wstar_in(ngrid),hfmax(ngrid),zmax(ngrid)
[496]	63	REAL, INTENT(IN) :: pts(ngrid),ph(ngrid,nlay)
[566]	64	REAL, INTENT(IN) :: z_out(n_out)
[529]	65
	66	! Outputs:
	67	! --------
	68
[605]	69	REAL, INTENT(OUT) :: T_out(ngrid,n_out),u_out(ngrid,n_out)
	70	REAL Teta_out(ngrid,n_out)
[529]	71	REAL, INTENT(OUT) :: ustar(ngrid), tstar(ngrid)
	72	REAL wstar(ngrid)
[496]	73	REAL, INTENT(OUT) :: L_mo(ngrid)
	74
	75	! Local:
	76	! ------
	77
[566]	78	INTEGER ig,k,n
[496]	79	REAL karman,nu
	80	DATA karman,nu/.41,0.001/
	81	SAVE karman,nu
	82
	83	! Local(2):
	84	! ---------
	85
	86	REAL zout
	87	REAL rib(ngrid) ! Bulk Richardson number
	88	REAL fm(ngrid) ! stability function for momentum
	89	REAL fh(ngrid) ! stability function for heat
	90	REAL z1z0,z1z0t ! ratios z1/z0 and z1/z0T
	91	! phim = 1+betamzeta or (1-bmzeta)**am
	92	! phih = alphah + betahzeta or alphah(1.-bhzeta)**ah
	93	REAL betam, betah, alphah, bm, bh, lambda
	94	! ah and am are assumed to be -0.25 and -0.5 respectively
	95	REAL cdn(ngrid),chn(ngrid) ! neutral momentum and heat drag coefficient
	96	REAL pz0t ! initial thermal roughness length. (local)
	97	REAL ric ! critical richardson number
	98	REAL reynolds(ngrid) ! reynolds number for UBL
	99	REAL prandtl(ngrid) ! prandtl number for UBL
	100	REAL pz0tcomp(ngrid) ! computed z0t
	101	REAL ite
	102	REAL residual,zcd0,z1
	103	REAL pcdv(ngrid),pcdh(ngrid)
	104	REAL zu2(ngrid) ! Large-scale wind at first layer
	105	REAL pbl_teta(ngrid) ! mixed-layer potential temperature
	106	INTEGER pbl_height_index(ngrid) ! index of nearest vertical grid point for zmax
	107	REAL dteta(ngrid,nlay),x(ngrid) ! potential temperature gradient and z/zi
	108	REAL dvhf(ngrid),dvvv(ngrid) ! dimensionless vertical heat flux and
	109	! dimensionless vertical velocity variance
	110	REAL vhf(ngrid),vvv(ngrid) ! vertical heat flux and vertical velocity variance
	111	INTEGER ii(1)
	112	! temporary
	113	INTEGER dimout
	114
	115	!------------------------------------------------------------------------
	116	!------------------------------------------------------------------------
	117	! PART I : RICHARDSON/REYNOLDS/THERMAL_ROUGHNESS/STABILITY_COEFFICIENTS
	118	!------------------------------------------------------------------------
	119	!------------------------------------------------------------------------
	120
[566]	121	DO n=1,n_out
	122
[496]	123	c Initialisation :
	124
	125	L_mo(:)=0.
	126	ustar(:)=0.
	127	tstar(:)=0.
[566]	128	zout=z_out(n)
[496]	129	reynolds(:)=0.
	130	pz0t = 0.
	131	pz0tcomp(:) = 0.1*pz0(:)
	132	rib(:)=0.
	133	pcdv(:)=0.
	134	pcdh(:)=0.
	135
	136	! this formulation assumes alphah=1., implying betah=betam
	137	! We use Dyer et al. parameters, as they cover a broad range of Richardson numbers :
	138
	139	bm=16. !UBL
	140	bh=16. !UBL
	141	alphah=1.
	142	betam=5. !SBL
	143	betah=5. !SBL
	144	lambda=(sqrt(bh/bm))/alphah
	145	ric=betah/(betam**2)
	146	DO ig=1,ngrid
	147	ite=0.
	148	residual=abs(pz0tcomp(ig)-pz0t)
[636]	149	! zu2(ig)=MAX(pu(ig,1)pu(ig,1)+pv(ig,1)pv(ig,1)
	150	! & ,(0.3wstar_in(ig))*2)
[648]	151	zu2(ig)=pu(ig,1)pu(ig,1)+pv(ig,1)pv(ig,1)
[636]	152	& + (log(1.+0.7wstar_in(ig) + 2.3wstar_in(ig)2))2
[496]	153
	154	DO WHILE((residual .gt. 0.01*pz0(ig)) .and. (ite .lt. 10.))
	155
	156	pz0t=pz0tcomp(ig)
	157	IF (zu2(ig) .ne. 0.) THEN
	158	! Richardson number formulation proposed by D.E. England et al. (1995)
[605]	159	rib(ig) = (pg/pts(ig))
[636]	160	& sqrt(zzlay(ig,1)pz0(ig))
	161	& (((log(zzlay(ig,1)/pz0(ig)))*2)/(log(zzlay(ig,1)/pz0t)))
[496]	162	& *(ph(ig,1)-pts(ig))/zu2(ig)
	163	ELSE
	164	print*,'warning, infinite Richardson at surface'
	165	print*,pu(ig,1),pv(ig,1)
	166	rib(ig) = ric
	167	ENDIF
	168
[636]	169	z1z0=zzlay(ig,1)/pz0(ig)
	170	z1z0t=zzlay(ig,1)/pz0t
[496]	171
	172	cdn(ig)=karman/log(z1z0)
	173	cdn(ig)=cdn(ig)*cdn(ig)
	174	chn(ig)=cdn(ig)*log(z1z0)/log(z1z0t)
	175
	176	! STABLE BOUNDARY LAYER :
	177	IF (rib(ig) .gt. 0.) THEN
	178	! From D.E. England et al. (95)
	179	prandtl(ig)=1.
	180	if(rib(ig) .lt. ric) then
	181	! Assuming alphah=1. and bh=bm for stable conditions :
	182	fm(ig)=((ric-rib(ig))/ric)**2
	183	fh(ig)=((ric-rib(ig))/ric)**2
	184	else
	185	! For Ri>Ric, we consider Ri->Infinity => no turbulent mixing at surface
	186	fm(ig)=0.
	187	fh(ig)=0.
	188	endif
	189	! UNSTABLE BOUNDARY LAYER :
	190	ELSE
	191	! From D.E. England et al. (95)
	192	fm(ig)=sqrt(1.-lambdabmrib(ig))
	193	fh(ig)=(1./alphah)((1.-lambdabhrib(ig))0.5)
	194	& (1.-lambdabmrib(ig))**0.25
	195	prandtl(ig)=alphah((1.-lambdabmrib(ig))*0.25)/
	196	& ((1.-lambdabhrib(ig))**0.5)
	197	ENDIF
	198
	199	reynolds(ig)=karman*sqrt(fm(ig))
	200	& *sqrt(zu2(ig))
	201	& pz0(ig)/(log(z1z0)nu)
	202	pz0tcomp(ig)=pz0(ig)exp(-karman7.3*
	203	& (reynolds(ig)*0.25)(prandtl(ig)**0.5))
	204	residual = abs(pz0t-pz0tcomp(ig))
	205	ite = ite+1
[605]	206
[496]	207	ENDDO ! of while
	208	pz0t=0.
	209
	210	! Drag computation:
	211
	212	pcdv(ig)=cdn(ig)*fm(ig)
	213	pcdh(ig)=chn(ig)*fh(ig)
	214
	215	ENDDO ! of ngrid
	216
	217	!------------------------------------------------------------------------
	218	!------------------------------------------------------------------------
	219	! PART II : USTAR/TSTAR/U_OUT/TETA_OUT COMPUTATIONS
	220	!------------------------------------------------------------------------
	221	!------------------------------------------------------------------------
	222
	223
	224	! u* theta* computation
[636]	225	! and Monin Obukhov length:
[496]	226
	227	DO ig=1,ngrid
	228	IF (rib(ig) .ge. ric) THEN
	229	ustar(ig)=0.
	230	tstar(ig)=0.
[636]	231	L_mo(ig)=0.
[496]	232	ELSE
	233	ustar(ig)=sqrt(pcdv(ig))
	234	& sqrt(pu(ig,1)pu(ig,1)+pv(ig,1)*pv(ig,1))
	235	tstar(ig)=-pcdh(ig)*(pts(ig)-ph(ig,1))
	236	& /sqrt(pcdv(ig))
[636]	237	L_mo(ig)=pts(ig)(ustar(ig)2)/(pgkarman*tstar(ig)) !as defined here, L is positive for SBL, negative for UBL
[496]	238	ENDIF
	239	ENDDO
	240
	241	! Monin Obukhov length:
	242
[636]	243	! DO ig=1,ngrid
	244	! IF (rib(ig) .ge. ric) THEN
	245	! L_mo(ig)=0.
	246	! ELSE
	247	! L_mo(ig)=pts(ig)(ustar(ig)2)/(pgkarman*tstar(ig)) !as defined here, L is positive for SBL, negative for UBL
	248	! ENDIF
	249	! ENDDO
[496]	250
	251	! Interpolation:
	252
	253	DO ig=1,ngrid
	254	IF(zout .lt. pz0tcomp(ig)) THEN
[566]	255	u_out(ig,n)=0.
	256	Teta_out(ig,n)=pts(ig)
[496]	257	ELSEIF (L_mo(ig) .gt. 0.) THEN
[636]	258	! u_out(ig,n)=(ustar(ig)/karman)*log(zout/pz0(ig)) +
	259	! & 5.(ustar(ig)/(karmanL_mo(ig)))*(zout-pz0(ig))
	260	! Teta_out(ig,n)=pts(ig)+(tstar(ig)/(prandtl(ig)*karman))
	261	! & *log(zout/pz0tcomp(ig)) +
	262	! & 5.(tstar(ig)/(prandtl(ig)karman*L_mo(ig)))
	263	! & *(zout-pz0tcomp(ig))
	264	IF ((zout/L_mo(ig).lt.ric).and.(pz0(ig).lt.ric)) THEN
	265	! & ((zout/L_mo(ig).gt.ric).and.(pz0(ig).gt.ric)) ) THEN
	266
	267	u_out(ig,n)=(ustar(ig)/karman)*(
	268	& log((ric-pz0(ig)/L_mo(ig))/(ric-zout/L_mo(ig)))
	269	& +log(zout/pz0(ig))
	270	& )
	271	ELSEIF ((zout/L_mo(ig).gt.ric).and.(pz0(ig).gt.ric)) THEN
	272
	273	u_out(ig,n)=(ustar(ig)/karman)*(
	274	& log((zout-ricL_mo(ig))/(pz0(ig)-ricL_mo(ig)))
	275	& +log(pz0(ig)/zout)
	276	& )
	277
	278	ELSEIF ((zout/L_mo(ig).gt.ric).and.(pz0(ig).lt.ric)) THEN
	279	!integral of the stability function does not converge
	280
	281
	282	u_out(ig,n)=sqrt(pu(ig,1)2 + pv(ig,1)2)
	283
	284
	285	ENDIF
	286	IF((zout/L_mo(ig).lt.ric).and.(pz0tcomp(ig).lt.ric)) THEN
	287	! & ((zout/L_mo(ig).gt.ric).and.(pz0tcomp(ig).gt.ric)) ) THEN
	288
	289	Teta_out(ig,n)=pts(ig)+(tstar(ig)/(prandtl(ig)karman))(
	290	& log((ric-pz0tcomp(ig)/L_mo(ig))/(ric-zout/L_mo(ig)))
	291	& +log(zout/pz0tcomp(ig))
	292	& )
	293	ELSEIF ((zout/L_mo(ig).gt.ric).and.(pz0tcomp(ig).gt.ric)) THEN
	294
	295	Teta_out(ig,n)=pts(ig)+(tstar(ig)/(prandtl(ig)karman))(
	296	& log((zout-ricL_mo(ig))/(pz0tcomp(ig)-ricL_mo(ig)))
	297	& +log(pz0tcomp(ig)/zout)
	298	& )
	299
	300	ELSEIF ((zout/L_mo(ig).gt.ric).and.(pz0tcomp(ig).lt.ric)) THEN
	301	!integral of the stability function does not converge
	302
	303	Teta_out(ig,n)=ph(ig,1)
	304
	305	ENDIF
	306
[496]	307	ELSEIF (L_mo(ig) .lt. 0.) THEN
	308
	309	IF(L_mo(ig) .gt. -1000.) THEN
	310
[566]	311	u_out(ig,n)=(ustar(ig)/karman)*(
[496]	312	& 2.atan((1.-16.zout/L_mo(ig))**0.25)
	313	& -2.atan((1.-16.pz0(ig)/L_mo(ig))**0.25)
	314	& -2.log(1.+(1.-16.zout/L_mo(ig))**0.25)
	315	& +2.log(1.+(1.-16.pz0(ig)/L_mo(ig))**0.25)
	316	& - log(1.+sqrt(1.-16.*zout/L_mo(ig)))
	317	& + log(1.+sqrt(1.-16.*pz0(ig)/L_mo(ig)))
	318	& + log(zout/pz0(ig))
	319	& )
	320
[605]	321	Teta_out(ig,n)=MAX(pts(ig)+(tstar(ig)/(prandtl(ig)karman))(
[496]	322	& 2.log(1.+sqrt(1.-16.pz0tcomp(ig)/L_mo(ig)))
	323	& -2.log(1.+sqrt(1.-16.zout/L_mo(ig)))
	324	& + log(zout/pz0tcomp(ig))
[605]	325	& ),ph(ig,1))
[496]	326
	327	ELSE
	328
	329	! We have to treat the case where L is very large and negative,
	330	! corresponding to a very nearly stable atmosphere (but not quite !)
	331	! i.e. teta* <0 and almost zero. We choose the cutoff
	332	! corresponding to L_mo < -1000 and use a 3rd order taylor expansion. The difference
	333	! between the two expression for z/L = -1/1000 is -1.54324*10^-9
	334	! (we do that to avoid using r*4 precision, otherwise, we get -inf values)
	335
[566]	336	u_out(ig,n)=(ustar(ig)/karman)*(
[496]	337	& (4./L_mo(ig))*(zout-pz0(ig))
	338	& + (20./(L_mo(ig))*2)(zout2-pz0(ig)2)
	339	& + (160./(L_mo(ig))*3)(zout3-pz0(ig)3)
	340	& + log(zout/pz0(ig))
	341	& )
	342
[605]	343	Teta_out(ig,n)=MAX(pts(ig)+(tstar(ig)/(prandtl(ig)karman))(
[496]	344	& (8./L_mo(ig))*(zout-pz0tcomp(ig))
	345	& + (48./(L_mo(ig))*2)(zout2-pz0tcomp(ig)2)
	346	& + (1280./(3.(L_mo(ig))3))(zout3-pz0tcomp(ig)3)
	347	& + log(zout/pz0tcomp(ig))
[605]	348	& ),ph(ig,1))
[496]	349
	350	ENDIF
	351	ELSE
[636]	352	u_out(ig,n)=sqrt(pu(ig,1)2 + pv(ig,1)2)
	353	Teta_out(ig,n)=ph(ig,1)
[496]	354	ENDIF
	355	IF(zout .lt. pz0(ig)) THEN
[566]	356	u_out(ig,n)=0.
[496]	357	ENDIF
[636]	358
	359	! Final check
	360	IF (L_mo(ig) .gt. 0) THEN
	361	IF (Teta_out(ig,n) .gt. ph(ig,1)) THEN
	362	Teta_out(ig,n)=ph(ig,1)
	363	ELSEIF (Teta_out(ig,n) .lt. pts(ig)) THEN
	364	Teta_out(ig,n)=pts(ig)
	365	ENDIF
	366	ELSEIF (L_mo(ig) .lt. 0) THEN
	367	IF (Teta_out(ig,n) .lt. ph(ig,1)) THEN
	368	Teta_out(ig,n)=ph(ig,1)
	369	ELSEIF (Teta_out(ig,n) .gt. pts(ig)) THEN
	370	Teta_out(ig,n)=pts(ig)
	371	ENDIF
	372	ENDIF
	373
	374	IF (u_out(ig,n) .gt. sqrt(pu(ig,1)2 + pv(ig,1)2)) THEN
	375	u_out(ig,n)=sqrt(pu(ig,1)2 + pv(ig,1)2)
	376	ENDIF
	377
[496]	378	ENDDO
	379
	380	! when using convective adjustment without thermals, a vertical potential temperature
	381	! profile is assumed up to the thermal roughness length. Hence, theoretically, theta
	382	! interpolated at any height in the surface layer is theta at the first level.
	383
	384	IF ((.not.calltherm).and.(calladj)) THEN
[605]	385	Teta_out(:,n)=ph(:,1)
	386	u_out(:,n)=(sqrt(cdn(:))sqrt(pu(:,1)pu(:,1)+pv(:,1)*pv(:,1))
	387	& /karman)*log(zout/pz0(:))
[496]	388	ENDIF
[566]	389	T_out(:,n) = Teta_out(:,n)*(exp(
[529]	390	& (zout/zzlay(:,1))*(log(pplay(:,1)/ps))
	391	& )
	392	& )**rcp
	393
[566]	394	ENDDO !of n=1,n_out
[529]	395
[496]	396	!------------------------------------------------------------------------
	397	!------------------------------------------------------------------------
	398	! PART III : WSTAR COMPUTATION
	399	!------------------------------------------------------------------------
	400	!------------------------------------------------------------------------
	401
	402	! Detection of the mixed-layer potential temperature
	403	! ------------
	404
	405	! Nearest index for the pbl height
	406
[566]	407	IF (calltherm) THEN
	408
[496]	409	pbl_height_index(:)=1
	410
	411	DO k=1,nlay-1
	412	DO ig=1, ngrid
[529]	413	IF (abs(zmax(ig)-zzlay(ig,k)) .lt.
	414	& abs(zmax(ig)-zzlay(ig,pbl_height_index(ig)))) THEN
[496]	415	pbl_height_index(ig)=k
	416	ENDIF
	417	ENDDO
	418	ENDDO
	419
	420	! Potential temperature gradient
	421
	422	dteta(:,nlay)=0.
	423	DO k=1,nlay-1
	424	DO ig=1, ngrid
[529]	425	dteta(ig,k) = (ph(ig,k+1)-ph(ig,k))/(zzlay(ig,k+1)-zzlay(ig,k))
[496]	426	ENDDO
	427	ENDDO
	428
	429	! Computation of the pbl mixed layer temperature
	430
	431	DO ig=1, ngrid
	432	ii=MINLOC(abs(dteta(ig,1:pbl_height_index(ig))))
	433	pbl_teta(ig) = ph(ig,ii(1))
	434	ENDDO
	435
[566]	436	! Recompute wstar
[496]	437	! We follow Spiga et. al 2010 (QJRMS)
	438	! ------------
	439
	440	DO ig=1, ngrid
	441	IF (zmax(ig) .gt. 0.) THEN
	442	wstar(ig)=(pgzmax(ig)hfmax(ig)/pbl_teta(ig))**(1./3.)
	443	ELSE
	444	wstar(ig)=0.
	445	ENDIF
	446	ENDDO
	447
	448	!------------------------------------------------------------------------
	449	!------------------------------------------------------------------------
	450	! PART IV : VERTICAL_VELOCITY_VARIANCE/VERTICAL_TURBULENT_FLUX PROFILES
	451	!------------------------------------------------------------------------
	452	!------------------------------------------------------------------------
	453
	454	! We follow Spiga et. al 2010 (QJRMS)
	455	! ------------
	456
	457	DO ig=1, ngrid
	458	IF (zmax(ig) .gt. 0.) THEN
	459	x(ig) = zout/zmax(ig)
	460	ELSE
	461	x(ig) = 999.
	462	ENDIF
	463	ENDDO
	464
	465	DO ig=1, ngrid
	466	! dimensionless vertical heat flux
	467	IF (x(ig) .le. 0.3) THEN
	468	dvhf(ig) = ((-3.85/log(x(ig)))+0.07*log(x(ig)))
	469	& exp(-4.61x(ig))
	470	ELSEIF (x(ig) .le. 1.) THEN
	471	dvhf(ig) = -1.52*x(ig) + 1.24
	472	ELSE
	473	dvhf(ig) = 0.
	474	ENDIF
	475	! dimensionless vertical velocity variance
	476	IF (x(ig) .le. 1.) THEN
	477	dvvv(ig) = 2.05(x(ig)(2./3.))(1.-0.64x(ig))*2
	478	ELSE
	479	dvvv(ig) = 0.
	480	ENDIF
	481	ENDDO
	482
	483	vhf(:) = dvhf(:)*hfmax(:)
	484	vvv(:) = dvvv(:)(wstar(:))*2
	485
[566]	486	ENDIF ! of if calltherm
[496]	487
[636]	488	call WRITEDIAGFI(ngrid,'rib_pbl',
	489	& 'Richardson in pbl parameter','m/s',2,rib)
	490	call WRITEDIAGFI(ngrid,'cdn_pbl',
	491	& 'neutral momentum coef','m/s',2,cdn)
	492	call WRITEDIAGFI(ngrid,'fm_pbl',
	493	& 'momentum stab function','m/s',2,fm)
	494	call WRITEDIAGFI(ngrid,'uv',
	495	& 'wind norm first layer','m/s',2,sqrt(zu2))
	496	call WRITEDIAGFI(ngrid,'uvtrue',
	497	& 'wind norm first layer','m/s',2,sqrt(pu(:,1)2+pv(:,1)2))
	498	call WRITEDIAGFI(ngrid,'chn_pbl',
	499	& 'neutral momentum coef','m/s',2,chn)
	500	call WRITEDIAGFI(ngrid,'fh_pbl',
	501	& 'momentum stab function','m/s',2,fh)
	502	call WRITEDIAGFI(ngrid,'B_pbl',
	503	& 'flottabilité','m/',2,(ph(:,1)-pts(:))/pts(:))
	504	call WRITEDIAGFI(ngrid,'zot_pbl',
	505	& 'flottabilité','ms',2,pz0tcomp)
	506	call WRITEDIAGFI(ngrid,'zz1','flottabilité','m',2,zzlay(:,1))
	507
[496]	508	RETURN
	509	END

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

Download in other formats:

Original Format