Context Navigation

thermcell_plume.F90 @ 2065

Last change on this file since 2065 was 2065, checked in by aboissinot, 6 years ago

In thermcell_plume, replace a useless test on zalpha by a test on zw2m to avoid a possible
division by zero and remove useless variable zbuoybis.

In physiq_mod, save variable f0.

File size: 23.5 KB

Rev	Line
[2060]	1	!
	2	!
	3	!
	4	SUBROUTINE thermcell_plume(itap,ngrid,klev,ptimestep,ztv, &
	5	& zthl,po,zl,rhobarz,zlev,pplev,pphi,zpspsk, &
	6	& alim_star,alim_star_tot,lalim,f0,detr_star, &
	7	& entr_star,f_star,ztva,ztla,zqla,zqta,zha, &
	8	& zw2,w_est,ztva_est,zqsatth,lmix,lmix_bis, &
	9	& lmin,lev_out,lunout1,igout)
	10
	11
	12	!==============================================================================
	13	! thermcell_plume: calcule les valeurs de qt, thetal et w dans l ascendance
	14	!==============================================================================
	15
	16	USE IOIPSL, ONLY : getin
	17	USE ioipsl_getin_p_mod, ONLY : getin_p
	18	USE print_control_mod, ONLY: prt_level
	19	USE watercommon_h, ONLY: RLvCp, RETV, watersat
	20	USE thermcell_mod
	21
	22	IMPLICIT NONE
	23
	24
	25	!==============================================================================
	26	! Declaration
	27	!==============================================================================
	28
	29	! inputs:
	30	! -------
	31
	32	INTEGER itap
	33	INTEGER ngrid
	34	INTEGER klev
	35	INTEGER lunout1
	36	INTEGER igout
	37	INTEGER lev_out ! niveau pour les print
	38
	39	REAL ptimestep ! time step
	40	REAL ztv(ngrid,klev) ! TRPV environment
	41	REAL zthl(ngrid,klev) ! TP environment
	42	REAL po(ngrid,klev) ! qt environment
	43	REAL zl(ngrid,klev) ! ql environment
	44	REAL rhobarz(ngrid,klev) ! levels density
	45	REAL zlev(ngrid,klev+1) ! levels altitude
	46	REAL pplev(ngrid,klev+1) ! levels pressure
	47	REAL pphi(ngrid,klev) ! geopotential
	48	REAL zpspsk(ngrid,klev) ! Exner function
	49
	50	! outputs:
	51	! --------
	52
	53	INTEGER lmin(ngrid) ! plume base level (first unstable level)
	54	INTEGER lalim(ngrid) ! higher alimentation level
	55	INTEGER lmix(ngrid) ! maximum vertical speed level
	56	INTEGER lmix_bis(ngrid) ! maximum vertical speed level (modified)
	57
	58	REAL alim_star(ngrid,klev) ! normalized alimentation
	59	REAL alim_star_tot(ngrid) ! integrated alimentation
	60
	61	REAL f0(ngrid) ! previous time step mass flux norm
	62
	63	REAL detr_star(ngrid,klev) ! normalized detrainment
	64	REAL entr_star(ngrid,klev) ! normalized entrainment
	65	REAL f_star(ngrid,klev+1) ! normalized mass flux
	66
	67	REAL ztva(ngrid,klev) ! TRPV plume (after mixing)
	68	REAL ztva_est(ngrid,klev) ! TRPV plume (before mixing)
	69	REAL ztla(ngrid,klev) ! TP plume
	70	REAL zqla(ngrid,klev) ! ql plume (after mixing)
	71	REAL zqta(ngrid,klev) ! qt plume
	72	REAL zha(ngrid,klev) ! TRPV plume
	73
	74	REAL w_est(ngrid,klev+1) ! updraft square vertical speed (before mixing)
	75	REAL zw2(ngrid,klev+1) ! updraft square vertical speed (after mixing)
	76
	77	REAL zqsatth(ngrid,klev) ! saturation vapor pressure (after mixing)
	78
	79	! local:
	80	! ------
	81
	82	INTEGER ig, l, k
	83	INTEGER lt
	84	INTEGER lm
	85
	86	REAL zqla_est(ngrid,klev) ! ql plume (before mixing)
	87	REAL zta_est(ngrid,klev) ! TR plume (before mixing)
	88	REAL zbuoy(ngrid,klev) ! B plume
	89	REAL zbuoyjam(ngrid,klev) ! B plume (modified)
	90
	91	REAL ztemp(ngrid) ! temperature for saturation vapor pressure computation in plume
	92	REAL zqsat(ngrid) ! saturation vapor pressure (before mixing)
	93	REAL zdz ! layers height
	94
	95	REAL zalpha !
	96	REAL zdqt(ngrid,klev) !
	97	REAL zbetalpha !
	98	REAL zdw2 !
	99	REAL zdw2bis !
	100	REAL zw2fact !
	101	REAL zw2factbis !
	102	REAL zw2m !
	103	REAL zdzbis !
	104	REAL d_temp(ngrid) !
	105	REAL coefzlmel !
	106	REAL zdz2 !
	107	REAL zdz3 !
	108	REAL lmel !
	109	REAL zlmel !
	110	REAL zlmelup !
	111	REAL zlmeldwn !
	112	REAL zlt !
	113	REAL zltdwn !
	114	REAL zltup ! useless here
	115
	116	LOGICAL active(ngrid) ! if the plume is active at ig,l (speed and incoming mass flux > 0 or l=lmin)
	117	LOGICAL activetmp(ngrid) ! if the plus is active at ig,l (active=true and outgoing mass flux > 0)
	118	LOGICAL, SAVE :: first = .true. ! if it is the first time step
	119
	120	!$OMP THREADPRIVATE(first)
	121
	122	!==============================================================================
	123	! Initialization
	124	!==============================================================================
	125
	126	zbetalpha = betalpha / (1. + betalpha)
	127
	128	ztva(:,:) = ztv(:,:) ! ztva is set to the virtual potential temperature withour latent heat release
	129	ztva_est(:,:) = ztva(:,:) ! ztva_est is set to the virtual potential temperature withour latent heat release
	130	ztla(:,:) = zthl(:,:) ! ztla is set to the potential temperature
	131	zqta(:,:) = po(:,:) ! zqta is set to qt
	132	zqla(:,:) = 0. ! zqla is set to ql
	133	zqla_est(:,:) = 0. ! zqla_est is set to ql
	134	zha(:,:) = ztva(:,:) ! zha is set to the plume virtual potential temperature withour latent heat release
	135
	136	zqsat(:) = 0.
	137	zqsatth(:,:) = 0.
	138
	139	w_est(:,:) = 0.
	140	zw2(:,:) = 0.
	141
	142	zbuoy(:,:) = 0.
	143	zbuoyjam(:,:) = 0.
	144
	145	f_star(:,:) = 0.
	146	detr_star(:,:) = 0.
	147	entr_star(:,:) = 0.
	148	alim_star(:,:) = 0.
	149	alim_star_tot(:) = 0.
	150
	151	lmix(:) = 1
	152	lmix_bis(:) = 2
	153	lalim(:) = 1
	154	lmin(:) = linf
	155
	156	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	157	! Pour activer un contraste de temperature a la base du panache
	158	! AB : It is used only in the thermcell_alim_init subroutine.
	159	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	160	d_temp(:) = 0.
	161
	162	!==============================================================================
	163	! 0. Calcul de l'alimentation
	164	!==============================================================================
	165
	166	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	167	! AB : Convective plumes can go off from every layer above the linf-th and
	168	! where pressure is lesser than pres_limit (cf. thermcell_plume).
	169	! The second constraint is added to avoid the parametrization occurs too
	170	! high when the low atmosphere is stable.
	171	! However, once there is a triggered plume, it can rise as high as its
	172	! velocity allows it (it can overshoot).
	173	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	174	DO ig=1,ngrid
	175	active(ig) = .false.
	176	l = linf
	177	DO WHILE ((.not.active(ig)) .and. pplev(ig,l+1).gt.pres_limit .and. l.lt.klev)
	178	IF (ztv(ig,l).gt.ztv(ig,l+1)) THEN
	179	active(ig) = .true.
	180	lmin(ig) = l
	181	ENDIF
	182	l = l + 1
	183	ENDDO
	184	ENDDO
	185
	186	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[2065]	187	! AB : On pourrait n'appeler thermcell_alim que si la plume est active
[2060]	188	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	189	CALL thermcell_alim(iflag_thermals_alim,ngrid,klev,ztv,d_temp, &
[2065]	190	& zlev,alim_star,lalim,lmin)
[2060]	191
	192	!==============================================================================
	193	! 1. Calcul dans la premiere couche
	194	!==============================================================================
	195
	196	DO ig=1,ngrid
	197	IF (active(ig)) THEN
	198
	199	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	200	! AB : plume takes the environment features for every layer below lmin.
	201	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	202	DO l=1,lmin(ig)
	203	ztla(ig,l) = zthl(ig,l)
	204	zqta(ig,l) = po(ig,l)
	205	zqla(ig,l) = zl(ig,l)
	206	ENDDO
	207
	208	l = lmin(ig)
	209	f_star(ig,l+1) = alim_star(ig,l)
	210
	211	zw2(ig,l+1) = 2. * RG * (zlev(ig,l+1) - zlev(ig,l)) &
	212	& * (ztv(ig,l) - ztv(ig,l+1)) / ztv(ig,l+1)
	213
	214	w_est(ig,l+1) = zw2(ig,l+1)
	215
	216	ENDIF
	217	ENDDO
	218
	219	!==============================================================================
	220	! 2. Boucle de calcul de la vitesse verticale dans le thermique
	221	!==============================================================================
	222
	223	DO l=2,klev-1
	224
	225	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	226	! AB : we decide here if the plume is still active or not. When the plume's
	227	! first level is reached, we set active to "true". Otherwise, it is given
[2065]	228	! by zw2, f_star, alim_star and entr_star.
[2060]	229	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	230	DO ig=1,ngrid
	231	IF (l==lmin(ig)+1) THEN
	232	active(ig) = .true.
	233	ENDIF
	234
	235	active(ig) = active(ig) &
	236	& .and. zw2(ig,l)>1.e-10 &
	237	& .and. f_star(ig,l)+alim_star(ig,l)+entr_star(ig,l)>1.e-10
	238	ENDDO
	239
	240	ztemp(:) = zpspsk(:,l) * ztla(:,l-1)
[2065]	241
[2060]	242	DO ig=1,ngrid
	243	CALL watersat(ztemp(ig), pplev(ig,l), zqsat(ig))
	244	ENDDO
	245
	246	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	247	! AB : we compute thermodynamical values and speed in the plume in the layer l
	248	! without mixing with environment.
	249	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	250
	251	DO ig=1,ngrid
	252	IF (active(ig)) THEN
	253	zqla_est(ig,l) = max(0.,zqta(ig,l-1)-zqsat(ig)) ! zqla_est set to ql plume
	254	ztva_est(ig,l) = ztla(ig,l-1)zpspsk(ig,l)+RLvCpzqla_est(ig,l) ! ztva_est set to TR plume
	255	zta_est(ig,l) = ztva_est(ig,l) ! zta_est set to TR plume
	256	ztva_est(ig,l) = ztva_est(ig,l)/zpspsk(ig,l) ! ztva_est set to TRP plume
	257	ztva_est(ig,l) = ztva_est(ig,l)(1.+RETV(zqta(ig,l-1) & ! ztva_est set to TRPV plume
	258	& - zqla_est(ig,l))-zqla_est(ig,l))
	259
	260	zbuoy(ig,l) = RG * (ztva_est(ig,l)-ztv(ig,l)) / ztv(ig,l)
	261	zdz = zlev(ig,l+1) - zlev(ig,l)
	262
	263	!==============================================================================
	264	! 3. Calcul de la flotabilite modifiee par melange avec l'air au dessus
	265	!==============================================================================
	266
	267	lmel = fact_thermals_ed_dz * zlev(ig,l)
	268	zlmel = zlev(ig,l) + lmel
	269	lt = l + 1
	270	zlt = zlev(ig,lt)
	271	zdz2 = zlev(ig,lt) - zlev(ig,l)
	272
	273	DO while (lmel.gt.zdz2)
	274	lt = lt + 1
	275	zlt = zlev(ig,lt)
	276	zdz2 = zlev(ig,lt) - zlev(ig,l)
	277	ENDDO
	278
	279	! IF (lt-l.gt.1) THEN
	280	! print *, 'WARNING: lt is greater than l+1!'
	281	! print *, 'l,lt', l, lt
	282	! ENDIF
	283
	284	zdz3 = zlev(ig,lt+1) - zlt
	285	zltdwn = zlev(ig,lt) - zdz3 / 2
	286	zlmelup = zlmel + (zdz / 2)
	287	coefzlmel = Min(1.,(zlmelup - zltdwn) / zdz)
	288
	289	zbuoyjam(ig,l) = 1.* RG * (coefzlmel * &
	290	& (ztva_est(ig,l) - ztv(ig,lt)) / ztv(ig,lt) &
	291	& + (1. - coefzlmel) * &
	292	& (ztva_est(ig,l) - ztv(ig,lt-1)) / ztv(ig,lt-1)) &
	293	& + 0. * zbuoy(ig,l)
	294
	295	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	296	! AB : initial formulae
	297	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	298	! zw2fact = fact_epsilon * 2. * zdz / (1. + betalpha)
	299	! zdw2 = afact * zbuoy(ig,l) / fact_epsilon
	300	! zdw2bis = afact * zbuoy(ig,l-1) / fact_epsilon
	301	! w_est(ig,l+1) = Max(0.0001,exp(-zw2fact)*(w_est(ig,l)-zdw2)+zdw2)
	302	! w_est(ig,l+1) = Max(0.0001,exp(-zw2fact)*(w_est(ig,l)-zdw2bis)+zdw2)
	303	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[2065]	304	! AB : own derivation for w_est (Rio 2010 formula with e=d=0)
[2060]	305	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	306	zw2fact = 2. * fact_epsilon * zdz
	307	zdw2 = 2. * afact * zbuoy(ig,l) * zdz
	308	w_est(ig,l+1) = Max(0., exp(-zw2fact) * w_est(ig,l) + zdw2)
	309
	310	! IF (w_est(ig,l+1).le.0.) THEN
	311	! print *, 'WARNING: w_est is negative!'
	312	! print *, 'l,w_est', l+1, w_est(ig,l+1)
	313	! ENDIF
	314	ENDIF
	315	ENDDO
	316
	317	!==============================================================================
	318	! 4. Calcul de l'entrainement et du detrainement
	319	!==============================================================================
	320
	321	DO ig=1,ngrid
	322	IF (active(ig)) THEN
	323
	324	zdz = zlev(ig,l+1) - zlev(ig,l)
[2065]	325	zalpha = f0(ig) * f_star(ig,l) / sqrt(w_est(ig,l+1)) / rhobarz(ig,l)
[2060]	326
	327	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[2065]	328	! AB : The next test is added to avoid a division by zero when w_est becomes
	329	! negative.
	330	! Indeed, entr and detr computed here are of no importance because w_est
	331	! <= 0 means it will be the last layer reached by the plume and then they
	332	! will be reset in thermcell_flux.
[2060]	333	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[2065]	334	IF (w_est(ig,l+1).eq.0.) THEN
	335	zw2m = 1.
[2060]	336	ELSE
[2065]	337	zw2m = w_est(ig,l+1)
[2060]	338	ENDIF
	339
	340	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[2065]	341	! AB : The next test is added to avoid a division by zero if there is no water
	342	! in the environment.
[2060]	343	! In the case where there is no water in the env. but water in the plume
	344	! (ascending from depth) we set the effect on detrainment equal to zero
	345	! but at the next time step, po will be positive thanks to the mixing and
	346	! then the physical effect of the water gradient will be taken on board.
	347	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	348	IF (po(ig,l).lt.1.e-6) THEN
[2065]	349	! print *, 'WARNING: po=0 in layer',l,'!'
[2060]	350	! print *, 'po,zqta', po(ig,l), zqta(ig,l-1)
	351	zdqt(ig,l) = 0.0
	352	ELSE
	353	zdqt(ig,l) = max(zqta(ig,l-1)-po(ig,l),0.) / po(ig,l)
	354	ENDIF
	355
	356	!------------------------------------------------------------------------------
	357	! Detrainment
	358	!------------------------------------------------------------------------------
	359
	360	detr_star(ig,l) = f_star(ig,l) * zdz * ( &
	361	& mix0 * 0.1 / (zalpha + 0.001) &
	362	& + MAX(detr_min, &
	363	& -afact * zbetalpha * zbuoyjam(ig,l) / zw2m &
	364	& + detr_q_coef(zdqt(ig,l)/zw2m)*detr_q_power) )
	365
	366	! IF (detr_star(ig,l).lt.0.) THEN
	367	! print *, 'WARNING: detrainment is negative!'
	368	! print *, 'l,detr', l, detr_star(ig,l)
	369	! ENDIF
	370
	371	!------------------------------------------------------------------------------
	372	! Entrainment
	373	!------------------------------------------------------------------------------
	374
	375	entr_star(ig,l) = f_star(ig,l) * zdz * ( &
	376	& mix0 * 0.1 / (zalpha+0.001) &
	377	& + MAX(entr_min, &
	378	& zbetalpha * afact * zbuoyjam(ig,l) / zw2m &
	379	& - zbetalpha * fact_epsilon) )
	380
	381	! IF (entr_star(ig,l).lt.0.) THEN
	382	! print *, 'WARNING: entrainment is negative!'
	383	! print *, 'l,entr', l, entr_star(ig,l)
	384	! ENDIF
	385
	386	!------------------------------------------------------------------------------
	387	! Alimentation and entrainment
	388	!------------------------------------------------------------------------------
	389
	390	IF (l.lt.lalim(ig)) THEN
	391	alim_star(ig,l) = max(alim_star(ig,l),entr_star(ig,l))
	392	entr_star(ig,l) = 0.
	393	ENDIF
	394
	395	!------------------------------------------------------------------------------
	396	! Mass flux
	397	!------------------------------------------------------------------------------
	398
	399	f_star(ig,l+1) = f_star(ig,l) + alim_star(ig,l) &
	400	& + entr_star(ig,l) - detr_star(ig,l)
	401
	402	! IF (f_star(ig,l+1).le.0.) THEN
	403	! print *, 'WARNING: mass flux is negative!'
	404	! print *, 'l,f_star', l+1, f_star(ig,l+1)
	405	! ENDIF
	406
	407	ENDIF
	408	ENDDO
	409
	410	!==============================================================================
	411	! 5. Calcul de la vitesse verticale en melangeant Tl et qt du thermique
	412	!==============================================================================
	413
	414	activetmp(:) = active(:) .and. f_star(:,l+1)>1.e-10
	415
	416	!------------------------------------------------------------------------------
	417	! Calcul du melange avec l'environnement
	418	!------------------------------------------------------------------------------
	419
	420	DO ig=1,ngrid
	421	IF (activetmp(ig)) THEN
[2065]	422	ztla(ig,l) = (f_star(ig,l) * ztla(ig,l-1) & ! ztla is set to TP in plume (mixed)
	423	& + (alim_star(ig,l) + entr_star(ig,l)) * zthl(ig,l)) &
	424	& / (f_star(ig,l+1) + detr_star(ig,l))
	425	zqta(ig,l) = (f_star(ig,l) * zqta(ig,l-1) + & ! zqta is set to qt in plume (mixed)
	426	& + (alim_star(ig,l) + entr_star(ig,l)) * po(ig,l)) &
	427	& / (f_star(ig,l+1) + detr_star(ig,l))
[2060]	428	ENDIF
	429	ENDDO
	430
	431	ztemp(:) = zpspsk(:,l) * ztla(:,l)
	432
	433	DO ig=1,ngrid
	434	IF (activetmp(ig)) THEN
	435	CALL watersat(ztemp(ig), pplev(ig,l), zqsatth(ig,l))
	436	ENDIF
	437	ENDDO
	438
	439	!------------------------------------------------------------------------------
	440	! Calcul de la vitesse verticale zw2 apres le melange
	441	!------------------------------------------------------------------------------
	442
	443	DO ig=1,ngrid
	444	IF (activetmp(ig)) THEN
	445	zqla(ig,l) = max(0.,zqta(ig,l)-zqsatth(ig,l)) ! zqla is set to ql plume (mixed)
	446	ztva(ig,l) = ztla(ig,l) * zpspsk(ig,l)+RLvCp*zqla(ig,l) ! ztva is set to TR plume (mixed)
	447	ztva(ig,l) = ztva(ig,l) / zpspsk(ig,l) ! ztva is set to TRP plume (mixed)
	448	zha(ig,l) = ztva(ig,l) ! zha is set to TRP plume (mixed)
	449	ztva(ig,l) = ztva(ig,l) * (1. + RETV*(zqta(ig,l)-zqla(ig,l)) & ! ztva is set to TRPV plume (mixed)
	450	& - zqla(ig,l))
	451
	452	zbuoy(ig,l) = RG * (ztva(ig,l) - ztv(ig,l)) / ztv(ig,l)
	453	zdz = zlev(ig,l+1) - zlev(ig,l)
	454
	455	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[2065]	456	! AB : initial formula
[2060]	457	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	458	! zw2fact = fact_epsilon * 2. * zdz / (1. + betalpha)
	459	! zdw2 = afact * zbuoy(ig,l) / fact_epsilon
	460	! zdw2bis = afact * zbuoy(ig,l-1) / fact_epsilon
	461	! zw2(ig,l+1) = Max(0.0001,exp(-zw2fact)*(zw2(ig,l)-zdw2)+zdw2)
	462	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[2065]	463	! AB : own derivation for zw2 (Rio 2010 formula)
[2060]	464	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	465	zw2fact = 2. * (fact_epsilon * zdz + entr_star(ig,l) / f_star(ig,l))
	466	zdw2 = 2. * afact * zbuoy(ig,l) * zdz
	467	zw2(ig,l+1) = Max(0., exp(-zw2fact) * zw2(ig,l) + zdw2)
	468
	469	! IF (zw2(ig,l+1).le.0.) THEN
	470	! print *, 'WARNING: zw2 is negative!'
	471	! print *, 'l,zw2', l+1, zw2(ig,l+1)
	472	! ENDIF
	473	ENDIF
	474	ENDDO
	475
	476	ENDDO
	477
	478	!==============================================================================
	479	! 6. New computation of alim_star_tot
	480	!==============================================================================
	481
	482	DO ig=1,ngrid
[2065]	483	alim_star_tot(ig) = 0.
[2060]	484	ENDDO
	485
	486	DO ig=1,ngrid
	487	DO l=1,lalim(ig)-1
	488	alim_star_tot(ig) = alim_star_tot(ig) + alim_star(ig,l)
	489	ENDDO
	490	ENDDO
	491
	492	#undef wrgrads_thermcell
	493	#ifdef wrgrads_thermcell
	494	call wrgradsfi(1,klev,entr_star(igout,1:klev),'esta ','esta ')
	495	call wrgradsfi(1,klev,detr_star(igout,1:klev),'dsta ','dsta ')
	496	call wrgradsfi(1,klev,zbuoy(igout,1:klev) ,'buoy ','buoy ')
	497	call wrgradsfi(1,klev,zdqt(igout,1:klev) ,'dqt ','dqt ')
	498	call wrgradsfi(1,klev,w_est(igout,1:klev) ,'w_est ','w_est ')
	499	call wrgradsfi(1,klev,w_est(igout,2:klev+1) ,'w_es2 ','w_es2 ')
	500	call wrgradsfi(1,klev,zw2(igout,1:klev) ,'zw2A ','zw2A ')
	501	#endif
	502
	503
	504	RETURN
	505	END

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

Download in other formats:

Original Format