Context Navigation

clmain.F @ 757

Last change on this file since 757 was 686, checked in by lmdzadmin, 19 years ago
Appel du slab lorsque nisurf=is_sic IM
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 75.6 KB

Line
1	!
2	! $Header$
3	!
4	c
5	c
6	SUBROUTINE clmain(dtime,itap,date0,pctsrf,pctsrf_new,
7	. t,q,u,v,
8	. jour, rmu0, co2_ppm,
9	. ok_veget, ocean, npas, nexca, ts,
10	. soil_model,cdmmax, cdhmax,
11	. ksta, ksta_ter, ok_kzmin, ftsoil,qsol,
12	cIM BAD . paprs,pplay,radsol,snow,qsurf,evap,albe,alblw,
13	. paprs,pplay,snow,qsurf,evap,albe,alblw,
14	. fluxlat,
15	. rain_f, snow_f, solsw, sollw, sollwdown, fder,
16	. rlon, rlat, cufi, cvfi, rugos,
17	. debut, lafin, agesno,rugoro,
18	. d_t,d_q,d_u,d_v,d_ts,
19	. flux_t,flux_q,flux_u,flux_v,cdragh,cdragm,
20	. q2,
21	. dflux_t,dflux_q,
22	. zcoefh,zu1,zv1, t2m, q2m, u10m, v10m,
23	cIM cf. AM : pbl
24	. pblh,capCL,oliqCL,cteiCL,pblT,
25	. therm,trmb1,trmb2,trmb3,plcl,
26	. fqcalving,ffonte, run_off_lic_0,
27	cIM "slab" ocean
28	. flux_o, flux_g, tslab, seaice)
29	cAA . itr, tr, flux_surf, d_tr)
30	cAA REM:
31	cAA-----
32	cAA Tout ce qui a trait au traceurs est dans phytrac maintenant
33	cAA pour l'instant le calcul de la couche limite pour les traceurs
34	cAA se fait avec cltrac et ne tient pas compte de la differentiation
35	cAA des sous-fraction de sol.
36	cAA REM bis :
37	cAA----------
38	cAA Pour pouvoir extraire les coefficient d'echanges et le vent
39	cAA dans la premiere couche, 3 champs supplementaires ont ete crees
40	cAA zcoefh,zu1 et zv1. Pour l'instant nous avons moyenne les valeurs
41	cAA de ces trois champs sur les 4 subsurfaces du modele. Dans l'avenir
42	cAA si les informations des subsurfaces doivent etre prises en compte
43	cAA il faudra sortir ces memes champs en leur ajoutant une dimension,
44	cAA c'est a dire nbsrf (nbre de subsurface).
45	USE ioipsl
46	USE interface_surf
47	IMPLICIT none
48	c======================================================================
49	c Auteur(s) Z.X. Li (LMD/CNRS) date: 19930818
50	c Objet: interface de "couche limite" (diffusion verticale)
51	c Arguments:
52	c dtime----input-R- interval du temps (secondes)
53	c itap-----input-I- numero du pas de temps
54	c date0----input-R- jour initial
55	c t--------input-R- temperature (K)
56	c q--------input-R- vapeur d'eau (kg/kg)
57	c u--------input-R- vitesse u
58	c v--------input-R- vitesse v
59	c ts-------input-R- temperature du sol (en Kelvin)
60	c paprs----input-R- pression a intercouche (Pa)
61	c pplay----input-R- pression au milieu de couche (Pa)
62	c radsol---input-R- flux radiatif net (positif vers le sol) en W/m**2
63	c rlat-----input-R- latitude en degree
64	c rugos----input-R- longeur de rugosite (en m)
65	c cufi-----input-R- resolution des mailles en x (m)
66	c cvfi-----input-R- resolution des mailles en y (m)
67	c
68	c d_t------output-R- le changement pour "t"
69	c d_q------output-R- le changement pour "q"
70	c d_u------output-R- le changement pour "u"
71	c d_v------output-R- le changement pour "v"
72	c d_ts-----output-R- le changement pour "ts"
73	c flux_t---output-R- flux de chaleur sensible (CpT) J/m2/s (W/m2)
74	c (orientation positive vers le bas)
75	c flux_q---output-R- flux de vapeur d'eau (kg/m**2/s)
76	c flux_u---output-R- tension du vent X: (kg m/s)/(m**2 s) ou Pascal
77	c flux_v---output-R- tension du vent Y: (kg m/s)/(m**2 s) ou Pascal
78	c dflux_t derive du flux sensible
79	c dflux_q derive du flux latent
80	cIM "slab" ocean
81	c flux_g---output-R- flux glace (pour OCEAN='slab ')
82	c flux_o---output-R- flux ocean (pour OCEAN='slab ')
83	c tslab-in/output-R temperature du slab ocean (en Kelvin) ! uniqmnt pour slab
84	c seaice---output-R- glace de mer (kg/m2) (pour OCEAN='slab ')
85	ccc
86	c ffonte----Flux thermique utilise pour fondre la neige
87	c fqcalving-Flux d'eau "perdue" par la surface et necessaire pour limiter la
88	c hauteur de neige, en kg/m2/s
89	cAA on rajoute en output yu1 et yv1 qui sont les vents dans
90	cAA la premiere couche
91	cAA ces 4 variables sont maintenant traites dans phytrac
92	c itr--------input-I- nombre de traceurs
93	c tr---------input-R- q. de traceurs
94	c flux_surf--input-R- flux de traceurs a la surface
95	c d_tr-------output-R tendance de traceurs
96	cIM cf. AM : PBL
97	c trmb1-------deep_cape
98	c trmb2--------inhibition
99	c trmb3-------Point Omega
100	c Cape(klon)-------Cape du thermique
101	c EauLiq(klon)-------Eau liqu integr du thermique
102	c ctei(klon)-------Critere d'instab d'entrainmt des nuages de CL
103	c lcl------- Niveau de condensation
104	c pblh------- HCL
105	c pblT------- T au nveau HCL
106	c======================================================================
107	#include "dimensions.h"
108	#include "dimphy.h"
109	#include "indicesol.h"
110	c$$$ PB ajout pour soil
111	#include "dimsoil.h"
112	#include "iniprint.h"
113	#include "compbl.h"
114	c
115	REAL dtime
116	real date0
117	integer itap
118	REAL t(klon,klev), q(klon,klev)
119	REAL u(klon,klev), v(klon,klev)
120	cIM 230604 BAD REAL radsol(klon) ???
121	REAL paprs(klon,klev+1), pplay(klon,klev)
122	REAL rlon(klon), rlat(klon), cufi(klon), cvfi(klon)
123	REAL d_t(klon, klev), d_q(klon, klev)
124	REAL d_u(klon, klev), d_v(klon, klev)
125	REAL flux_t(klon,klev, nbsrf), flux_q(klon,klev, nbsrf)
126	REAL dflux_t(klon), dflux_q(klon)
127	cIM "slab" ocean
128	REAL flux_o(klon), flux_g(klon)
129	REAL y_flux_o(klon), y_flux_g(klon)
130	REAL tslab(klon), ytslab(klon)
131	REAL seaice(klon), y_seaice(klon)
132	c
133	REAL amn, amx
134	cIM cf JLD
135	REAL y_fqcalving(klon), y_ffonte(klon)
136	REAL fqcalving(klon,nbsrf), ffonte(klon,nbsrf)
137	REAL run_off_lic_0(klon), y_run_off_lic_0(klon)
138
139	REAL flux_u(klon,klev, nbsrf), flux_v(klon,klev, nbsrf)
140	REAL rugmer(klon), agesno(klon,nbsrf),rugoro(klon)
141	REAL cdragh(klon), cdragm(klon)
142	integer jour ! jour de l'annee en cours
143	real rmu0(klon) ! cosinus de l'angle solaire zenithal
144	REAL co2_ppm ! taux CO2 atmosphere
145	LOGICAL debut, lafin, ok_veget
146	character*6 ocean
147	integer npas, nexca
148	cAA INTEGER itr
149	cAA REAL tr(klon,klev,nbtr)
150	cAA REAL d_tr(klon,klev,nbtr)
151	cAA REAL flux_surf(klon,nbtr)
152	c
153	REAL pctsrf(klon,nbsrf)
154	REAL ts(klon,nbsrf)
155	REAL d_ts(klon,nbsrf)
156	REAL snow(klon,nbsrf)
157	REAL qsurf(klon,nbsrf)
158	REAL evap(klon,nbsrf)
159	REAL albe(klon,nbsrf)
160	REAL alblw(klon,nbsrf)
161	c$$$ PB
162	REAL fluxlat(klon,nbsrf)
163	C
164	real rain_f(klon), snow_f(klon)
165	REAL fder(klon)
166	cIM cf. JLD REAL sollw(klon), solsw(klon), sollwdown(klon)
167	REAL sollw(klon,nbsrf), solsw(klon,nbsrf), sollwdown(klon)
168	REAL rugos(klon,nbsrf)
169	C la nouvelle repartition des surfaces sortie de l'interface
170	REAL pctsrf_new(klon,nbsrf)
171	cAA
172	REAL zcoefh(klon,klev)
173	REAL zu1(klon)
174	REAL zv1(klon)
175	cAA
176	c$$$ PB ajout pour soil
177	LOGICAL soil_model
178	cIM ajout seuils cdrm, cdrh
179	REAL cdmmax, cdhmax
180	cIM: 261103
181	REAL ksta, ksta_ter
182	LOGICAL ok_kzmin
183	cIM: 261103
184	REAL ftsoil(klon,nsoilmx,nbsrf)
185	REAL ytsoil(klon,nsoilmx)
186	REAL qsol(klon)
187	c======================================================================
188	EXTERNAL clqh, clvent, coefkz, calbeta, cltrac
189	c======================================================================
190	REAL yts(klon), yrugos(klon), ypct(klon), yz0_new(klon)
191	REAL yalb(klon)
192	REAL yalblw(klon)
193	REAL yu1(klon), yv1(klon)
194	real ysnow(klon), yqsurf(klon), yagesno(klon), yqsol(klon)
195	real yrain_f(klon), ysnow_f(klon)
196	real ysollw(klon), ysolsw(klon), ysollwdown(klon)
197	real yfder(klon), ytaux(klon), ytauy(klon)
198	REAL yrugm(klon), yrads(klon),yrugoro(klon)
199	c$$$ PB
200	REAL yfluxlat(klon)
201	C
202	REAL y_d_ts(klon)
203	REAL y_d_t(klon, klev), y_d_q(klon, klev)
204	REAL y_d_u(klon, klev), y_d_v(klon, klev)
205	REAL y_flux_t(klon,klev), y_flux_q(klon,klev)
206	REAL y_flux_u(klon,klev), y_flux_v(klon,klev)
207	REAL y_dflux_t(klon), y_dflux_q(klon)
208	REAL ycoefh(klon,klev), ycoefm(klon,klev)
209	REAL yu(klon,klev), yv(klon,klev)
210	REAL yt(klon,klev), yq(klon,klev)
211	REAL ypaprs(klon,klev+1), ypplay(klon,klev), ydelp(klon,klev)
212	cAA REAL ytr(klon,klev,nbtr)
213	cAA REAL y_d_tr(klon,klev,nbtr)
214	cAA REAL yflxsrf(klon,nbtr)
215	c
216	c
217	LOGICAL ok_nonloc
218	PARAMETER (ok_nonloc=.FALSE.)
219	REAL ycoefm0(klon,klev), ycoefh0(klon,klev)
220
221	cIM 081204 hcl_Anne ? BEG
222	real yzlay(klon,klev),yzlev(klon,klev+1),yteta(klon,klev)
223	real ykmm(klon,klev+1),ykmn(klon,klev+1)
224	real ykmq(klon,klev+1)
225	real yq2(klon,klev+1),q2(klon,klev+1,nbsrf)
226	real q2diag(klon,klev+1)
227	cIM 081204 real yustar(klon),y_cd_m(klon),y_cd_h(klon)
228	cIM 081204 hcl_Anne ? END
229	c
230	#include "YOMCST.h"
231	#include "YOETHF.h"
232	#include "FCTTRE.h"
233	REAL u1lay(klon), v1lay(klon)
234	REAL delp(klon,klev)
235	INTEGER i, k, nsrf
236	cAA INTEGER it
237	INTEGER ni(klon), knon, j
238	c Introduction d'une variable "pourcentage potentiel" pour tenir compte
239	c des eventuelles apparitions et/ou disparitions de la glace de mer
240	REAL pctsrf_pot(klon,nbsrf)
241
242	c======================================================================
243	REAL zx_alf1, zx_alf2 !valeur ambiante par extrapola.
244	c======================================================================
245	c
246	c maf pour sorties IOISPL en cas de debugagage
247	c
248	CHARACTER*80 cldebug
249	SAVE cldebug
250	CHARACTER*8 cl_surf(nbsrf)
251	SAVE cl_surf
252	INTEGER nhoridbg, nidbg
253	SAVE nhoridbg, nidbg
254	INTEGER ndexbg(iim*(jjm+1))
255	REAL zx_lon(iim,jjm+1), zx_lat(iim,jjm+1), zjulian
256	REAL tabindx(klon)
257	REAL debugtab(iim,jjm+1)
258	LOGICAL first_appel
259	SAVE first_appel
260	DATA first_appel/.true./
261	LOGICAL debugindex
262	SAVE debugindex
263	DATA debugindex/.false./
264	integer idayref
265	#include "temps.h"
266	REAL t2m(klon,nbsrf), q2m(klon,nbsrf)
267	REAL u10m(klon,nbsrf), v10m(klon,nbsrf)
268	c
269	REAL yt2m(klon), yq2m(klon), yu10m(klon)
270	REAL yustar(klon)
271	c -- LOOP
272	REAL yu10mx(klon)
273	REAL yu10my(klon)
274	REAL ywindsp(klon)
275	c -- LOOP
276	c
277	REAL yt10m(klon), yq10m(klon)
278	cIM cf. AM : pbl, hbtm2 (Comme les autres diagnostics on cumule ds physic ce qui
279	c permet de sortir les grdeurs par sous surface)
280	REAL pblh(klon,nbsrf)
281	REAL plcl(klon,nbsrf)
282	REAL capCL(klon,nbsrf)
283	REAL oliqCL(klon,nbsrf)
284	REAL cteiCL(klon,nbsrf)
285	REAL pblT(klon,nbsrf)
286	REAL therm(klon,nbsrf)
287	REAL trmb1(klon,nbsrf)
288	REAL trmb2(klon,nbsrf)
289	REAL trmb3(klon,nbsrf)
290	REAL ypblh(klon)
291	REAL ylcl(klon)
292	REAL ycapCL(klon)
293	REAL yoliqCL(klon)
294	REAL ycteiCL(klon)
295	REAL ypblT(klon)
296	REAL ytherm(klon)
297	REAL ytrmb1(klon)
298	REAL ytrmb2(klon)
299	REAL ytrmb3(klon)
300	REAL y_cd_h(klon), y_cd_m(klon)
301	c REAL ygamt(klon,2:klev) ! contre-gradient pour temperature
302	c REAL ygamq(klon,2:klev) ! contre-gradient pour humidite
303	REAL uzon(klon), vmer(klon)
304	REAL tair1(klon), qair1(klon), tairsol(klon)
305	REAL psfce(klon), patm(klon)
306	c
307	REAL qairsol(klon), zgeo1(klon)
308	REAL rugo1(klon)
309	c
310	LOGICAL zxli ! utiliser un jeu de fonctions simples
311	PARAMETER (zxli=.FALSE.)
312	c
313	REAL zt, zqs, zdelta, zcor
314	REAL t_coup
315	PARAMETER(t_coup=273.15)
316	C
317	character (len = 20) :: modname = 'clmain'
318	LOGICAL check
319	PARAMETER (check=.false.)
320
321
322	c initialisation Anne
323	ytherm(:) = 0.
324	C
325	if (check) THEN
326	write(,) modname,' klon=',klon
327	call flush(6)
328	endif
329	IF (debugindex .and. first_appel) THEN
330	first_appel=.false.
331	!
332	! initialisation sorties netcdf
333	!
334	idayref = day_ini
335	CALL ymds2ju(annee_ref, 1, idayref, 0.0, zjulian)
336	CALL gr_fi_ecrit(1,klon,iim,jjm+1,rlon,zx_lon)
337	DO i = 1, iim
338	zx_lon(i,1) = rlon(i+1)
339	zx_lon(i,jjm+1) = rlon(i+1)
340	ENDDO
341	CALL gr_fi_ecrit(1,klon,iim,jjm+1,rlat,zx_lat)
342	cldebug='sous_index'
343	CALL histbeg(cldebug, iim,zx_lon(:,1),jjm+1,zx_lat(1,:),
344	$ 1,iim,1,jjm
345	$ +1, itau_phy,zjulian,dtime,nhoridbg,nidbg)
346	! no vertical axis
347	cl_surf(1)='ter'
348	cl_surf(2)='lic'
349	cl_surf(3)='oce'
350	cl_surf(4)='sic'
351	DO nsrf=1,nbsrf
352	CALL histdef(nidbg, cl_surf(nsrf),cl_surf(nsrf), "-",iim,
353	$ jjm+1,nhoridbg, 1, 1, 1, -99, 32, "inst", dtime,dtime)
354	END DO
355	CALL histend(nidbg)
356	CALL histsync(nidbg)
357	ENDIF
358
359	DO k = 1, klev ! epaisseur de couche
360	DO i = 1, klon
361	delp(i,k) = paprs(i,k)-paprs(i,k+1)
362	ENDDO
363	ENDDO
364	DO i = 1, klon ! vent de la premiere couche
365	ccc zx_alf1 = (paprs(i,1)-pplay(i,2))/(pplay(i,1)-pplay(i,2))
366	zx_alf1 = 1.0
367	zx_alf2 = 1.0 - zx_alf1
368	u1lay(i) = u(i,1)zx_alf1 + u(i,2)zx_alf2
369	v1lay(i) = v(i,1)zx_alf1 + v(i,2)zx_alf2
370	ENDDO
371	c
372	c initialisation:
373	c
374	DO i = 1, klon
375	rugmer(i) = 0.0
376	cdragh(i) = 0.0
377	cdragm(i) = 0.0
378	dflux_t(i) = 0.0
379	dflux_q(i) = 0.0
380	zu1(i) = 0.0
381	zv1(i) = 0.0
382	ENDDO
383	ypct = 0.0
384	yts = 0.0
385	ysnow = 0.0
386	yqsurf = 0.0
387	yalb = 0.0
388	yalblw = 0.0
389	yrain_f = 0.0
390	ysnow_f = 0.0
391	yfder = 0.0
392	ytaux = 0.0
393	ytauy = 0.0
394	ysolsw = 0.0
395	ysollw = 0.0
396	ysollwdown = 0.0
397	yrugos = 0.0
398	yu1 = 0.0
399	yv1 = 0.0
400	yrads = 0.0
401	ypaprs = 0.0
402	ypplay = 0.0
403	ydelp = 0.0
404	yu = 0.0
405	yv = 0.0
406	yt = 0.0
407	yq = 0.0
408	pctsrf_new = 0.0
409	y_flux_u = 0.0
410	y_flux_v = 0.0
411	C$$ PB
412	y_dflux_t = 0.0
413	y_dflux_q = 0.0
414	ytsoil = 999999.
415	yrugoro = 0.
416	c -- LOOP
417	yu10mx = 0.0
418	yu10my = 0.0
419	ywindsp = 0.0
420	c -- LOOP
421	DO nsrf = 1, nbsrf
422	DO i = 1, klon
423	d_ts(i,nsrf) = 0.0
424	ENDDO
425	END DO
426	C§§§ PB
427	yfluxlat=0.
428	flux_t = 0.
429	flux_q = 0.
430	flux_u = 0.
431	flux_v = 0.
432	DO k = 1, klev
433	DO i = 1, klon
434	d_t(i,k) = 0.0
435	d_q(i,k) = 0.0
436	c$$$ flux_t(i,k) = 0.0
437	c$$$ flux_q(i,k) = 0.0
438	d_u(i,k) = 0.0
439	d_v(i,k) = 0.0
440	c$$$ flux_u(i,k) = 0.0
441	c$$$ flux_v(i,k) = 0.0
442	zcoefh(i,k) = 0.0
443	ENDDO
444	ENDDO
445	cAA IF (itr.GE.1) THEN
446	cAA DO it = 1, itr
447	cAA DO k = 1, klev
448	cAA DO i = 1, klon
449	cAA d_tr(i,k,it) = 0.0
450	cAA ENDDO
451	cAA ENDDO
452	cAA ENDDO
453	cAA ENDIF
454
455	c
456	c Boucler sur toutes les sous-fractions du sol:
457	c
458	C Initialisation des "pourcentages potentiels". On considere ici qu'on
459	C peut avoir potentiellementdela glace sur tout le domaine oceanique
460	C (a affiner)
461
462	pctsrf_pot = pctsrf
463	pctsrf_pot(:,is_oce) = 1. - zmasq(:)
464	pctsrf_pot(:,is_sic) = 1. - zmasq(:)
465
466	DO 99999 nsrf = 1, nbsrf
467
468	c chercher les indices:
469	DO j = 1, klon
470	ni(j) = 0
471	ENDDO
472	knon = 0
473	DO i = 1, klon
474
475	C pour determiner le domaine a traiter on utilise les surfaces "potentielles"
476	C
477	IF (pctsrf_pot(i,nsrf).GT.epsfra) THEN
478	knon = knon + 1
479	ni(knon) = i
480	ENDIF
481	ENDDO
482	c
483	if (check) THEN
484	write(,)'CLMAIN, nsrf, knon =',nsrf, knon
485	call flush(6)
486	endif
487	c
488	c variables pour avoir une sortie IOIPSL des INDEX
489	c
490	IF (debugindex) THEN
491	tabindx(:)=0.
492	c tabindx(1:knon)=(/FLOAT(i),i=1:knon/)
493	DO i=1,knon
494	tabindx(1:knon)=FLOAT(i)
495	END DO
496	debugtab(:,:)=0.
497	ndexbg(:)=0
498	CALL gath2cpl(tabindx,debugtab,klon,knon,iim,jjm,ni)
499	CALL histwrite(nidbg,cl_surf(nsrf),itap,debugtab,iim*(jjm+1)
500	$ ,ndexbg)
501	ENDIF
502	IF (knon.EQ.0) GOTO 99999
503	DO j = 1, knon
504	i = ni(j)
505	ypct(j) = pctsrf(i,nsrf)
506	yts(j) = ts(i,nsrf)
507	cIM "slab" ocean
508	c PRINT *, 'tslab = ', i, tslab(i)
509	ytslab(i) = tslab(i)
510	y_seaice(i) = seaice(i)
511	c
512	ysnow(j) = snow(i,nsrf)
513	yqsurf(j) = qsurf(i,nsrf)
514	yalb(j) = albe(i,nsrf)
515	yalblw(j) = alblw(i,nsrf)
516	yrain_f(j) = rain_f(i)
517	ysnow_f(j) = snow_f(i)
518	yagesno(j) = agesno(i,nsrf)
519	yfder(j) = fder(i)
520	ytaux(j) = flux_u(i,1,nsrf)
521	ytauy(j) = flux_v(i,1,nsrf)
522	ysolsw(j) = solsw(i,nsrf)
523	ysollw(j) = sollw(i,nsrf)
524	ysollwdown(j) = sollwdown(i)
525	yrugos(j) = rugos(i,nsrf)
526	yrugoro(j) = rugoro(i)
527	yu1(j) = u1lay(i)
528	yv1(j) = v1lay(i)
529	yrads(j) = ysolsw(j)+ ysollw(j)
530	ypaprs(j,klev+1) = paprs(i,klev+1)
531	y_run_off_lic_0(j) = run_off_lic_0(i)
532	c -- LOOP
533	yu10mx(j) = u10m(i,nsrf)
534	yu10my(j) = v10m(i,nsrf)
535	ywindsp(j) = SQRT(yu10mx(j)yu10mx(j) + yu10my(j)yu10my(j) )
536	c -- LOOP
537	END DO
538	C
539	C IF bucket model for continent, copy soil water content
540	IF ( nsrf .eq. is_ter .and. .not. ok_veget ) THEN
541	DO j = 1, knon
542	i = ni(j)
543	yqsol(j) = qsol(i)
544	END DO
545	ELSE
546	yqsol(:)=0.
547	ENDIF
548	c$$$ PB ajour pour soil
549	DO k = 1, nsoilmx
550	DO j = 1, knon
551	i = ni(j)
552	ytsoil(j,k) = ftsoil(i,k,nsrf)
553	END DO
554	END DO
555	DO k = 1, klev
556	DO j = 1, knon
557	i = ni(j)
558	ypaprs(j,k) = paprs(i,k)
559	ypplay(j,k) = pplay(i,k)
560	ydelp(j,k) = delp(i,k)
561	yu(j,k) = u(i,k)
562	yv(j,k) = v(i,k)
563	yt(j,k) = t(i,k)
564	yq(j,k) = q(i,k)
565	ENDDO
566	ENDDO
567	c
568	c
569	c calculer Cdrag et les coefficients d'echange
570	CALL coefkz(nsrf, knon, ypaprs, ypplay,
571	cIM 261103
572	. ksta, ksta_ter,
573	cIM 261103
574	. yts, yrugos, yu, yv, yt, yq,
575	. yqsurf,
576	. ycoefm, ycoefh)
577	cIM 081204 BEG
578	cCR test
579	if (iflag_pbl.eq.1) then
580	cIM 081204 END
581	CALL coefkz2(nsrf, knon, ypaprs, ypplay,yt,
582	. ycoefm0, ycoefh0)
583	DO k = 1, klev
584	DO i = 1, knon
585	ycoefm(i,k) = MAX(ycoefm(i,k),ycoefm0(i,k))
586	ycoefh(i,k) = MAX(ycoefh(i,k),ycoefh0(i,k))
587	ENDDO
588	ENDDO
589	endif
590	c
591	cIM cf JLD : on seuille ycoefm et ycoefh
592	if (nsrf.eq.is_oce) then
593	do j=1,knon
594	c ycoefm(j,1)=min(ycoefm(j,1),1.1E-3)
595	ycoefm(j,1)=min(ycoefm(j,1),cdmmax)
596	c ycoefh(j,1)=min(ycoefh(j,1),1.1E-3)
597	ycoefh(j,1)=min(ycoefh(j,1),cdhmax)
598	enddo
599	endif
600
601	c
602	cIM: 261103
603	if (ok_kzmin) THEN
604	cIM cf FH: 201103 BEG
605	c Calcul d'une diffusion minimale pour les conditions tres stables.
606	call coefkzmin(knon,ypaprs,ypplay,yu,yv,yt,yq,ycoefm
607	. ,ycoefm0,ycoefh0)
608	c call dump2d(iim,jjm-1,ycoefm(2:klon-1,2), 'KZ ')
609	c call dump2d(iim,jjm-1,ycoefm0(2:klon-1,2),'KZMIN ')
610
611	if ( 1.eq.1 ) then
612	DO k = 1, klev
613	DO i = 1, knon
614	ycoefm(i,k) = MAX(ycoefm(i,k),ycoefm0(i,k))
615	ycoefh(i,k) = MAX(ycoefh(i,k),ycoefh0(i,k))
616	ENDDO
617	ENDDO
618	endif
619	cIM cf FH: 201103 END
620	endif !ok_kzmin
621	cIM: 261103
622
623
624	IF (iflag_pbl.ge.3) then
625
626	cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
627	c MELLOR ET YAMADA adapte a Mars Richard Fournier et Frederic Hourdin
628	cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
629
630	yzlay(1:knon,1)=
631	. RDyt(1:knon,1)/(0.5(ypaprs(1:knon,1)+ypplay(1:knon,1)))
632	. *(ypaprs(1:knon,1)-ypplay(1:knon,1))/RG
633	do k=2,klev
634	yzlay(1:knon,k)=
635	. yzlay(1:knon,k-1)+RD0.5(yt(1:knon,k-1)+yt(1:knon,k))
636	. /ypaprs(1:knon,k)*(ypplay(1:knon,k-1)-ypplay(1:knon,k))/RG
637	enddo
638	do k=1,klev
639	yteta(1:knon,k)=
640	. yt(1:knon,k)(ypaprs(1:knon,1)/ypplay(1:knon,k))*rkappa
641	. (1.+0.61yq(1:knon,k))
642	enddo
643	yzlev(1:knon,1)=0.
644	yzlev(1:knon,klev+1)=2.*yzlay(1:knon,klev)-yzlay(1:knon,klev-1)
645	do k=2,klev
646	yzlev(1:knon,k)=0.5*(yzlay(1:knon,k)+yzlay(1:knon,k-1))
647	enddo
648	DO k = 1, klev+1
649	DO j = 1, knon
650	i = ni(j)
651	yq2(j,k)=q2(i,k,nsrf)
652	enddo
653	enddo
654
655
656	c Bug introduit volontairement pour converger avec les resultats
657	c du papier sur les thermiques.
658	if (1.eq.1) then
659	y_cd_m(1:knon) = ycoefm(1:knon,1)
660	y_cd_h(1:knon) = ycoefh(1:knon,1)
661	else
662	y_cd_h(1:knon) = ycoefm(1:knon,1)
663	y_cd_m(1:knon) = ycoefh(1:knon,1)
664	endif
665	call ustarhb(knon,yu,yv,y_cd_m, yustar)
666
667	if (prt_level > 9) THEN
668	WRITE(lunout,*)'USTAR = ',yustar
669	ENDIF
670
671	c iflag_pbl peut etre utilise comme longuer de melange
672
673	if (iflag_pbl.ge.11) then
674	call vdif_kcay(knon,dtime,rg,rd,ypaprs,yt
675	s ,yzlev,yzlay,yu,yv,yteta
676	s ,y_cd_m,yq2,q2diag,ykmm,ykmn,yustar,
677	s iflag_pbl)
678	else
679	call yamada4(knon,dtime,rg,rd,ypaprs,yt
680	s ,yzlev,yzlay,yu,yv,yteta
681	s ,y_cd_m,yq2,ykmm,ykmn,ykmq,yustar,
682	s iflag_pbl)
683	endif
684
685	ycoefm(1:knon,1)=y_cd_m(1:knon)
686	ycoefh(1:knon,1)=y_cd_h(1:knon)
687	ycoefm(1:knon,2:klev)=ykmm(1:knon,2:klev)
688	ycoefh(1:knon,2:klev)=ykmn(1:knon,2:klev)
689
690
691	ENDIF
692
693	cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
694	c calculer la diffusion des vitesses "u" et "v"
695	cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
696
697	CALL clvent(knon,dtime,yu1,yv1,ycoefm,yt,yu,ypaprs,ypplay,ydelp,
698	s y_d_u,y_flux_u)
699	CALL clvent(knon,dtime,yu1,yv1,ycoefm,yt,yv,ypaprs,ypplay,ydelp,
700	s y_d_v,y_flux_v)
701
702	c pour le couplage
703	ytaux = y_flux_u(:,1)
704	ytauy = y_flux_v(:,1)
705
706	c FH modif sur le cdrag temperature
707	c$$$PB : déplace dans clcdrag
708	c$$$ do i=1,knon
709	c$$$ ycoefh(i,1)=ycoefm(i,1)*0.8
710	c$$$ enddo
711
712	cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
713	c calculer la diffusion de "q" et de "h"
714	cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
715	CALL clqh(dtime, itap, date0,jour, debut,lafin,
716	e rlon, rlat, cufi, cvfi,
717	e knon, nsrf, ni, pctsrf,
718	e soil_model, ytsoil,yqsol,
719	e ok_veget, ocean, npas, nexca,
720	e rmu0, co2_ppm, yrugos, yrugoro,
721	e yu1, yv1, ycoefh,
722	e yt,yq,yts,ypaprs,ypplay,
723	e ydelp,yrads,yalb, yalblw, ysnow, yqsurf,
724	e yrain_f, ysnow_f, yfder, ytaux, ytauy,
725	c -- LOOP
726	e ywindsp,
727	c -- LOOP
728	c$$$ e ysollw, ysolsw,
729	e ysollw, ysollwdown, ysolsw,yfluxlat,
730	s pctsrf_new, yagesno,
731	s y_d_t, y_d_q, y_d_ts, yz0_new,
732	s y_flux_t, y_flux_q, y_dflux_t, y_dflux_q,
733	s y_fqcalving,y_ffonte,y_run_off_lic_0,
734	cIM "slab" ocean
735	s y_flux_o, y_flux_g, ytslab, y_seaice)
736	c
737	c calculer la longueur de rugosite sur ocean
738	yrugm=0.
739	IF (nsrf.EQ.is_oce) THEN
740	DO j = 1, knon
741	yrugm(j) = 0.018ycoefm(j,1) (yu1(j)2+yv1(j)2)/RG
742	$ + 0.1114e-6 / sqrt(ycoefm(j,1) (yu1(j)2+yv1(j)2))
743	yrugm(j) = MAX(1.5e-05,yrugm(j))
744	ENDDO
745	ENDIF
746	DO j = 1, knon
747	y_dflux_t(j) = y_dflux_t(j) * ypct(j)
748	y_dflux_q(j) = y_dflux_q(j) * ypct(j)
749	yu1(j) = yu1(j) * ypct(j)
750	yv1(j) = yv1(j) * ypct(j)
751	ENDDO
752	c
753	DO k = 1, klev
754	DO j = 1, knon
755	i = ni(j)
756	ycoefh(j,k) = ycoefh(j,k) * ypct(j)
757	ycoefm(j,k) = ycoefm(j,k) * ypct(j)
758	y_d_t(j,k) = y_d_t(j,k) * ypct(j)
759	y_d_q(j,k) = y_d_q(j,k) * ypct(j)
760	C§§§ PB
761	flux_t(i,k,nsrf) = y_flux_t(j,k)
762	flux_q(i,k,nsrf) = y_flux_q(j,k)
763	flux_u(i,k,nsrf) = y_flux_u(j,k)
764	flux_v(i,k,nsrf) = y_flux_v(j,k)
765	c$$$ PB y_flux_t(j,k) = y_flux_t(j,k) * ypct(j)
766	c$$$ PB y_flux_q(j,k) = y_flux_q(j,k) * ypct(j)
767	y_d_u(j,k) = y_d_u(j,k) * ypct(j)
768	y_d_v(j,k) = y_d_v(j,k) * ypct(j)
769	c$$$ PB y_flux_u(j,k) = y_flux_u(j,k) * ypct(j)
770	c$$$ PB y_flux_v(j,k) = y_flux_v(j,k) * ypct(j)
771	ENDDO
772	ENDDO
773
774
775	evap(:,nsrf) = - flux_q(:,1,nsrf)
776	c
777	albe(:, nsrf) = 0.
778	alblw(:, nsrf) = 0.
779	snow(:, nsrf) = 0.
780	qsurf(:, nsrf) = 0.
781	rugos(:, nsrf) = 0.
782	fluxlat(:,nsrf) = 0.
783	DO j = 1, knon
784	i = ni(j)
785	d_ts(i,nsrf) = y_d_ts(j)
786	albe(i,nsrf) = yalb(j)
787	alblw(i,nsrf) = yalblw(j)
788	snow(i,nsrf) = ysnow(j)
789	qsurf(i,nsrf) = yqsurf(j)
790	rugos(i,nsrf) = yz0_new(j)
791	fluxlat(i,nsrf) = yfluxlat(j)
792	c$$$ pb rugmer(i) = yrugm(j)
793	IF (nsrf .EQ. is_oce) then
794	rugmer(i) = yrugm(j)
795	rugos(i,nsrf) = yrugm(j)
796	endif
797	cIM cf JLD ??
798	agesno(i,nsrf) = yagesno(j)
799	fqcalving(i,nsrf) = y_fqcalving(j)
800	ffonte(i,nsrf) = y_ffonte(j)
801	cdragh(i) = cdragh(i) + ycoefh(j,1)
802	cdragm(i) = cdragm(i) + ycoefm(j,1)
803	dflux_t(i) = dflux_t(i) + y_dflux_t(j)
804	dflux_q(i) = dflux_q(i) + y_dflux_q(j)
805	zu1(i) = zu1(i) + yu1(j)
806	zv1(i) = zv1(i) + yv1(j)
807	END DO
808	IF ( nsrf .eq. is_ter ) THEN
809	DO j = 1, knon
810	i = ni(j)
811	qsol(i) = yqsol(j)
812	END DO
813	END IF
814	IF ( nsrf .eq. is_lic ) THEN
815	DO j = 1, knon
816	i = ni(j)
817	run_off_lic_0(i) = y_run_off_lic_0(j)
818	END DO
819	END IF
820	c$$$ PB ajout pour soil
821	ftsoil(:,:,nsrf) = 0.
822	DO k = 1, nsoilmx
823	DO j = 1, knon
824	i = ni(j)
825	ftsoil(i, k, nsrf) = ytsoil(j,k)
826	END DO
827	END DO
828	c
829	#ifdef CRAY
830	DO k = 1, klev
831	DO j = 1, knon
832	i = ni(j)
833	#else
834	DO j = 1, knon
835	i = ni(j)
836	DO k = 1, klev
837	#endif
838	d_t(i,k) = d_t(i,k) + y_d_t(j,k)
839	d_q(i,k) = d_q(i,k) + y_d_q(j,k)
840	c$$$ PB flux_t(i,k) = flux_t(i,k) + y_flux_t(j,k)
841	c$$$ flux_q(i,k) = flux_q(i,k) + y_flux_q(j,k)
842	d_u(i,k) = d_u(i,k) + y_d_u(j,k)
843	d_v(i,k) = d_v(i,k) + y_d_v(j,k)
844	c$$$ PB flux_u(i,k) = flux_u(i,k) + y_flux_u(j,k)
845	c$$$ flux_v(i,k) = flux_v(i,k) + y_flux_v(j,k)
846	zcoefh(i,k) = zcoefh(i,k) + ycoefh(j,k)
847	ENDDO
848	ENDDO
849	c
850	c
851	#undef T2m
852	#define T2m
853	#ifdef T2m
854	ccc diagnostic t,q a 2m et u, v a 10m
855	c
856	DO j=1, knon
857	i = ni(j)
858	uzon(j) = yu(j,1) + y_d_u(j,1)
859	vmer(j) = yv(j,1) + y_d_v(j,1)
860	tair1(j) = yt(j,1) + y_d_t(j,1)
861	qair1(j) = yq(j,1) + y_d_q(j,1)
862	zgeo1(j) = RD * tair1(j) / (0.5*(ypaprs(j,1)+ypplay(j,1)))
863	& * (ypaprs(j,1)-ypplay(j,1))
864	tairsol(j) = yts(j) + y_d_ts(j)
865	rugo1(j) = yrugos(j)
866	IF(nsrf.EQ.is_oce) THEN
867	rugo1(j) = rugos(i,nsrf)
868	ENDIF
869	psfce(j)=ypaprs(j,1)
870	patm(j)=ypplay(j,1)
871	c
872	qairsol(j) = yqsurf(j)
873	ENDDO
874	c
875	CALL stdlevvar(klon, knon, nsrf, zxli,
876	& uzon, vmer, tair1, qair1, zgeo1,
877	& tairsol, qairsol, rugo1, psfce, patm,
878	cIM & yt2m, yq2m, yu10m)
879	& yt2m, yq2m, yt10m, yq10m, yu10m, yustar)
880	cIM 081204 END
881	c
882	c
883	DO j=1, knon
884	i = ni(j)
885	t2m(i,nsrf)=yt2m(j)
886
887	c
888	q2m(i,nsrf)=yq2m(j)
889	c
890	c u10m, v10m : composantes du vent a 10m sans spirale de Ekman
891	u10m(i,nsrf)=(yu10m(j) * uzon(j))/sqrt(uzon(j)2+vmer(j)2)
892	v10m(i,nsrf)=(yu10m(j) * vmer(j))/sqrt(uzon(j)2+vmer(j)2)
893	c
894	ENDDO
895	c
896	cIM cf AM : pbl, HBTM
897	DO i = 1, knon
898	y_cd_h(i) = ycoefh(i,1)
899	y_cd_m(i) = ycoefm(i,1)
900	ENDDO
901	c print*,'appel hbtm2'
902	CALL HBTM(knon, ypaprs, ypplay,
903	. yt2m,yt10m,yq2m,yq10m,yustar,
904	. y_flux_t,y_flux_q,yu,yv,yt,yq,
905	. ypblh,ycapCL,yoliqCL,ycteiCL,ypblT,
906	. ytherm,ytrmb1,ytrmb2,ytrmb3,ylcl)
907	c print*,'fin hbtm2'
908	c
909	DO j=1, knon
910	i = ni(j)
911	pblh(i,nsrf) = ypblh(j)
912	plcl(i,nsrf) = ylcl(j)
913	capCL(i,nsrf) = ycapCL(j)
914	oliqCL(i,nsrf) = yoliqCL(j)
915	cteiCL(i,nsrf) = ycteiCL(j)
916	pblT(i,nsrf) = ypblT(j)
917	therm(i,nsrf) = ytherm(j)
918	trmb1(i,nsrf) = ytrmb1(j)
919	trmb2(i,nsrf) = ytrmb2(j)
920	trmb3(i,nsrf) = ytrmb3(j)
921	ENDDO
922	c
923	#else
924	cIM cf AM : pbl, HBTM
925	DO j=1, knon
926	i = ni(j)
927	pblh(i,nsrf) = 0.
928	plcl(i,nsrf) = 0.
929	capCL(i,nsrf) = 0.
930	oliqCL(i,nsrf) = 0.
931	cteiCL(i,nsrf) = 0.
932	pblT(i,nsrf) = 0.
933	therm(i,nsrf) = 0.
934	trmb1(i,nsrf) = 0.
935	trmb2(i,nsrf) = 0.
936	trmb3(i,nsrf) = 0.
937	ENDDO
938	DO j=1, knon
939	i = ni(j)
940	t2m(i,nsrf)=0.
941	q2m(i,nsrf)=0.
942	u10m(i,nsrf)=0.
943	v10m(i,nsrf)=0.
944	ENDDO
945	#endif
946
947	do j=1,knon
948	do k=1,klev+1
949	i=ni(j)
950	q2(i,k,nsrf)=yq2(j,k)
951	enddo
952	enddo
953	cIM "slab" ocean
954	IF(OCEAN.EQ.'slab '.OR.OCEAN.EQ.'force ') THEN
955	IF (nsrf.EQ.is_oce) THEN
956	DO j = 1, knon
957	c on projette sur la grille globale
958	i = ni(j)
959	IF(pctsrf_new(i,is_oce).GT.epsfra) THEN
960	flux_o(i) = y_flux_o(j)
961	ELSE
962	flux_o(i) = 0.
963	ENDIF
964	ENDDO
965	ENDIF
966	c
967	IF (nsrf.EQ.is_sic) THEN
968	DO j = 1, knon
969	i = ni(j)
970	cIM 230604 on pondere lorsque l'on fait le bilan au sol : flux_g(i) = y_flux_g(j)*ypct(j)
971	IF(pctsrf_new(i,is_sic).GT.epsfra) THEN
972	flux_g(i) = y_flux_g(j)
973	ELSE
974	flux_g(i) = 0.
975	ENDIF
976	ENDDO
977	ENDIF !nsrf.EQ.is_sic
978	ENDIF !OCEAN
979	c
980	IF(OCEAN.EQ.'slab ') THEN
981	cIM appel interfoce_slab sur is_oce IF(nsrf.EQ.is_oce) then
982	cIM appel interfoce_slab sur is_sic
983	IF(nsrf.EQ.is_sic) then !appel sic
984	tslab(1:klon) = ytslab(1:klon)
985	seaice(1:klon) = y_seaice(1:klon)
986	ENDIF !nsrf
987	ENDIF !OCEAN
988	99999 CONTINUE
989	C
990	C On utilise les nouvelles surfaces
991	C A rajouter: conservation de l'albedo
992	C
993	rugos(:,is_oce) = rugmer
994	pctsrf = pctsrf_new
995
996	RETURN
997	END
998	SUBROUTINE clqh(dtime,itime, date0,jour,debut,lafin,
999	e rlon, rlat, cufi, cvfi,
1000	e knon, nisurf, knindex, pctsrf,
1001	$ soil_model,tsoil,qsol,
1002	e ok_veget, ocean, npas, nexca,
1003	e rmu0, co2_ppm, rugos, rugoro,
1004	e u1lay,v1lay,coef,
1005	e t,q,ts,paprs,pplay,
1006	e delp,radsol,albedo,alblw,snow,qsurf,
1007	e precip_rain, precip_snow, fder, taux, tauy,
1008	c -- LOOP
1009	e ywindsp,
1010	c -- LOOP
1011	$ sollw, sollwdown, swnet,fluxlat,
1012	s pctsrf_new, agesno,
1013	s d_t, d_q, d_ts, z0_new,
1014	s flux_t, flux_q,dflux_s,dflux_l,
1015	s fqcalving,ffonte,run_off_lic_0,
1016	cIM "slab" ocean
1017	s flux_o,flux_g,tslab,seaice)
1018
1019	USE interface_surf
1020
1021	IMPLICIT none
1022	c======================================================================
1023	c Auteur(s): Z.X. Li (LMD/CNRS) date: 19930818
1024	c Objet: diffusion verticale de "q" et de "h"
1025	c======================================================================
1026	#include "dimensions.h"
1027	#include "dimphy.h"
1028	#include "YOMCST.h"
1029	#include "YOETHF.h"
1030	#include "FCTTRE.h"
1031	#include "indicesol.h"
1032	#include "dimsoil.h"
1033	#include "iniprint.h"
1034
1035	c Arguments:
1036	INTEGER knon
1037	REAL dtime ! intervalle du temps (s)
1038	real date0
1039	REAL u1lay(klon) ! vitesse u de la 1ere couche (m/s)
1040	REAL v1lay(klon) ! vitesse v de la 1ere couche (m/s)
1041	REAL coef(klon,klev) ! le coefficient d'echange (m**2/s)
1042	c multiplie par le cisaillement du
1043	c vent (dV/dz); la premiere valeur
1044	c indique la valeur de Cdrag (sans unite)
1045	REAL t(klon,klev) ! temperature (K)
1046	REAL q(klon,klev) ! humidite specifique (kg/kg)
1047	REAL ts(klon) ! temperature du sol (K)
1048	REAL evap(klon) ! evaporation au sol
1049	REAL paprs(klon,klev+1) ! pression a inter-couche (Pa)
1050	REAL pplay(klon,klev) ! pression au milieu de couche (Pa)
1051	REAL delp(klon,klev) ! epaisseur de couche en pression (Pa)
1052	REAL radsol(klon) ! ray. net au sol (Solaire+IR) W/m2
1053	REAL albedo(klon) ! albedo de la surface
1054	REAL alblw(klon)
1055	REAL snow(klon) ! hauteur de neige
1056	REAL qsurf(klon) ! humidite de l'air au dessus de la surface
1057	real precip_rain(klon), precip_snow(klon)
1058	REAL agesno(klon)
1059	REAL rugoro(klon)
1060	REAL run_off_lic_0(klon)! runof glacier au pas de temps precedent
1061	integer jour ! jour de l'annee en cours
1062	real rmu0(klon) ! cosinus de l'angle solaire zenithal
1063	real rugos(klon) ! rugosite
1064	integer knindex(klon)
1065	real pctsrf(klon,nbsrf)
1066	real rlon(klon), rlat(klon), cufi(klon), cvfi(klon)
1067	logical ok_veget
1068	REAL co2_ppm ! taux CO2 atmosphere
1069	character*6 ocean
1070	integer npas, nexca
1071	c -- LOOP
1072	REAL yu10mx(klon)
1073	REAL yu10my(klon)
1074	REAL ywindsp(klon)
1075	c -- LOOP
1076
1077
1078	c
1079	REAL d_t(klon,klev) ! incrementation de "t"
1080	REAL d_q(klon,klev) ! incrementation de "q"
1081	REAL d_ts(klon) ! incrementation de "ts"
1082	REAL flux_t(klon,klev) ! (diagnostic) flux de la chaleur
1083	c sensible, flux de Cp*T, positif vers
1084	c le bas: j/(m**2 s) c.a.d.: W/m2
1085	REAL flux_q(klon,klev) ! flux de la vapeur d'eau:kg/(m**2 s)
1086	REAL dflux_s(klon) ! derivee du flux sensible dF/dTs
1087	REAL dflux_l(klon) ! derivee du flux latent dF/dTs
1088	cIM cf JLD
1089	c Flux thermique utiliser pour fondre la neige
1090	REAL ffonte(klon)
1091	c Flux d'eau "perdue" par la surface et nécessaire pour que limiter la
1092	c hauteur de neige, en kg/m2/s
1093	REAL fqcalving(klon)
1094	cIM "slab" ocean
1095	REAL tslab(klon) !temperature du slab ocean (K) (OCEAN='slab ')
1096	REAL seaice(klon) ! glace de mer en kg/m2
1097	REAL flux_o(klon) ! flux entre l'ocean et l'atmosphere W/m2
1098	REAL flux_g(klon) ! flux entre l'ocean et la glace de mer W/m2
1099	c
1100	c======================================================================
1101	REAL t_grnd ! temperature de rappel pour glace de mer
1102	PARAMETER (t_grnd=271.35)
1103	REAL t_coup
1104	PARAMETER(t_coup=273.15)
1105	c======================================================================
1106	INTEGER i, k
1107	REAL zx_cq(klon,klev)
1108	REAL zx_dq(klon,klev)
1109	REAL zx_ch(klon,klev)
1110	REAL zx_dh(klon,klev)
1111	REAL zx_buf1(klon)
1112	REAL zx_buf2(klon)
1113	REAL zx_coef(klon,klev)
1114	REAL local_h(klon,klev) ! enthalpie potentielle
1115	REAL local_q(klon,klev)
1116	REAL local_ts(klon)
1117	REAL psref(klon) ! pression de reference pour temperature potent.
1118	REAL zx_pkh(klon,klev), zx_pkf(klon,klev)
1119	c======================================================================
1120	c contre-gradient pour la vapeur d'eau: (kg/kg)/metre
1121	REAL gamq(klon,2:klev)
1122	c contre-gradient pour la chaleur sensible: Kelvin/metre
1123	REAL gamt(klon,2:klev)
1124	REAL z_gamaq(klon,2:klev), z_gamah(klon,2:klev)
1125	REAL zdelz
1126	c======================================================================
1127	#include "compbl.h"
1128	c======================================================================
1129	c Rajout pour l'interface
1130	integer itime
1131	integer nisurf
1132	logical debut, lafin
1133	real zlev1(klon)
1134	real fder(klon), taux(klon), tauy(klon)
1135	real temp_air(klon), spechum(klon)
1136	real epot_air(klon), ccanopy(klon)
1137	real tq_cdrag(klon), petAcoef(klon), peqAcoef(klon)
1138	real petBcoef(klon), peqBcoef(klon)
1139	real sollw(klon), sollwdown(klon), swnet(klon), swdown(klon)
1140	real p1lay(klon)
1141	c$$$C PB ajout pour soil
1142	LOGICAL soil_model
1143	REAL tsoil(klon, nsoilmx)
1144	REAL qsol(klon)
1145
1146	! Parametres de sortie
1147	real fluxsens(klon), fluxlat(klon)
1148	real tsol_rad(klon), tsurf_new(klon), alb_new(klon)
1149	real emis_new(klon), z0_new(klon)
1150	real pctsrf_new(klon,nbsrf)
1151	c JLD
1152	real zzpk
1153	C
1154	character (len = 20) :: modname = 'Debut clqh'
1155	LOGICAL check
1156	PARAMETER (check=.false.)
1157	C
1158	if (check) THEN
1159	write(,) modname,' nisurf=',nisurf
1160	call flush(6)
1161	endif
1162	c
1163	if (check) THEN
1164	WRITE(,)' qsurf (min, max)'
1165	$ , minval(qsurf(1:knon)), maxval(qsurf(1:knon))
1166	call flush(6)
1167	ENDIF
1168	C
1169	C
1170	if (iflag_pbl.eq.1) then
1171	do k = 3, klev
1172	do i = 1, knon
1173	gamq(i,k)= 0.0
1174	gamt(i,k)= -1.0e-03
1175	enddo
1176	enddo
1177	do i = 1, knon
1178	gamq(i,2) = 0.0
1179	gamt(i,2) = -2.5e-03
1180	enddo
1181	else
1182	do k = 2, klev
1183	do i = 1, knon
1184	gamq(i,k) = 0.0
1185	gamt(i,k) = 0.0
1186	enddo
1187	enddo
1188	endif
1189
1190	DO i = 1, knon
1191	psref(i) = paprs(i,1) !pression de reference est celle au sol
1192	local_ts(i) = ts(i)
1193	ENDDO
1194	DO k = 1, klev
1195	DO i = 1, knon
1196	zx_pkh(i,k) = (psref(i)/paprs(i,k))**RKAPPA
1197	zx_pkf(i,k) = (psref(i)/pplay(i,k))**RKAPPA
1198	local_h(i,k) = RCPD * t(i,k) * zx_pkf(i,k)
1199	local_q(i,k) = q(i,k)
1200	ENDDO
1201	ENDDO
1202	c
1203	c Convertir les coefficients en variables convenables au calcul:
1204	c
1205	c
1206	DO k = 2, klev
1207	DO i = 1, knon
1208	zx_coef(i,k) = coef(i,k)*RG/(pplay(i,k-1)-pplay(i,k))
1209	. (paprs(i,k)2/(t(i,k)+t(i,k-1))/RD)**2
1210	zx_coef(i,k) = zx_coef(i,k) * dtime*RG
1211	ENDDO
1212	ENDDO
1213	c
1214	c Preparer les flux lies aux contre-gardients
1215	c
1216	DO k = 2, klev
1217	DO i = 1, knon
1218	zdelz = RD * (t(i,k-1)+t(i,k))/2.0 / RG /paprs(i,k)
1219	. *(pplay(i,k-1)-pplay(i,k))
1220	z_gamaq(i,k) = gamq(i,k) * zdelz
1221	z_gamah(i,k) = gamt(i,k) * zdelz RCPD zx_pkh(i,k)
1222	ENDDO
1223	ENDDO
1224	DO i = 1, knon
1225	zx_buf1(i) = zx_coef(i,klev) + delp(i,klev)
1226	zx_cq(i,klev) = (local_q(i,klev)*delp(i,klev)
1227	. -zx_coef(i,klev)*z_gamaq(i,klev))/zx_buf1(i)
1228	zx_dq(i,klev) = zx_coef(i,klev) / zx_buf1(i)
1229	c
1230	zzpk=(pplay(i,klev)/psref(i))**RKAPPA
1231	zx_buf2(i) = zzpk*delp(i,klev) + zx_coef(i,klev)
1232	zx_ch(i,klev) = (local_h(i,klev)zzpkdelp(i,klev)
1233	. -zx_coef(i,klev)*z_gamah(i,klev))/zx_buf2(i)
1234	zx_dh(i,klev) = zx_coef(i,klev) / zx_buf2(i)
1235	ENDDO
1236	DO k = klev-1, 2 , -1
1237	DO i = 1, knon
1238	zx_buf1(i) = delp(i,k)+zx_coef(i,k)
1239	. +zx_coef(i,k+1)*(1.-zx_dq(i,k+1))
1240	zx_cq(i,k) = (local_q(i,k)*delp(i,k)
1241	. +zx_coef(i,k+1)*zx_cq(i,k+1)
1242	. +zx_coef(i,k+1)*z_gamaq(i,k+1)
1243	. -zx_coef(i,k)*z_gamaq(i,k))/zx_buf1(i)
1244	zx_dq(i,k) = zx_coef(i,k) / zx_buf1(i)
1245	c
1246	zzpk=(pplay(i,k)/psref(i))**RKAPPA
1247	zx_buf2(i) = zzpk*delp(i,k)+zx_coef(i,k)
1248	. +zx_coef(i,k+1)*(1.-zx_dh(i,k+1))
1249	zx_ch(i,k) = (local_h(i,k)zzpkdelp(i,k)
1250	. +zx_coef(i,k+1)*zx_ch(i,k+1)
1251	. +zx_coef(i,k+1)*z_gamah(i,k+1)
1252	. -zx_coef(i,k)*z_gamah(i,k))/zx_buf2(i)
1253	zx_dh(i,k) = zx_coef(i,k) / zx_buf2(i)
1254	ENDDO
1255	ENDDO
1256	C
1257	C nouvelle formulation JL Dufresne
1258	C
1259	C q1 = zx_cq(i,1) + zx_dq(i,1) * Flux_Q(i,1) * dt
1260	C h1 = zx_ch(i,1) + zx_dh(i,1) * Flux_H(i,1) * dt
1261	C
1262	DO i = 1, knon
1263	zx_buf1(i) = delp(i,1) + zx_coef(i,2)*(1.-zx_dq(i,2))
1264	zx_cq(i,1) = (local_q(i,1)*delp(i,1)
1265	. +zx_coef(i,2)*(z_gamaq(i,2)+zx_cq(i,2)))
1266	. /zx_buf1(i)
1267	zx_dq(i,1) = -1. * RG / zx_buf1(i)
1268	c
1269	zzpk=(pplay(i,1)/psref(i))**RKAPPA
1270	zx_buf2(i) = zzpkdelp(i,1) + zx_coef(i,2)(1.-zx_dh(i,2))
1271	zx_ch(i,1) = (local_h(i,1)zzpkdelp(i,1)
1272	. +zx_coef(i,2)*(z_gamah(i,2)+zx_ch(i,2)))
1273	. /zx_buf2(i)
1274	zx_dh(i,1) = -1. * RG / zx_buf2(i)
1275	ENDDO
1276
1277	C Appel a interfsurf (appel generique) routine d'interface avec la surface
1278
1279	c initialisation
1280	petAcoef =0.
1281	peqAcoef = 0.
1282	petBcoef =0.
1283	peqBcoef = 0.
1284	p1lay =0.
1285
1286	c do i = 1, knon
1287	petAcoef(1:knon) = zx_ch(1:knon,1)
1288	peqAcoef(1:knon) = zx_cq(1:knon,1)
1289	petBcoef(1:knon) = zx_dh(1:knon,1)
1290	peqBcoef(1:knon) = zx_dq(1:knon,1)
1291	tq_cdrag(1:knon) =coef(1:knon,1)
1292	temp_air(1:knon) =t(1:knon,1)
1293	epot_air(1:knon) =local_h(1:knon,1)
1294	spechum(1:knon)=q(1:knon,1)
1295	p1lay(1:knon) = pplay(1:knon,1)
1296	zlev1(1:knon) = delp(1:knon,1)
1297	c swnet = swdown * (1. - albedo)
1298	c
1299	cIM swdown=flux SW incident sur terres
1300	cIM swdown=flux SW net sur les autres surfaces
1301	cIM swdown(1:knon) = swnet(1:knon)
1302	if(nisurf.eq.is_ter) THEN
1303	swdown(1:knon) = swnet(1:knon)/(1-albedo(1:knon))
1304	else
1305	swdown(1:knon) = swnet(1:knon)
1306	endif
1307	c enddo
1308	ccanopy = co2_ppm
1309
1310	CALL interfsurf(itime, dtime, date0, jour, rmu0,
1311	e klon, iim, jjm, nisurf, knon, knindex, pctsrf,
1312	e rlon, rlat, cufi, cvfi,
1313	e debut, lafin, ok_veget, soil_model, nsoilmx,tsoil, qsol,
1314	e zlev1, u1lay, v1lay, temp_air, spechum, epot_air, ccanopy,
1315	e tq_cdrag, petAcoef, peqAcoef, petBcoef, peqBcoef,
1316	e precip_rain, precip_snow, sollw, sollwdown, swnet, swdown,
1317	e fder, taux, tauy,
1318	c -- LOOP
1319	e ywindsp,
1320	c -- LOOP
1321	e rugos, rugoro,
1322	e albedo, snow, qsurf,
1323	e ts, p1lay, psref, radsol,
1324	e ocean, npas, nexca, zmasq,
1325	s evap, fluxsens, fluxlat, dflux_l, dflux_s,
1326	s tsol_rad, tsurf_new, alb_new, alblw, emis_new, z0_new,
1327	s pctsrf_new, agesno,fqcalving,ffonte, run_off_lic_0,
1328	cIM "slab" ocean
1329	s flux_o, flux_g, tslab, seaice)
1330
1331
1332	do i = 1, knon
1333	flux_t(i,1) = fluxsens(i)
1334	flux_q(i,1) = - evap(i)
1335	d_ts(i) = tsurf_new(i) - ts(i)
1336	albedo(i) = alb_new(i)
1337	enddo
1338
1339	c==== une fois on a zx_h_ts, on peut faire l'iteration ========
1340	DO i = 1, knon
1341	local_h(i,1) = zx_ch(i,1) + zx_dh(i,1)flux_t(i,1)dtime
1342	local_q(i,1) = zx_cq(i,1) + zx_dq(i,1)flux_q(i,1)dtime
1343	ENDDO
1344	DO k = 2, klev
1345	DO i = 1, knon
1346	local_q(i,k) = zx_cq(i,k) + zx_dq(i,k)*local_q(i,k-1)
1347	local_h(i,k) = zx_ch(i,k) + zx_dh(i,k)*local_h(i,k-1)
1348	ENDDO
1349	ENDDO
1350	c======================================================================
1351	c== flux_q est le flux de vapeur d'eau: kg/(m**2 s) positive vers bas
1352	c== flux_t est le flux de cpt (energie sensible): j/(m**2 s)
1353	DO k = 2, klev
1354	DO i = 1, knon
1355	flux_q(i,k) = (zx_coef(i,k)/RG/dtime)
1356	. * (local_q(i,k)-local_q(i,k-1)+z_gamaq(i,k))
1357	flux_t(i,k) = (zx_coef(i,k)/RG/dtime)
1358	. * (local_h(i,k)-local_h(i,k-1)+z_gamah(i,k))
1359	. / zx_pkh(i,k)
1360	ENDDO
1361	ENDDO
1362	c======================================================================
1363	C Calcul tendances
1364	DO k = 1, klev
1365	DO i = 1, knon
1366	d_t(i,k) = local_h(i,k)/zx_pkf(i,k)/RCPD - t(i,k)
1367	d_q(i,k) = local_q(i,k) - q(i,k)
1368	ENDDO
1369	ENDDO
1370	c
1371
1372	RETURN
1373	END
1374	SUBROUTINE clvent(knon,dtime, u1lay,v1lay,coef,t,ven,
1375	e paprs,pplay,delp,
1376	s d_ven,flux_v)
1377	IMPLICIT none
1378	c======================================================================
1379	c Auteur(s): Z.X. Li (LMD/CNRS) date: 19930818
1380	c Objet: diffusion vertical de la vitesse "ven"
1381	c======================================================================
1382	c Arguments:
1383	c dtime----input-R- intervalle du temps (en second)
1384	c u1lay----input-R- vent u de la premiere couche (m/s)
1385	c v1lay----input-R- vent v de la premiere couche (m/s)
1386	c coef-----input-R- le coefficient d'echange (m**2/s) multiplie par
1387	c le cisaillement du vent (dV/dz); la premiere
1388	c valeur indique la valeur de Cdrag (sans unite)
1389	c t--------input-R- temperature (K)
1390	c ven------input-R- vitesse horizontale (m/s)
1391	c paprs----input-R- pression a inter-couche (Pa)
1392	c pplay----input-R- pression au milieu de couche (Pa)
1393	c delp-----input-R- epaisseur de couche (Pa)
1394	c
1395	c
1396	c d_ven----output-R- le changement de "ven"
1397	c flux_v---output-R- (diagnostic) flux du vent: (kg m/s)/(m**2 s)
1398	c======================================================================
1399	#include "dimensions.h"
1400	#include "dimphy.h"
1401	#include "iniprint.h"
1402	INTEGER knon
1403	REAL dtime
1404	REAL u1lay(klon), v1lay(klon)
1405	REAL coef(klon,klev)
1406	REAL t(klon,klev), ven(klon,klev)
1407	REAL paprs(klon,klev+1), pplay(klon,klev), delp(klon,klev)
1408	REAL d_ven(klon,klev)
1409	REAL flux_v(klon,klev)
1410	c======================================================================
1411	#include "YOMCST.h"
1412	c======================================================================
1413	INTEGER i, k
1414	REAL zx_cv(klon,2:klev)
1415	REAL zx_dv(klon,2:klev)
1416	REAL zx_buf(klon)
1417	REAL zx_coef(klon,klev)
1418	REAL local_ven(klon,klev)
1419	REAL zx_alf1(klon), zx_alf2(klon)
1420	c======================================================================
1421	DO k = 1, klev
1422	DO i = 1, knon
1423	local_ven(i,k) = ven(i,k)
1424	ENDDO
1425	ENDDO
1426	c======================================================================
1427	DO i = 1, knon
1428	ccc zx_alf1(i) = (paprs(i,1)-pplay(i,2))/(pplay(i,1)-pplay(i,2))
1429	zx_alf1(i) = 1.0
1430	zx_alf2(i) = 1.0 - zx_alf1(i)
1431	zx_coef(i,1) = coef(i,1)
1432	. * (1.0+SQRT(u1lay(i)2+v1lay(i)2))
1433	. * pplay(i,1)/(RD*t(i,1))
1434	zx_coef(i,1) = zx_coef(i,1) * dtime*RG
1435	ENDDO
1436	c======================================================================
1437	DO k = 2, klev
1438	DO i = 1, knon
1439	zx_coef(i,k) = coef(i,k)*RG/(pplay(i,k-1)-pplay(i,k))
1440	. (paprs(i,k)2/(t(i,k)+t(i,k-1))/RD)**2
1441	zx_coef(i,k) = zx_coef(i,k) * dtime*RG
1442	ENDDO
1443	ENDDO
1444	c======================================================================
1445	DO i = 1, knon
1446	zx_buf(i) = delp(i,1) + zx_coef(i,1)*zx_alf1(i)+zx_coef(i,2)
1447	zx_cv(i,2) = local_ven(i,1)*delp(i,1) / zx_buf(i)
1448	zx_dv(i,2) = (zx_coef(i,2)-zx_alf2(i)*zx_coef(i,1))
1449	. /zx_buf(i)
1450	ENDDO
1451	DO k = 3, klev
1452	DO i = 1, knon
1453	zx_buf(i) = delp(i,k-1) + zx_coef(i,k)
1454	. + zx_coef(i,k-1)*(1.-zx_dv(i,k-1))
1455	zx_cv(i,k) = (local_ven(i,k-1)*delp(i,k-1)
1456	. +zx_coef(i,k-1)*zx_cv(i,k-1) )/zx_buf(i)
1457	zx_dv(i,k) = zx_coef(i,k)/zx_buf(i)
1458	ENDDO
1459	ENDDO
1460	DO i = 1, knon
1461	local_ven(i,klev) = ( local_ven(i,klev)*delp(i,klev)
1462	. +zx_coef(i,klev)*zx_cv(i,klev) )
1463	. / ( delp(i,klev) + zx_coef(i,klev)
1464	. -zx_coef(i,klev)*zx_dv(i,klev) )
1465	ENDDO
1466	DO k = klev-1, 1, -1
1467	DO i = 1, knon
1468	local_ven(i,k) = zx_cv(i,k+1) + zx_dv(i,k+1)*local_ven(i,k+1)
1469	ENDDO
1470	ENDDO
1471	c======================================================================
1472	c== flux_v est le flux de moment angulaire (positif vers bas)
1473	c== dont l'unite est: (kg m/s)/(m**2 s)
1474	DO i = 1, knon
1475	flux_v(i,1) = zx_coef(i,1)/(RG*dtime)
1476	. (local_ven(i,1)zx_alf1(i)
1477	. +local_ven(i,2)*zx_alf2(i))
1478	ENDDO
1479	DO k = 2, klev
1480	DO i = 1, knon
1481	flux_v(i,k) = zx_coef(i,k)/(RG*dtime)
1482	. * (local_ven(i,k)-local_ven(i,k-1))
1483	ENDDO
1484	ENDDO
1485	c
1486	DO k = 1, klev
1487	DO i = 1, knon
1488	d_ven(i,k) = local_ven(i,k) - ven(i,k)
1489	ENDDO
1490	ENDDO
1491	c
1492	RETURN
1493	END
1494	SUBROUTINE coefkz(nsrf, knon, paprs, pplay,
1495	cIM 261103
1496	. ksta, ksta_ter,
1497	cIM 261103
1498	. ts, rugos,
1499	. u,v,t,q,
1500	. qsurf,
1501	. pcfm, pcfh)
1502	IMPLICIT none
1503	c======================================================================
1504	c Auteur(s) F. Hourdin, M. Forichon, Z.X. Li (LMD/CNRS) date: 19930922
1505	c (une version strictement identique a l'ancien modele)
1506	c Objet: calculer le coefficient du frottement du sol (Cdrag) et les
1507	c coefficients d'echange turbulent dans l'atmosphere.
1508	c Arguments:
1509	c nsrf-----input-I- indicateur de la nature du sol
1510	c knon-----input-I- nombre de points a traiter
1511	c paprs----input-R- pression a chaque intercouche (en Pa)
1512	c pplay----input-R- pression au milieu de chaque couche (en Pa)
1513	c ts-------input-R- temperature du sol (en Kelvin)
1514	c rugos----input-R- longeur de rugosite (en m)
1515	c u--------input-R- vitesse u
1516	c v--------input-R- vitesse v
1517	c t--------input-R- temperature (K)
1518	c q--------input-R- vapeur d'eau (kg/kg)
1519	c
1520	c itop-----output-I- numero de couche du sommet de la couche limite
1521	c pcfm-----output-R- coefficients a calculer (vitesse)
1522	c pcfh-----output-R- coefficients a calculer (chaleur et humidite)
1523	c======================================================================
1524	#include "dimensions.h"
1525	#include "dimphy.h"
1526	#include "YOMCST.h"
1527	#include "indicesol.h"
1528	#include "iniprint.h"
1529	c
1530	c Arguments:
1531	c
1532	INTEGER knon, nsrf
1533	REAL ts(klon)
1534	REAL paprs(klon,klev+1), pplay(klon,klev)
1535	REAL u(klon,klev), v(klon,klev), t(klon,klev), q(klon,klev)
1536	REAL rugos(klon)
1537	c
1538	REAL pcfm(klon,klev), pcfh(klon,klev)
1539	INTEGER itop(klon)
1540	c
1541	c Quelques constantes et options:
1542	c
1543	REAL cepdu2, ckap, cb, cc, cd, clam
1544	PARAMETER (cepdu2 =(0.1)**2)
1545	PARAMETER (CKAP=0.4)
1546	PARAMETER (cb=5.0)
1547	PARAMETER (cc=5.0)
1548	PARAMETER (cd=5.0)
1549	PARAMETER (clam=160.0)
1550	REAL ratqs ! largeur de distribution de vapeur d'eau
1551	PARAMETER (ratqs=0.05)
1552	LOGICAL richum ! utilise le nombre de Richardson humide
1553	PARAMETER (richum=.TRUE.)
1554	REAL ric ! nombre de Richardson critique
1555	PARAMETER(ric=0.4)
1556	REAL prandtl
1557	PARAMETER (prandtl=0.4)
1558	REAL kstable ! diffusion minimale (situation stable)
1559	! GKtest
1560	! PARAMETER (kstable=1.0e-10)
1561	REAL ksta, ksta_ter
1562	cIM: 261103 REAL kstable_ter, kstable_sinon
1563	cIM: 211003 cf GK PARAMETER (kstable_ter = 1.0e-6)
1564	cIM: 261103 PARAMETER (kstable_ter = 1.0e-8)
1565	cIM: 261103 PARAMETER (kstable_ter = 1.0e-10)
1566	cIM: 261103 PARAMETER (kstable_sinon = 1.0e-10)
1567	! fin GKtest
1568	REAL mixlen ! constante controlant longueur de melange
1569	PARAMETER (mixlen=35.0)
1570	INTEGER isommet ! le sommet de la couche limite
1571	PARAMETER (isommet=klev)
1572	LOGICAL tvirtu ! calculer Ri d'une maniere plus performante
1573	PARAMETER (tvirtu=.TRUE.)
1574	LOGICAL opt_ec ! formule du Centre Europeen dans l'atmosphere
1575	PARAMETER (opt_ec=.FALSE.)
1576
1577	#include "compbl.h"
1578	c
1579	c Variables locales:
1580	c
1581	INTEGER i, k, kk !IM 120704
1582	REAL zgeop(klon,klev)
1583	REAL zmgeom(klon)
1584	REAL zri(klon)
1585	REAL zl2(klon)
1586
1587	REAL u1(klon), v1(klon), t1(klon), q1(klon), z1(klon)
1588	REAL pcfm1(klon), pcfh1(klon)
1589	c
1590	REAL zdphi, zdu2, ztvd, ztvu, zcdn
1591	REAL zscf
1592	REAL zt, zq, zdelta, zcvm5, zcor, zqs, zfr, zdqs
1593	REAL z2geomf, zalh2, zalm2, zscfh, zscfm
1594	REAL t_coup
1595	PARAMETER (t_coup=273.15)
1596	cIM
1597	LOGICAL check
1598	PARAMETER (check=.false.)
1599	c
1600	c contre-gradient pour la chaleur sensible: Kelvin/metre
1601	REAL gamt(2:klev)
1602	real qsurf(klon)
1603	c
1604	LOGICAL appel1er
1605	SAVE appel1er
1606	c
1607	c Fonctions thermodynamiques et fonctions d'instabilite
1608	REAL fsta, fins, x
1609	LOGICAL zxli ! utiliser un jeu de fonctions simples
1610	PARAMETER (zxli=.FALSE.)
1611	c
1612	#include "YOETHF.h"
1613	#include "FCTTRE.h"
1614	fsta(x) = 1.0 / (1.0+10.0x(1+8.0*x))
1615	fins(x) = SQRT(1.0-18.0*x)
1616	c
1617	DATA appel1er /.TRUE./
1618	c
1619	IF (appel1er) THEN
1620	if (prt_level > 9) THEN
1621	WRITE(lunout,*)'coefkz, opt_ec:', opt_ec
1622	WRITE(lunout,*)'coefkz, richum:', richum
1623	IF (richum) WRITE(lunout,*)'coefkz, ratqs:', ratqs
1624	WRITE(lunout,*)'coefkz, isommet:', isommet
1625	WRITE(lunout,*)'coefkz, tvirtu:', tvirtu
1626	appel1er = .FALSE.
1627	endif
1628	ENDIF
1629	c
1630	c Initialiser les sorties
1631	c
1632	DO k = 1, klev
1633	DO i = 1, knon
1634	pcfm(i,k) = 0.0
1635	pcfh(i,k) = 0.0
1636	ENDDO
1637	ENDDO
1638	DO i = 1, knon
1639	itop(i) = 0
1640	ENDDO
1641
1642	c
1643	c Prescrire la valeur de contre-gradient
1644	c
1645	if (iflag_pbl.eq.1) then
1646	DO k = 3, klev
1647	gamt(k) = -1.0E-03
1648	ENDDO
1649	gamt(2) = -2.5E-03
1650	else
1651	DO k = 2, klev
1652	gamt(k) = 0.0
1653	ENDDO
1654	ENDIF
1655	cIM cf JLD/ GKtest
1656	IF ( nsrf .NE. is_oce ) THEN
1657	cIM 261103 kstable = kstable_ter
1658	kstable = ksta_ter
1659	ELSE
1660	cIM 261103 kstable = kstable_sinon
1661	kstable = ksta
1662	ENDIF
1663	cIM cf JLD/ GKtest fin
1664	c
1665	c Calculer les geopotentiels de chaque couche
1666	c
1667	DO i = 1, knon
1668	zgeop(i,1) = RD * t(i,1) / (0.5*(paprs(i,1)+pplay(i,1)))
1669	. * (paprs(i,1)-pplay(i,1))
1670	ENDDO
1671	DO k = 2, klev
1672	DO i = 1, knon
1673	zgeop(i,k) = zgeop(i,k-1)
1674	. + RD * 0.5*(t(i,k-1)+t(i,k)) / paprs(i,k)
1675	. * (pplay(i,k-1)-pplay(i,k))
1676	ENDDO
1677	ENDDO
1678	c
1679	c Calculer le frottement au sol (Cdrag)
1680	c
1681	DO i = 1, knon
1682	u1(i) = u(i,1)
1683	v1(i) = v(i,1)
1684	t1(i) = t(i,1)
1685	q1(i) = q(i,1)
1686	z1(i) = zgeop(i,1)
1687	ENDDO
1688	c
1689	CALL clcdrag(klon, knon, nsrf, zxli,
1690	$ u1, v1, t1, q1, z1,
1691	$ ts, qsurf, rugos,
1692	$ pcfm1, pcfh1)
1693	cIM $ ts, qsurf, rugos,
1694	C
1695	DO i = 1, knon
1696	pcfm(i,1)=pcfm1(i)
1697	pcfh(i,1)=pcfh1(i)
1698	ENDDO
1699	c
1700	c Calculer les coefficients turbulents dans l'atmosphere
1701	c
1702	DO i = 1, knon
1703	itop(i) = isommet
1704	ENDDO
1705
1706
1707	DO k = 2, isommet
1708	DO i = 1, knon
1709	zdu2=MAX(cepdu2,(u(i,k)-u(i,k-1))**2
1710	. +(v(i,k)-v(i,k-1))**2)
1711	zmgeom(i)=zgeop(i,k)-zgeop(i,k-1)
1712	zdphi =zmgeom(i) / 2.0
1713	zt = (t(i,k)+t(i,k-1)) * 0.5
1714	zq = (q(i,k)+q(i,k-1)) * 0.5
1715	c
1716	c calculer Qs et dQs/dT:
1717	c
1718	IF (thermcep) THEN
1719	zdelta = MAX(0.,SIGN(1.,RTT-zt))
1720	zcvm5 = R5LESRLVTT/RCPD/(1.0+RVTMP2zq)*(1.-zdelta)
1721	. + R5IESRLSTT/RCPD/(1.0+RVTMP2zq)*zdelta
1722	zqs = R2ES * FOEEW(zt,zdelta) / pplay(i,k)
1723	zqs = MIN(0.5,zqs)
1724	zcor = 1./(1.-RETV*zqs)
1725	zqs = zqs*zcor
1726	zdqs = FOEDE(zt,zdelta,zcvm5,zqs,zcor)
1727	ELSE
1728	IF (zt .LT. t_coup) THEN
1729	zqs = qsats(zt) / pplay(i,k)
1730	zdqs = dqsats(zt,zqs)
1731	ELSE
1732	zqs = qsatl(zt) / pplay(i,k)
1733	zdqs = dqsatl(zt,zqs)
1734	ENDIF
1735	ENDIF
1736	c
1737	c calculer la fraction nuageuse (processus humide):
1738	c
1739	zfr = (zq+ratqszq-zqs) / (2.0ratqs*zq)
1740	zfr = MAX(0.0,MIN(1.0,zfr))
1741	IF (.NOT.richum) zfr = 0.0
1742	c
1743	c calculer le nombre de Richardson:
1744	c
1745	IF (tvirtu) THEN
1746	ztvd =( t(i,k)
1747	. + zdphi/RCPD/(1.+RVTMP2*zq)
1748	. ( (1.-zfr) + zfr(1.+RLVTT*zqs/RD/zt)/(1.+zdqs) )
1749	. )(1.+RETVq(i,k))
1750	ztvu =( t(i,k-1)
1751	. - zdphi/RCPD/(1.+RVTMP2*zq)
1752	. ( (1.-zfr) + zfr(1.+RLVTT*zqs/RD/zt)/(1.+zdqs) )
1753	. )(1.+RETVq(i,k-1))
1754	zri(i) =zmgeom(i)(ztvd-ztvu)/(zdu20.5*(ztvd+ztvu))
1755	zri(i) = zri(i)
1756	. + zmgeom(i)zmgeom(i)/RGgamt(k)
1757	. (paprs(i,k)/101325.0)*RKAPPA
1758	. /(zdu20.5(ztvd+ztvu))
1759	c
1760	ELSE ! calcul de Ridchardson compatible LMD5
1761	c
1762	zri(i) =(RCPD*(t(i,k)-t(i,k-1))
1763	. -RD0.5(t(i,k)+t(i,k-1))/paprs(i,k)
1764	. *(pplay(i,k)-pplay(i,k-1))
1765	. )zmgeom(i)/(zdu20.5RCPD(t(i,k-1)+t(i,k)))
1766	zri(i) = zri(i) +
1767	. zmgeom(i)zmgeom(i)gamt(k)/RG
1768	cSB . /(paprs(i,k)/101325.0)**RKAPPA
1769	. (paprs(i,k)/101325.0)*RKAPPA
1770	. /(zdu20.5(t(i,k-1)+t(i,k)))
1771	ENDIF
1772	c
1773	c finalement, les coefficients d'echange sont obtenus:
1774	c
1775	zcdn=SQRT(zdu2) / zmgeom(i) * RG
1776	c
1777	IF (opt_ec) THEN
1778	z2geomf=zgeop(i,k-1)+zgeop(i,k)
1779	zalm2=(0.5ckap/RGz2geomf
1780	. /(1.+0.5ckap/rg/clamz2geomf))**2
1781	zalh2=(0.5ckap/rgz2geomf
1782	. /(1.+0.5ckap/RG/(clamSQRT(1.5cd))z2geomf))**2
1783	IF (zri(i).LT.0.0) THEN ! situation instable
1784	zscf = ((zgeop(i,k)/zgeop(i,k-1))(1./3.)-1.)3
1785	. / (zmgeom(i)/RG)**3 / (zgeop(i,k-1)/RG)
1786	zscf = SQRT(-zri(i)*zscf)
1787	zscfm = 1.0 / (1.0+3.0cbcczalm2zscf)
1788	zscfh = 1.0 / (1.0+3.0cbcczalh2zscf)
1789	pcfm(i,k)=zcdnzalm2(1.-2.0cbzri(i)*zscfm)
1790	pcfh(i,k)=zcdnzalh2(1.-3.0cbzri(i)*zscfh)
1791	ELSE ! situation stable
1792	zscf=SQRT(1.+cd*zri(i))
1793	pcfm(i,k)=zcdnzalm2/(1.+2.0cb*zri(i)/zscf)
1794	pcfh(i,k)=zcdnzalh2/(1.+3.0cbzri(i)zscf)
1795	ENDIF
1796	ELSE
1797	zl2(i)=(mixlen*MAX(0.0,(paprs(i,k)-paprs(i,itop(i)+1))
1798	. /(paprs(i,2)-paprs(i,itop(i)+1)) ))**2
1799	pcfm(i,k)=sqrt(max(zcdnzcdn(ric-zri(i))/ric, kstable))
1800	pcfm(i,k)= zl2(i)* pcfm(i,k)
1801	pcfh(i,k) = pcfm(i,k) /prandtl ! h et m different
1802	ENDIF
1803	ENDDO
1804	ENDDO
1805	c
1806	c Au-dela du sommet, pas de diffusion turbulente:
1807	c
1808	DO i = 1, knon
1809	IF (itop(i)+1 .LE. klev) THEN
1810	DO k = itop(i)+1, klev
1811	pcfh(i,k) = 0.0
1812	pcfm(i,k) = 0.0
1813	ENDDO
1814	ENDIF
1815	ENDDO
1816	c
1817	RETURN
1818	END
1819
1820	SUBROUTINE coefkz2(nsrf, knon, paprs, pplay,t,
1821	. pcfm, pcfh)
1822	IMPLICIT none
1823	c======================================================================
1824	c J'introduit un peu de diffusion sauf dans les endroits
1825	c ou une forte inversion est presente
1826	c On peut dire qu'il represente la convection peu profonde
1827	c
1828	c Arguments:
1829	c nsrf-----input-I- indicateur de la nature du sol
1830	c knon-----input-I- nombre de points a traiter
1831	c paprs----input-R- pression a chaque intercouche (en Pa)
1832	c pplay----input-R- pression au milieu de chaque couche (en Pa)
1833	c t--------input-R- temperature (K)
1834	c
1835	c pcfm-----output-R- coefficients a calculer (vitesse)
1836	c pcfh-----output-R- coefficients a calculer (chaleur et humidite)
1837	c======================================================================
1838	#include "dimensions.h"
1839	#include "dimphy.h"
1840	#include "YOMCST.h"
1841	#include "indicesol.h"
1842	#include "iniprint.h"
1843	c
1844	c Arguments:
1845	c
1846	INTEGER knon, nsrf
1847	REAL paprs(klon,klev+1), pplay(klon,klev)
1848	REAL t(klon,klev)
1849	c
1850	REAL pcfm(klon,klev), pcfh(klon,klev)
1851	c
1852	c Quelques constantes et options:
1853	c
1854	REAL prandtl
1855	PARAMETER (prandtl=0.4)
1856	REAL kstable
1857	PARAMETER (kstable=0.002)
1858	ccc PARAMETER (kstable=0.001)
1859	REAL mixlen ! constante controlant longueur de melange
1860	PARAMETER (mixlen=35.0)
1861	REAL seuil ! au-dela l'inversion est consideree trop faible
1862	PARAMETER (seuil=-0.02)
1863	ccc PARAMETER (seuil=-0.04)
1864	ccc PARAMETER (seuil=-0.06)
1865	ccc PARAMETER (seuil=-0.09)
1866	c
1867	c Variables locales:
1868	c
1869	INTEGER i, k, invb(knon)
1870	REAL zl2(knon)
1871	REAL zdthmin(knon), zdthdp
1872	c
1873	c Initialiser les sorties
1874	c
1875	DO k = 1, klev
1876	DO i = 1, knon
1877	pcfm(i,k) = 0.0
1878	pcfh(i,k) = 0.0
1879	ENDDO
1880	ENDDO
1881	c
1882	c Chercher la zone d'inversion forte
1883	c
1884	DO i = 1, knon
1885	invb(i) = klev
1886	zdthmin(i)=0.0
1887	ENDDO
1888	DO k = 2, klev/2-1
1889	DO i = 1, knon
1890	zdthdp = (t(i,k)-t(i,k+1))/(pplay(i,k)-pplay(i,k+1))
1891	. - RD * 0.5*(t(i,k)+t(i,k+1))/RCPD/paprs(i,k+1)
1892	zdthdp = zdthdp * 100.0
1893	IF (pplay(i,k).GT.0.8*paprs(i,1) .AND.
1894	. zdthdp.LT.zdthmin(i) ) THEN
1895	zdthmin(i) = zdthdp
1896	invb(i) = k
1897	ENDIF
1898	ENDDO
1899	ENDDO
1900	c
1901	c Introduire une diffusion:
1902	c
1903	DO k = 2, klev
1904	DO i = 1, knon
1905	cIM cf FH/GK IF ( (nsrf.NE.is_oce) .OR. ! si ce n'est pas sur l'ocean
1906	cIM cf FH/GK . (invb(i).EQ.klev) .OR. ! s'il n'y a pas d'inversion
1907	!IM cf JLD/ GKtest TERkz2
1908	! IF ( (nsrf.EQ.is_ter) .OR. ! si on est sur la terre
1909	! fin GKtest
1910	IF ( (nsrf.EQ.is_oce) .AND. ! si on est sur ocean et si
1911	. ( (invb(i).EQ.klev) .OR. ! s'il n'y a pas d'inversion
1912	. (zdthmin(i).GT.seuil) ) )THEN ! si l'inversion est trop faible
1913	zl2(i)=(mixlen*MAX(0.0,(paprs(i,k)-paprs(i,klev+1))
1914	. /(paprs(i,2)-paprs(i,klev+1)) ))**2
1915	pcfm(i,k)= zl2(i)* kstable
1916	pcfh(i,k) = pcfm(i,k) /prandtl ! h et m different
1917	ENDIF
1918	ENDDO
1919	ENDDO
1920	c
1921	RETURN
1922	END
1923	SUBROUTINE calbeta(dtime,indice,knon,snow,qsol,
1924	. vbeta,vcal,vdif)
1925	IMPLICIT none
1926	c======================================================================
1927	c Auteur(s): Z.X. Li (LMD/CNRS) (adaptation du GCM du LMD)
1928	c date: 19940414
1929	c======================================================================
1930	c
1931	c Calculer quelques parametres pour appliquer la couche limite
1932	c ------------------------------------------------------------
1933	#include "dimensions.h"
1934	#include "dimphy.h"
1935	#include "YOMCST.h"
1936	#include "indicesol.h"
1937	#include "iniprint.h"
1938	REAL tau_gl ! temps de relaxation pour la glace de mer
1939	ccc PARAMETER (tau_gl=86400.0*30.0)
1940	PARAMETER (tau_gl=86400.0*5.0)
1941	REAL mx_eau_sol
1942	PARAMETER (mx_eau_sol=150.0)
1943	c
1944	REAL calsol, calsno, calice ! epaisseur du sol: 0.15 m
1945	PARAMETER (calsol=1.0/(2.5578E+06*0.15))
1946	PARAMETER (calsno=1.0/(2.3867E+06*0.15))
1947	PARAMETER (calice=1.0/(5.1444E+06*0.15))
1948	C
1949	INTEGER i
1950	c
1951	REAL dtime
1952	REAL snow(klon), qsol(klon)
1953	INTEGER indice, knon
1954	C
1955	REAL vbeta(klon)
1956	REAL vcal(klon)
1957	REAL vdif(klon)
1958	C
1959
1960	IF (indice.EQ.is_oce) THEN
1961	DO i = 1, knon
1962	vcal(i) = 0.0
1963	vbeta(i) = 1.0
1964	vdif(i) = 0.0
1965	ENDDO
1966	ENDIF
1967	c
1968	IF (indice.EQ.is_sic) THEN
1969	DO i = 1, knon
1970	vcal(i) = calice
1971	IF (snow(i) .GT. 0.0) vcal(i) = calsno
1972	vbeta(i) = 1.0
1973	vdif(i) = 1.0/tau_gl
1974	ccc vdif(i) = calice/tau_gl ! c'etait une erreur
1975	ENDDO
1976	ENDIF
1977	c
1978	IF (indice.EQ.is_ter) THEN
1979	DO i = 1, knon
1980	vcal(i) = calsol
1981	IF (snow(i) .GT. 0.0) vcal(i) = calsno
1982	vbeta(i) = MIN(2.0*qsol(i)/mx_eau_sol, 1.0)
1983	vdif(i) = 0.0
1984	ENDDO
1985	ENDIF
1986	c
1987	IF (indice.EQ.is_lic) THEN
1988	DO i = 1, knon
1989	vcal(i) = calice
1990	IF (snow(i) .GT. 0.0) vcal(i) = calsno
1991	vbeta(i) = 1.0
1992	vdif(i) = 0.0
1993	ENDDO
1994	ENDIF
1995	c
1996	RETURN
1997	END
1998	C======================================================================
1999	SUBROUTINE nonlocal(knon, paprs, pplay,
2000	. tsol,beta,u,v,t,q,
2001	. cd_h, cd_m, pcfh, pcfm, cgh, cgq)
2002	IMPLICIT none
2003	c======================================================================
2004	c Laurent Li (LMD/CNRS), le 30 septembre 1998
2005	c Couche limite non-locale. Adaptation du code du CCM3.
2006	c Code non teste, donc a ne pas utiliser.
2007	c======================================================================
2008	c Nonlocal scheme that determines eddy diffusivities based on a
2009	c diagnosed boundary layer height and a turbulent velocity scale.
2010	c Also countergradient effects for heat and moisture are included.
2011	c
2012	c For more information, see Holtslag, A.A.M., and B.A. Boville, 1993:
2013	c Local versus nonlocal boundary-layer diffusion in a global climate
2014	c model. J. of Climate, vol. 6, 1825-1842.
2015	c======================================================================
2016	#include "dimensions.h"
2017	#include "dimphy.h"
2018	#include "YOMCST.h"
2019	#include "iniprint.h"
2020	c
2021	c Arguments:
2022	c
2023	INTEGER knon ! nombre de points a calculer
2024	REAL tsol(klon) ! temperature du sol (K)
2025	REAL beta(klon) ! efficacite d'evaporation (entre 0 et 1)
2026	REAL paprs(klon,klev+1) ! pression a inter-couche (Pa)
2027	REAL pplay(klon,klev) ! pression au milieu de couche (Pa)
2028	REAL u(klon,klev) ! vitesse U (m/s)
2029	REAL v(klon,klev) ! vitesse V (m/s)
2030	REAL t(klon,klev) ! temperature (K)
2031	REAL q(klon,klev) ! vapeur d'eau (kg/kg)
2032	REAL cd_h(klon) ! coefficient de friction au sol pour chaleur
2033	REAL cd_m(klon) ! coefficient de friction au sol pour vitesse
2034	c
2035	INTEGER isommet
2036	PARAMETER (isommet=klev)
2037	REAL vk
2038	PARAMETER (vk=0.40)
2039	REAL ricr
2040	PARAMETER (ricr=0.4)
2041	REAL fak
2042	PARAMETER (fak=8.5)
2043	REAL fakn
2044	PARAMETER (fakn=7.2)
2045	REAL onet
2046	PARAMETER (onet=1.0/3.0)
2047	REAL t_coup
2048	PARAMETER(t_coup=273.15)
2049	REAL zkmin
2050	PARAMETER (zkmin=0.01)
2051	REAL betam
2052	PARAMETER (betam=15.0)
2053	REAL betah
2054	PARAMETER (betah=15.0)
2055	REAL betas
2056	PARAMETER (betas=5.0)
2057	REAL sffrac
2058	PARAMETER (sffrac=0.1)
2059	REAL binm
2060	PARAMETER (binm=betam*sffrac)
2061	REAL binh
2062	PARAMETER (binh=betah*sffrac)
2063	REAL ccon
2064	PARAMETER (ccon=faksffracvk)
2065	c
2066	REAL z(klon,klev)
2067	REAL pcfm(klon,klev), pcfh(klon,klev)
2068	c
2069	INTEGER i, k
2070	REAL zxt, zxq, zxu, zxv, zxmod, taux, tauy
2071	REAL zx_alf1, zx_alf2 ! parametres pour extrapolation
2072	REAL khfs(klon) ! surface kinematic heat flux [mK/s]
2073	REAL kqfs(klon) ! sfc kinematic constituent flux [m/s]
2074	REAL heatv(klon) ! surface virtual heat flux
2075	REAL ustar(klon)
2076	REAL rino(klon,klev) ! bulk Richardon no. from level to ref lev
2077	LOGICAL unstbl(klon) ! pts w/unstbl pbl (positive virtual ht flx)
2078	LOGICAL stblev(klon) ! stable pbl with levels within pbl
2079	LOGICAL unslev(klon) ! unstbl pbl with levels within pbl
2080	LOGICAL unssrf(klon) ! unstb pbl w/lvls within srf pbl lyr
2081	LOGICAL unsout(klon) ! unstb pbl w/lvls in outer pbl lyr
2082	LOGICAL check(klon) ! True=>chk if Richardson no.>critcal
2083	REAL pblh(klon)
2084	REAL cgh(klon,2:klev) ! counter-gradient term for heat [K/m]
2085	REAL cgq(klon,2:klev) ! counter-gradient term for constituents
2086	REAL cgs(klon,2:klev) ! counter-gradient star (cg/flux)
2087	REAL obklen(klon)
2088	REAL ztvd, ztvu, zdu2
2089	REAL therm(klon) ! thermal virtual temperature excess
2090	REAL phiminv(klon) ! inverse phi function for momentum
2091	REAL phihinv(klon) ! inverse phi function for heat
2092	REAL wm(klon) ! turbulent velocity scale for momentum
2093	REAL fak1(klon) ! kustarpblh
2094	REAL fak2(klon) ! kwmpblh
2095	REAL fak3(klon) ! fakn*wstr/wm
2096	REAL pblk(klon) ! level eddy diffusivity for momentum
2097	REAL pr(klon) ! Prandtl number for eddy diffusivities
2098	REAL zl(klon) ! zmzp / Obukhov length
2099	REAL zh(klon) ! zmzp / pblh
2100	REAL zzh(klon) ! (1-(zmzp/pblh))**2
2101	REAL wstr(klon) ! w*, convective velocity scale
2102	REAL zm(klon) ! current level height
2103	REAL zp(klon) ! current level height + one level up
2104	REAL zcor, zdelta, zcvm5, zxqs
2105	REAL fac, pblmin, zmzp, term
2106	c
2107	#include "YOETHF.h"
2108	#include "FCTTRE.h"
2109	c
2110	c Initialisation
2111	c
2112	DO i = 1, klon
2113	pcfh(i,1) = cd_h(i)
2114	pcfm(i,1) = cd_m(i)
2115	ENDDO
2116	DO k = 2, klev
2117	DO i = 1, klon
2118	pcfh(i,k) = zkmin
2119	pcfm(i,k) = zkmin
2120	cgs(i,k) = 0.0
2121	cgh(i,k) = 0.0
2122	cgq(i,k) = 0.0
2123	ENDDO
2124	ENDDO
2125	c
2126	c Calculer les hauteurs de chaque couche
2127	c
2128	DO i = 1, knon
2129	z(i,1) = RD * t(i,1) / (0.5*(paprs(i,1)+pplay(i,1)))
2130	. * (paprs(i,1)-pplay(i,1)) / RG
2131	ENDDO
2132	DO k = 2, klev
2133	DO i = 1, knon
2134	z(i,k) = z(i,k-1)
2135	. + RD * 0.5*(t(i,k-1)+t(i,k)) / paprs(i,k)
2136	. * (pplay(i,k-1)-pplay(i,k)) / RG
2137	ENDDO
2138	ENDDO
2139	c
2140	DO i = 1, knon
2141	IF (thermcep) THEN
2142	zdelta=MAX(0.,SIGN(1.,RTT-tsol(i)))
2143	zcvm5 = R5LESRLVTT(1.-zdelta) + R5IESRLSTTzdelta
2144	zcvm5 = zcvm5 / RCPD / (1.0+RVTMP2*q(i,1))
2145	zxqs= r2es * FOEEW(tsol(i),zdelta)/paprs(i,1)
2146	zxqs=MIN(0.5,zxqs)
2147	zcor=1./(1.-retv*zxqs)
2148	zxqs=zxqs*zcor
2149	ELSE
2150	IF (tsol(i).LT.t_coup) THEN
2151	zxqs = qsats(tsol(i)) / paprs(i,1)
2152	ELSE
2153	zxqs = qsatl(tsol(i)) / paprs(i,1)
2154	ENDIF
2155	ENDIF
2156	zx_alf1 = 1.0
2157	zx_alf2 = 1.0 - zx_alf1
2158	zxt = (t(i,1)+z(i,1)RG/RCPD/(1.+RVTMP2q(i,1)))
2159	. (1.+RETVq(i,1))*zx_alf1
2160	. + (t(i,2)+z(i,2)RG/RCPD/(1.+RVTMP2q(i,2)))
2161	. (1.+RETVq(i,2))*zx_alf2
2162	zxu = u(i,1)zx_alf1+u(i,2)zx_alf2
2163	zxv = v(i,1)zx_alf1+v(i,2)zx_alf2
2164	zxq = q(i,1)zx_alf1+q(i,2)zx_alf2
2165	zxmod = 1.0+SQRT(zxu2+zxv2)
2166	khfs(i) = (tsol(i)(1.+RETVq(i,1))-zxt) zxmodcd_h(i)
2167	kqfs(i) = (zxqs-zxq) zxmodcd_h(i) * beta(i)
2168	heatv(i) = khfs(i) + 0.61zxtkqfs(i)
2169	taux = zxu zxmodcd_m(i)
2170	tauy = zxv zxmodcd_m(i)
2171	ustar(i) = SQRT(taux2+tauy2)
2172	ustar(i) = MAX(SQRT(ustar(i)),0.01)
2173	ENDDO
2174	c
2175	DO i = 1, knon
2176	rino(i,1) = 0.0
2177	check(i) = .TRUE.
2178	pblh(i) = z(i,1)
2179	obklen(i) = -t(i,1)ustar(i)3/(RGvk*heatv(i))
2180	ENDDO
2181
2182	C
2183	C PBL height calculation:
2184	C Search for level of pbl. Scan upward until the Richardson number between
2185	C the first level and the current level exceeds the "critical" value.
2186	C
2187	fac = 100.0
2188	DO k = 1, isommet
2189	DO i = 1, knon
2190	IF (check(i)) THEN
2191	zdu2 = (u(i,k)-u(i,1))2+(v(i,k)-v(i,1))2+facustar(i)*2
2192	zdu2 = max(zdu2,1.0e-20)
2193	ztvd =(t(i,k)+z(i,k)0.5RG/RCPD/(1.+RVTMP2*q(i,k)))
2194	. (1.+RETVq(i,k))
2195	ztvu =(t(i,1)-z(i,k)0.5RG/RCPD/(1.+RVTMP2*q(i,1)))
2196	. (1.+RETVq(i,1))
2197	rino(i,k) = (z(i,k)-z(i,1))RG(ztvd-ztvu)
2198	. /(zdu20.5(ztvd+ztvu))
2199	IF (rino(i,k).GE.ricr) THEN
2200	pblh(i) = z(i,k-1) + (z(i,k-1)-z(i,k)) *
2201	. (ricr-rino(i,k-1))/(rino(i,k-1)-rino(i,k))
2202	check(i) = .FALSE.
2203	ENDIF
2204	ENDIF
2205	ENDDO
2206	ENDDO
2207
2208	C
2209	C Set pbl height to maximum value where computation exceeds number of
2210	C layers allowed
2211	C
2212	DO i = 1, knon
2213	if (check(i)) pblh(i) = z(i,isommet)
2214	ENDDO
2215	C
2216	C Improve estimate of pbl height for the unstable points.
2217	C Find unstable points (sensible heat flux is upward):
2218	C
2219	DO i = 1, knon
2220	IF (heatv(i) .GT. 0.) THEN
2221	unstbl(i) = .TRUE.
2222	check(i) = .TRUE.
2223	ELSE
2224	unstbl(i) = .FALSE.
2225	check(i) = .FALSE.
2226	ENDIF
2227	ENDDO
2228	C
2229	C For the unstable case, compute velocity scale and the
2230	C convective temperature excess:
2231	C
2232	DO i = 1, knon
2233	IF (check(i)) THEN
2234	phiminv(i) = (1.-binmpblh(i)/obklen(i))*onet
2235	wm(i)= ustar(i)*phiminv(i)
2236	therm(i) = heatv(i)*fak/wm(i)
2237	rino(i,1) = 0.0
2238	ENDIF
2239	ENDDO
2240	C
2241	C Improve pblh estimate for unstable conditions using the
2242	C convective temperature excess:
2243	C
2244	DO k = 1, isommet
2245	DO i = 1, knon
2246	IF (check(i)) THEN
2247	zdu2 = (u(i,k)-u(i,1))2+(v(i,k)-v(i,1))2+facustar(i)*2
2248	zdu2 = max(zdu2,1.0e-20)
2249	ztvd =(t(i,k)+z(i,k)0.5RG/RCPD/(1.+RVTMP2*q(i,k)))
2250	. (1.+RETVq(i,k))
2251	ztvu =(t(i,1)+therm(i)-z(i,k)0.5RG/RCPD/(1.+RVTMP2*q(i,1)))
2252	. (1.+RETVq(i,1))
2253	rino(i,k) = (z(i,k)-z(i,1))RG(ztvd-ztvu)
2254	. /(zdu20.5(ztvd+ztvu))
2255	IF (rino(i,k).GE.ricr) THEN
2256	pblh(i) = z(i,k-1) + (z(i,k-1)-z(i,k)) *
2257	. (ricr-rino(i,k-1))/(rino(i,k-1)-rino(i,k))
2258	check(i) = .FALSE.
2259	ENDIF
2260	ENDIF
2261	ENDDO
2262	ENDDO
2263	C
2264	C Set pbl height to maximum value where computation exceeds number of
2265	C layers allowed
2266	C
2267	DO i = 1, knon
2268	if (check(i)) pblh(i) = z(i,isommet)
2269	ENDDO
2270	C
2271	C Points for which pblh exceeds number of pbl layers allowed;
2272	C set to maximum
2273	C
2274	DO i = 1, knon
2275	IF (check(i)) pblh(i) = z(i,isommet)
2276	ENDDO
2277	C
2278	C PBL height must be greater than some minimum mechanical mixing depth
2279	C Several investigators have proposed minimum mechanical mixing depth
2280	C relationships as a function of the local friction velocity, u*. We
2281	C make use of a linear relationship of the form h = c u* where c=700.
2282	C The scaling arguments that give rise to this relationship most often
2283	C represent the coefficient c as some constant over the local coriolis
2284	C parameter. Here we make use of the experimental results of Koracin
2285	C and Berkowicz (1988) [BLM, Vol 43] for wich they recommend 0.07/f
2286	C where f was evaluated at 39.5 N and 52 N. Thus we use a typical mid
2287	C latitude value for f so that c = 0.07/f = 700.
2288	C
2289	DO i = 1, knon
2290	pblmin = 700.0*ustar(i)
2291	pblh(i) = MAX(pblh(i),pblmin)
2292	ENDDO
2293	C
2294	C pblh is now available; do preparation for diffusivity calculation:
2295	C
2296	DO i = 1, knon
2297	pblk(i) = 0.0
2298	fak1(i) = ustar(i)pblh(i)vk
2299	C
2300	C Do additional preparation for unstable cases only, set temperature
2301	C and moisture perturbations depending on stability.
2302	C
2303	IF (unstbl(i)) THEN
2304	zxt=(t(i,1)-z(i,1)0.5RG/RCPD/(1.+RVTMP2*q(i,1)))
2305	. (1.+RETVq(i,1))
2306	phiminv(i) = (1. - binmpblh(i)/obklen(i))*onet
2307	phihinv(i) = sqrt(1. - binh*pblh(i)/obklen(i))
2308	wm(i) = ustar(i)*phiminv(i)
2309	fak2(i) = wm(i)pblh(i)vk
2310	wstr(i) = (heatv(i)RGpblh(i)/zxt)**onet
2311	fak3(i) = fakn*wstr(i)/wm(i)
2312	ENDIF
2313	ENDDO
2314
2315	C Main level loop to compute the diffusivities and
2316	C counter-gradient terms:
2317	C
2318	DO 1000 k = 2, isommet
2319	C
2320	C Find levels within boundary layer:
2321	C
2322	DO i = 1, knon
2323	unslev(i) = .FALSE.
2324	stblev(i) = .FALSE.
2325	zm(i) = z(i,k-1)
2326	zp(i) = z(i,k)
2327	IF (zkmin.EQ.0.0 .AND. zp(i).GT.pblh(i)) zp(i) = pblh(i)
2328	IF (zm(i) .LT. pblh(i)) THEN
2329	zmzp = 0.5*(zm(i) + zp(i))
2330	zh(i) = zmzp/pblh(i)
2331	zl(i) = zmzp/obklen(i)
2332	zzh(i) = 0.
2333	IF (zh(i).LE.1.0) zzh(i) = (1. - zh(i))**2
2334	C
2335	C stblev for points zm < plbh and stable and neutral
2336	C unslev for points zm < plbh and unstable
2337	C
2338	IF (unstbl(i)) THEN
2339	unslev(i) = .TRUE.
2340	ELSE
2341	stblev(i) = .TRUE.
2342	ENDIF
2343	ENDIF
2344	ENDDO
2345	C
2346	C Stable and neutral points; set diffusivities; counter-gradient
2347	C terms zero for stable case:
2348	C
2349	DO i = 1, knon
2350	IF (stblev(i)) THEN
2351	IF (zl(i).LE.1.) THEN
2352	pblk(i) = fak1(i)zh(i)zzh(i)/(1. + betas*zl(i))
2353	ELSE
2354	pblk(i) = fak1(i)zh(i)zzh(i)/(betas + zl(i))
2355	ENDIF
2356	pcfm(i,k) = pblk(i)
2357	pcfh(i,k) = pcfm(i,k)
2358	ENDIF
2359	ENDDO
2360	C
2361	C unssrf, unstable within surface layer of pbl
2362	C unsout, unstable within outer layer of pbl
2363	C
2364	DO i = 1, knon
2365	unssrf(i) = .FALSE.
2366	unsout(i) = .FALSE.
2367	IF (unslev(i)) THEN
2368	IF (zh(i).lt.sffrac) THEN
2369	unssrf(i) = .TRUE.
2370	ELSE
2371	unsout(i) = .TRUE.
2372	ENDIF
2373	ENDIF
2374	ENDDO
2375	C
2376	C Unstable for surface layer; counter-gradient terms zero
2377	C
2378	DO i = 1, knon
2379	IF (unssrf(i)) THEN
2380	term = (1. - betamzl(i))*onet
2381	pblk(i) = fak1(i)zh(i)zzh(i)*term
2382	pr(i) = term/sqrt(1. - betah*zl(i))
2383	ENDIF
2384	ENDDO
2385	C
2386	C Unstable for outer layer; counter-gradient terms non-zero:
2387	C
2388	DO i = 1, knon
2389	IF (unsout(i)) THEN
2390	pblk(i) = fak2(i)zh(i)zzh(i)
2391	cgs(i,k) = fak3(i)/(pblh(i)*wm(i))
2392	cgh(i,k) = khfs(i)*cgs(i,k)
2393	pr(i) = phiminv(i)/phihinv(i) + ccon*fak3(i)/fak
2394	cgq(i,k) = kqfs(i)*cgs(i,k)
2395	ENDIF
2396	ENDDO
2397	C
2398	C For all unstable layers, set diffusivities
2399	C
2400	DO i = 1, knon
2401	IF (unslev(i)) THEN
2402	pcfm(i,k) = pblk(i)
2403	pcfh(i,k) = pblk(i)/pr(i)
2404	ENDIF
2405	ENDDO
2406	1000 continue ! end of level loop
2407
2408	RETURN
2409	END

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

Download in other formats:

Original Format