Context Navigation

surface_mod.F @ 149

Last change on this file since 149 was 149, checked in by sdubey, 11 years ago

Added few new routines to read NC files and compute diagnostics to r145.

Few routines of dry physics including radiation module, surface process and convective adjustment in new routine phyparam.f90. dynetat to read start files for dynamics.

check_conserve routine to compute conservation of quatities like mass, energy etc.etat0_heldsz.f90 for held-suarez test case initial conditions.

new Key time_style=lmd or dcmip to use day_step, ndays like in LMDZ

File size: 27.0 KB

Line
1	MODULE SURFACE_PROCESS
2
3	USE ICOSA
4	USE dimphys_mod
5	USE RADIATION
6	DATA lmixmin,emin_turb,karman/100.,1.e-8,.4/
7
8	contains
9
10	!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
11
12	subroutiNE vdif(ngrid,nlay,ptime,
13	$ ptimestep,pcapcal,pz0,
14	$ pplay,pplev,pzlay,pzlev,
15	$ pu,pv,ph,ptsrf,pemis,
16	$ pdufi,pdvfi,pdhfi,pfluxsrf,
17	$ pdudif,pdvdif,pdhdif,pdtsrf,pq2,pq2l,
18	$ lwrite)
19	IMPLICIT NONE
20
21	c=======================================================================
22	c
23	c Diffusion verticale
24	c Shema implicite
25	c On commence par rajouter au variables x la tendance physique
26	c et on resoult en fait:
27	c x(t+1) = x(t) + dt * (dx/dt)phys(t) + dt * (dx/dt)difv(t+1)
28	c
29	c !!! attention :
30	c pour utilisation sur une machine sans allocation dynamique de
31	c memoires (sur SUN par exemple) il faut que ngrid soit egal
32	c a ngrid.
33	c
34	c arguments:
35	c ----------
36	c
37	c entree:
38	c -------
39	c
40	c
41	c=======================================================================
42
43	c-----------------------------------------------------------------------
44	c declarations:
45	c -------------
46	c
47	c arguments:
48	c ----------
49
50	INTEGER ngrid,nlay
51	REAL ptime,ptimestep
52	REAL pplay(ngrid,nlay),pplev(ngrid,nlay+1)
53	REAL pzlay(ngrid,nlay),pzlev(ngrid,nlay+1)
54	REAL pu(ngrid,nlay),pv(ngrid,nlay),ph(ngrid,nlay)
55	REAL ptsrf(ngrid),pemis(ngrid)
56	REAL pdufi(ngrid,nlay),pdvfi(ngrid,nlay),pdhfi(ngrid,nlay)
57	REAL pfluxsrf(ngrid)
58	REAL pdudif(ngrid,nlay),pdvdif(ngrid,nlay),pdhdif(ngrid,nlay)
59	REAL pdtsrf(ngrid),pcapcal(ngrid),pz0(ngrid)
60	REAL pq2(ngrid,nlay+1),pq2l(ngrid,nlay+1)
61	LOGICAL lwrite
62	c
63	c local:
64	c ------
65
66	INTEGER ilev,ig,ilay,nlev
67	INTEGER unit,ierr,it1,it2,icount
68	SAVE icount
69	INTEGER cluvdb,putdat,putvdim,setname,setvdim
70	REAL z4st,zdplanck(ngrid),zu2
71	REAL zkv(ngrid,nlayermx+1),zkh(ngrid,nlayermx+1)
72	REAL zcdv(ngrid),zcdh(ngrid)
73	REAL zu(ngrid,nlayermx),zv(ngrid,nlayermx)
74	REAL zh(ngrid,nlayermx)
75	REAL ztsrf2(ngrid)
76	REAL z1(ngrid),z2(ngrid)
77	REAL za(ngrid,nlayermx),zb(ngrid,nlayermx)
78	REAL zb0(ngrid,nlayermx)
79	REAL zc(ngrid,nlayermx),zd(ngrid,nlayermx)
80	REAL zout_dyn(iim+1,jjm+1,nlayermx+1),zout_fi(ngrid,nlayermx+1)
81	REAL zcst1
82	REAL karman
83
84	EXTERNAL coefdifv
85	EXTERNAL SSUM
86	REAL SSUM
87	SAVE karman
88
89	DATA karman/0.4/
90	DATA icount/0/
91	c
92	c-----------------------------------------------------------------------
93	c initialisations:
94	c ----------------
95
96	nlev=nlay+1
97
98	IF(ngrid.NE.ngrid) THEN
99	PRINT*,'STOP dans coefdifv'
100	PRINT*,'probleme de dimensions :'
101	PRINT*,'ngrid =',ngrid
102	PRINT*,'ngrid =',ngrid
103	STOP
104	ENDIF
105
106	c computation of rhodz and dtrho/dz=dtrho*2 g/dp:
107	c with rho=p/RT=p/ (R Theta) (p/ps)**kappa
108	c ---------------------------------
109
110	DO ilay=1,nlay
111	DO ig=1,ngrid
112	za(ig,ilay)=
113	s (pplev(ig,ilay)-pplev(ig,ilay+1))/g
114	ENDDO
115	ENDDO
116
117	zcst1=4.gptimestep/(kappacpp)*2
118	DO ilev=2,nlev-1
119	DO ig=1,ngrid
120	zb0(ig,ilev)=pplev(ig,ilev)*
121	s (pplev(ig,1)/pplev(ig,ilev))**kappa /
122	s (ph(ig,ilev-1)+ph(ig,ilev))
123	zb0(ig,ilev)=zcst1zb0(ig,ilev)zb0(ig,ilev)/
124	s (pplay(ig,ilev-1)-pplay(ig,ilev))
125	ENDDO
126	ENDDO
127	DO ig=1,ngrid
128	zb0(ig,1)=ptimesteppplev(ig,1)/(kappacpp*ptsrf(ig))
129	ENDDO
130	IF(lwrite) THEN
131	ig=ngrid/2+1
132	PRINT*,'Pression (mbar) ,altitude (km),u,v,theta, rho dz'
133	DO ilay=1,nlay
134	WRITE(,) .01pplay(ig,ilay),.001pzlay(ig,ilay),
135	s pu(ig,ilay),pv(ig,ilay),ph(ig,ilay),za(ig,ilay)
136	ENDDO
137	PRINT*,'Pression (mbar) ,altitude (km),zb'
138	DO ilev=1,nlay
139	WRITE(,) .01pplev(ig,ilev),.001pzlev(ig,ilev),
140	s zb0(ig,ilev)
141	ENDDO
142	ENDIF
143
144	c-----------------------------------------------------------------------
145	c 2. ajout des tendances physiques:
146	c ------------------------------
147
148	DO ilev=1,nlay
149	DO ig=1,ngrid
150	zu(ig,ilev)=pu(ig,ilev)+pdufi(ig,ilev)*ptimestep
151	zv(ig,ilev)=pv(ig,ilev)+pdvfi(ig,ilev)*ptimestep
152	zh(ig,ilev)=ph(ig,ilev)+pdhfi(ig,ilev)*ptimestep
153	ENDDO
154	ENDDO
155
156	c-----------------------------------------------------------------------
157	c 3. calcul de cd :
158	c ----------------
159	c
160	CALL vdif_cd( ngrid,pz0,g,pzlay,pu,pv,ptsrf,ph,zcdv,zcdh)
161
162	c CALL my_25(ptimestep,g,pzlev,pzlay,pu,pv,ph,zcdv,
163	c a pq2,pq2l,zkv,zkh)
164
165	CALL vdif_k(ngrid,nlay,
166	s ptimestep,g,pzlev,pzlay,pz0,pu,pv,ph,zcdv,zkv,zkh,pq2,pq2l)
167
168	DO ig=1,ngrid
169	zu2=pu(ig,1)pu(ig,1)+pv(ig,1)pv(ig,1)
170	zcdv(ig)=zcdv(ig)*sqrt(zu2)
171	zcdh(ig)=zcdh(ig)*sqrt(zu2)
172	ENDDO
173
174	IF(lwrite) THEN
175	PRINT*
176	PRINT*,'Diagnostique diffusion verticale'
177	PRINT*,'coefficients Cd pour v et h'
178	PRINT*,zcdv(ngrid/2+1),zcdh(ngrid/2+1)
179	PRINT*,'coefficients K pour v et h'
180	DO ilev=1,nlay
181	PRINT*,zkv(ngrid/2+1,ilev),zkh(ngrid/2+1,ilev)
182	ENDDO
183	ENDIF
184
185	c-----------------------------------------------------------------------
186	c integration verticale pour u:
187	c -----------------------------
188	c
189	CALL multipl((nlay-1)*ngrid,zkv(1,2),zb0(1,2),zb(1,2))
190	CALL multipl(ngrid,zcdv,zb0,zb)
191	DO ig=1,ngrid
192	z1(ig)=1./(za(ig,nlay)+zb(ig,nlay))
193	zc(ig,nlay)=za(ig,nlay)zu(ig,nlay)z1(ig)
194	zd(ig,nlay)=zb(ig,nlay)*z1(ig)
195	ENDDO
196
197	DO ilay=nlay-1,1,-1
198	DO ig=1,ngrid
199	z1(ig)=1./(za(ig,ilay)+zb(ig,ilay)+
200	$ zb(ig,ilay+1)*(1.-zd(ig,ilay+1)))
201	zc(ig,ilay)=(za(ig,ilay)*zu(ig,ilay)+
202	$ zb(ig,ilay+1)zc(ig,ilay+1))z1(ig)
203	zd(ig,ilay)=zb(ig,ilay)*z1(ig)
204	ENDDO
205	ENDDO
206
207	DO ig=1,ngrid
208	zu(ig,1)=zc(ig,1)
209	ENDDO
210	DO ilay=2,nlay
211	DO ig=1,ngrid
212	zu(ig,ilay)=zc(ig,ilay)+zd(ig,ilay)*zu(ig,ilay-1)
213	ENDDO
214	ENDDO
215
216	c-----------------------------------------------------------------------
217	c integration verticale pour v:
218	c -----------------------------
219	c
220	DO ig=1,ngrid
221	z1(ig)=1./(za(ig,nlay)+zb(ig,nlay))
222	zc(ig,nlay)=za(ig,nlay)zv(ig,nlay)z1(ig)
223	zd(ig,nlay)=zb(ig,nlay)*z1(ig)
224	ENDDO
225
226	DO ilay=nlay-1,1,-1
227	DO ig=1,ngrid
228	z1(ig)=1./(za(ig,ilay)+zb(ig,ilay)+
229	$ zb(ig,ilay+1)*(1.-zd(ig,ilay+1)))
230	zc(ig,ilay)=(za(ig,ilay)*zv(ig,ilay)+
231	$ zb(ig,ilay+1)zc(ig,ilay+1))z1(ig)
232	zd(ig,ilay)=zb(ig,ilay)*z1(ig)
233	ENDDO
234	ENDDO
235
236	DO ig=1,ngrid
237	zv(ig,1)=zc(ig,1)
238	ENDDO
239	DO ilay=2,nlay
240	DO ig=1,ngrid
241	zv(ig,ilay)=zc(ig,ilay)+zd(ig,ilay)*zv(ig,ilay-1)
242	ENDDO
243	ENDDO
244
245	c-----------------------------------------------------------------------
246	c integration verticale pour h:
247	c -----------------------------
248	c
249	CALL multipl((nlay-1)*ngrid,zkh(1,2),zb0(1,2),zb(1,2))
250	CALL multipl(ngrid,zcdh,zb0,zb)
251	DO ig=1,ngrid
252	z1(ig)=1./(za(ig,nlay)+zb(ig,nlay))
253	zc(ig,nlay)=za(ig,nlay)zh(ig,nlay)z1(ig)
254	zd(ig,nlay)=zb(ig,nlay)*z1(ig)
255	ENDDO
256
257	DO ilay=nlay-1,1,-1
258	DO ig=1,ngrid
259	z1(ig)=1./(za(ig,ilay)+zb(ig,ilay)+
260	$ zb(ig,ilay+1)*(1.-zd(ig,ilay+1)))
261	zc(ig,ilay)=(za(ig,ilay)*zh(ig,ilay)+
262	$ zb(ig,ilay+1)zc(ig,ilay+1))z1(ig)
263	zd(ig,ilay)=zb(ig,ilay)*z1(ig)
264	ENDDO
265	ENDDO
266
267	c-----------------------------------------------------------------------
268	c rajout eventuel de planck dans le shema implicite:
269	c --------------------------------------------------
270
271	z4st=4.5.67e-8ptimestep
272	c z4st=0.
273	DO ig=1,ngrid
274	zdplanck(ig)=z4stpemis(ig)ptsrf(ig)ptsrf(ig)ptsrf(ig)
275	ENDDO
276
277	c-----------------------------------------------------------------------
278	c calcul le l'evolution de la temperature du sol':
279	c -----------------------------------------------
280
281	DO ig=1,ngrid
282	z1(ig)=pcapcal(ig)ptsrf(ig)+cppzb(ig,1)*zc(ig,1)
283	s +zdplanck(ig)ptsrf(ig)+ pfluxsrf(ig)ptimestep
284	z2(ig)= pcapcal(ig)+cppzb(ig,1)(1.-zd(ig,1))+zdplanck(ig)
285	ztsrf2(ig)=z1(ig)/z2(ig)
286	zh(ig,1)=zc(ig,1)+zd(ig,1)*ztsrf2(ig)
287	pdtsrf(ig)=(ztsrf2(ig)-ptsrf(ig))/ptimestep
288	ENDDO
289
290	c-----------------------------------------------------------------------
291	c integration verticale finale:
292	c -----------------------------
293
294	DO ilay=2,nlay
295	DO ig=1,ngrid
296	zh(ig,ilay)=zc(ig,ilay)+zd(ig,ilay)*zh(ig,ilay-1)
297	ENDDO
298	ENDDO
299
300	c-----------------------------------------------------------------------
301	c calcul final des tendances de la diffusion verticale:
302	c -----------------------------------------------------
303
304	DO ilev = 1, nlay
305	DO ig=1,ngrid
306	pdudif(ig,ilev)=( zu(ig,ilev)-
307	$ (pu(ig,ilev)+pdufi(ig,ilev)*ptimestep) )/ptimestep
308	pdvdif(ig,ilev)=( zv(ig,ilev)-
309	$ (pv(ig,ilev)+pdvfi(ig,ilev)*ptimestep) )/ptimestep
310	pdhdif(ig,ilev)=( zh(ig,ilev)-
311	$ (ph(ig,ilev)+pdhfi(ig,ilev)*ptimestep) )/ptimestep
312	ENDDO
313	ENDDO
314
315	IF(lwrite) THEN
316	PRINT*
317	PRINT*,'Diagnostique de la diffusion verticale'
318	PRINT*,'h avant et apres diffusion verticale'
319	PRINT*,ptsrf(ngrid/2+1),ztsrf2(ngrid/2+1)
320	DO 3110 ilev=1,nlay
321	PRINT*,ph(ngrid/2+1,ilev),zh(ngrid/2+1,ilev)
322	3110 CONTINUE
323	ENDIF
324	c---------------------------------------------------------------------
325	RETURN
326	END SUBROUTINE vdif
327	!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
328
329	SUBROUTINE convadj(ngrid,nlay,ptimestep,
330	S pplay,pplev,ppopsk,
331	$ pu,pv,ph,
332	$ pdufi,pdvfi,pdhfi,
333	$ pduadj,pdvadj,pdhadj)
334	IMPLICIT NONE
335
336	c=======================================================================
337	c
338	c ajustement convectif sec
339	c on peut ajouter les tendances pdhfi au profil pdh avant l'ajustement
340	c'
341	c=======================================================================
342
343	c-----------------------------------------------------------------------
344	c declarations:
345	c -------------
346	c arguments:
347	c ----------
348
349	INTEGER ngrid,nlay
350	REAL ptimestep
351	REAL ph(ngrid,nlay),pdhfi(ngrid,nlay),pdhadj(ngrid,nlay)
352	REAL pplay(ngrid,nlay),pplev(ngrid,nlay+1),ppopsk(ngrid,nlay)
353	REAL pu(ngrid,nlay),pdufi(ngrid,nlay),pduadj(ngrid,nlay)
354	REAL pv(ngrid,nlay),pdvfi(ngrid,nlay),pdvadj(ngrid,nlay)
355
356	c local:
357	c ------
358
359	INTEGER ig,i,l,l1,l2,jj
360	INTEGER jcnt, jadrs(ngrid)
361
362	REAL*8 sig(nlayermx+1),sdsig(nlayermx),dsig(nlayermx)
363	REAL*8 zu(ngrid,nlayermx),zv(ngrid,nlayermx)
364	REAL*8 zh(ngrid,nlayermx)
365	REAL*8 zu2(ngrid,nlayermx),zv2(ngrid,nlayermx)
366	REAL*8 zh2(ngrid,nlayermx)
367	REAL*8 zhm,zsm,zum,zvm,zalpha
368
369	LOGICAL vtest(ngrid),down
370
371	c
372	c-----------------------------------------------------------------------
373	c initialisation:
374	c ---------------
375	c
376	IF(ngrid.NE.ngrid) THEN
377	PRINT*
378	PRINT*,'STOP dans convadj'
379	PRINT*,'ngrid =',ngrid
380	PRINT*,'ngrid =',ngrid
381	ENDIF
382	c
383	c-----------------------------------------------------------------------
384	c detection des profils a modifier:
385	c ---------------------------------
386	c si le profil est a modifier
387	c (i.e. ph(niv_sup) < ph(niv_inf) )
388	c alors le tableau "vtest" est mis a .TRUE. ;
389	c sinon, il reste a sa valeur initiale (.FALSE.)
390	c cette operation est vectorisable
391	c On en profite pour copier la valeur initiale de "ph"
392	c dans le champ de travail "zh"
393
394
395	DO 1010 l=1,nlay
396	DO 1015 ig=1,ngrid
397	zh(ig,l)=ph(ig,l)+pdhfi(ig,l)*ptimestep
398	zu(ig,l)=pu(ig,l)+pdufi(ig,l)*ptimestep
399	zv(ig,l)=pv(ig,l)+pdvfi(ig,l)*ptimestep
400	1015 CONTINUE
401	1010 CONTINUE
402
403	zu2(:,:)=zu(:,:)
404	zv2(:,:)=zv(:,:)
405	zh2(:,:)=zh(:,:)
406
407	DO 1020 ig=1,ngrid
408	vtest(ig)=.FALSE.
409	1020 CONTINUE
410	c
411	DO 1040 l=2,nlay
412	DO 1060 ig=1,ngrid
413	CRAY vtest(ig)=CVMGM(.TRUE. , vtest(ig),
414	CRAY . zh2(ig,l)-zh2(ig,l-1))
415	IF(zh2(ig,l).LT.zh2(ig,l-1)) vtest(ig)=.TRUE.
416	1060 CONTINUE
417	1040 CONTINUE
418	c
419	CRAY CALL WHENNE(ngrid, vtest, 1, 0, jadrs, jcnt)
420	jcnt=0
421	DO 1070 ig=1,ngrid
422	IF(vtest(ig)) THEN
423	jcnt=jcnt+1
424	jadrs(jcnt)=ig
425	ENDIF
426	1070 CONTINUE
427
428
429	c-----------------------------------------------------------------------
430	c Ajustement des "jcnt" profils instables indices par "jadrs":
431	c ------------------------------------------------------------
432	c
433	DO 1080 jj = 1, jcnt
434	c
435	i = jadrs(jj)
436	c
437	c Calcul des niveaux sigma sur cette colonne
438	DO l=1,nlay+1
439	sig(l)=pplev(i,l)/pplev(i,1)
440	ENDDO
441	DO l=1,nlay
442	dsig(l)=sig(l)-sig(l+1)
443	sdsig(l)=ppopsk(i,l)*dsig(l)
444	ENDDO
445	l2 = 1
446	c
447	c -- boucle de sondage vers le haut
448	c
449	cins$ Loop
450	8000 CONTINUE
451	c
452	l2 = l2 + 1
453	c
454	cins$ Exit
455	IF (l2 .GT. nlay) Goto 8001
456	c
457	IF (zh2(i, l2) .LT. zh2(i, l2-1)) THEN
458	c
459	c -- l2 est le niveau le plus haut de la colonne instable
460	c
461	l1 = l2 - 1
462	l = l1
463	zsm = sdsig(l2)
464	zhm = zh2(i, l2)
465	c
466	c -- boucle de sondage vers le bas
467	c
468	cins$ Loop
469	8020 CONTINUE
470	c
471	zsm = zsm + sdsig(l)
472	zhm = zhm + sdsig(l) * (zh2(i, l) - zhm) / zsm
473	c
474	c -- doit on etendre la colonne vers le bas ?
475	c
476	c_EC (M1875) 20/6/87 : AND -> AND THEN
477	c
478	down = .FALSE.
479	IF (l1 .NE. 1) THEN !-- and then
480	IF (zhm .LT. zh2(i, l1-1)) THEN
481	down = .TRUE.
482	END IF
483	END IF
484	c
485	IF (down) THEN
486	c
487	l1 = l1 - 1
488	l = l1
489	c
490	ELSE
491	c
492	c -- peut on etendre la colonne vers le haut ?
493	c
494	cins$ Exit
495	IF (l2 .EQ. nlay) Goto 8021
496	c
497	cins$ Exit
498	IF (zh2(i, l2+1) .GE. zhm) Goto 8021
499	c
500	l2 = l2 + 1
501	l = l2
502	c
503	END IF
504	c
505	cins$ End Loop
506	GO TO 8020
507	8021 CONTINUE
508	c
509	c -- nouveau profil : constant (valeur moyenne)
510	c
511	zalpha=0.
512	zum=0.
513	zvm=0.
514	DO 1100 l = l1, l2
515	zalpha=zalpha+ABS(zh2(i,l)-zhm)*dsig(l)
516	zh2(i, l) = zhm
517	zum=zum+dsig(l)*zu(i,l)
518	zvm=zvm+dsig(l)*zv(i,l)
519	1100 CONTINUE
520	zalpha=zalpha/(zhm*(sig(l1)-sig(l2+1)))
521	zum=zum/(sig(l1)-sig(l2+1))
522	zvm=zvm/(sig(l1)-sig(l2+1))
523	IF(zalpha.GT.1.) THEN
524	PRINT*,'WARNING dans convadj zalpha=',zalpha
525	if(ig.eq.1) then
526	print*,'Au pole nord'
527	elseif (ig.eq.ngrid) then
528	print*,'Au pole sud'
529	else
530	print*,'Point i=',
531	. ig-((ig-1)/iim)*iim,'j=',(ig-1)/iim+1
532	endif
533	! STOP !problem with icosa pole
534	zalpha=1.
535	ELSE
536	c IF(zalpha.LT.0.) STOP'zalpha=0'
537	IF(zalpha.LT.1.e-5) zalpha=1.e-5
538	ENDIF
539	DO l=l1,l2
540	zu2(i,l)=zu2(i,l)+zalpha*(zum-zu2(i,l))
541	zv2(i,l)=zv2(i,l)+zalpha*(zvm-zv2(i,l))
542	ENDDO
543
544	l2 = l2 + 1
545	c
546	END IF
547	c
548	cins$ End Loop
549	GO TO 8000
550	8001 CONTINUE
551	c
552	1080 CONTINUE
553	c
554	DO 4000 l=1,nlay
555	DO 4020 ig=1,ngrid
556	pdhadj(ig,l)=(zh2(ig,l)-zh(ig,l))/ptimestep
557	pduadj(ig,l)=(zu2(ig,l)-zu(ig,l))/ptimestep
558	pdvadj(ig,l)=(zv2(ig,l)-zv(ig,l))/ptimestep
559	4020 CONTINUE
560	4000 CONTINUE
561	c
562	RETURN
563	END SUBROUTINE convadj
564	!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
565
566	SUBROUTINE soil(ngrid,nsoil,firstcall,ptherm_i,
567	s ptimestep,ptsrf,ptsoil,
568	s pcapcal,pfluxgrd)
569	IMPLICIT NONE
570
571	c=======================================================================
572	c
573	c Auteur: Frederic Hourdin 30/01/92
574	c -------
575	c
576	c objet: computation of : the soil temperature evolution
577	c ------ the surfacic heat capacity "Capcal"
578	c the surface conduction flux pcapcal
579	c
580	c
581	c Method: implicit time integration
582	c -------
583	c Consecutive ground temperatures are related by:
584	c T(k+1) = C(k) + D(k)*T(k) (1)
585	c the coefficients C and D are computed at the t-dt time-step.
586	c Routine structure:
587	c 1)new temperatures are computed using (1)
588	c 2)C and D coefficients are computed from the new temperature
589	c profile for the t+dt time-step
590	c 3)the coefficients A and B are computed where the diffusive
591	c fluxes at the t+dt time-step is given by
592	c Fdiff = A + B Ts(t+dt)
593	c or Fdiff = F0 + Capcal (Ts(t+dt)-Ts(t))/dt
594	c with F0 = A + B (Ts(t))
595	c Capcal = B*dt
596	c
597	c Interface:
598	c ----------
599	c
600	c Arguments:
601	c ----------
602	c ngird number of grid-points
603	c ptimestep physical timestep (s)
604	c pto(ngrid,nsoil) temperature at time-step t (K)
605	c ptn(ngrid,nsoil) temperature at time step t+dt (K)
606	c pcapcal(ngrid) specific heat (Wm-2s*K-1)
607	c pfluxgrd(ngrid) surface diffusive flux from ground (Wm-2)
608	c
609	c=======================================================================
610	c declarations:
611	c -------------
612	c-----------------------------------------------------------------------
613	c arguments
614	c ---------
615
616	INTEGER ngrid,nsoil
617	REAL ptimestep
618	REAL ptsrf(ngrid),ptsoil(ngrid,nsoilmx),ptherm_i(ngrid)
619	REAL pcapcal(ngrid),pfluxgrd(ngrid)
620	LOGICAL firstcall
621
622	c-----------------------------------------------------------------------
623	c local arrays
624	c ------------
625
626	INTEGER ig,jk
627	REAL za(ngrid),zb(ngrid)
628	REAL zdz2(nsoilmx),z1(ngrid)
629	REAL min_period,dalph_soil
630
631	c local saved variables:
632	c ----------------------
633	REAL dz1(nsoilmx),dz2(nsoilmx)
634	REAL zc(ngrid,nsoilmx),zd(ngrid,nsoilmx)
635	REAL lambda
636
637	!!!!!!!! SARVESH !!!!!!! SAVE ATTRIBUTE
638	!! SAVE dz1,dz2,zc,zd,lambda
639
640	c-----------------------------------------------------------------------
641	c Depthts:
642	c --------
643
644	REAL fz,rk,fz1,rk1,rk2
645	fz(rk)=fz1(dalph_soil*rk-1.)/(dalph_soil-1.)
646
647	IF (firstcall) THEN
648
649	c-----------------------------------------------------------------------
650	c ground levels
651	c grnd=z/l where l is the skin depth of the diurnal cycle:
652	c --------------------------------------------------------
653
654	min_period=20000.
655	dalph_soil=2.
656
657	OPEN(99,file='soil.def',status='old',form='formatted',err=9999)
658	READ(99,*) min_period
659	READ(99,*) dalph_soil
660	PRINT*,'Discretization for the soil model'
661	PRINT*,'First level e-folding depth',min_period,
662	s ' dalph',dalph_soil
663	CLOSE(99)
664	9999 CONTINUE
665
666	c la premiere couche represente un dixieme de cycle diurne
667	fz1=sqrt(min_period/3.14)
668
669	DO jk=1,nsoil
670	rk1=jk
671	rk2=jk-1
672	dz2(jk)=fz(rk1)-fz(rk2)
673	ENDDO
674	DO jk=1,nsoil-1
675	rk1=jk+.5
676	rk2=jk-.5
677	dz1(jk)=1./(fz(rk1)-fz(rk2))
678	ENDDO
679	lambda=fz(.5)*dz1(1)
680	PRINT*,'full layers, intermediate layers (secoonds)'
681	DO jk=1,nsoil
682	rk=jk
683	rk1=jk+.5
684	rk2=jk-.5
685	PRINT,fz(rk1)fz(rk2)*3.14,
686	s fz(rk)fz(rk)3.14
687	ENDDO
688
689	c Initialisations:
690	c ----------------
691
692	ELSE
693	c-----------------------------------------------------------------------
694	c Computation of the soil temperatures using the Cgrd and Dgrd
695	c coefficient computed at the previous time-step:
696	c -----------------------------------------------
697
698	c surface temperature
699	DO ig=1,ngrid
700	ptsoil(ig,1)=(lambda*zc(ig,1)+ptsrf(ig))/
701	s (lambda*(1.-zd(ig,1))+1.)
702	ENDDO
703
704	c other temperatures
705	DO jk=1,nsoil-1
706	DO ig=1,ngrid
707	ptsoil(ig,jk+1)=zc(ig,jk)+zd(ig,jk)*ptsoil(ig,jk)
708	ENDDO
709	ENDDO
710
711	ENDIF
712	c-----------------------------------------------------------------------
713	c Computation of the Cgrd and Dgrd coefficient for the next step:
714	c ---------------------------------------------------------------
715
716	DO jk=1,nsoil
717	zdz2(jk)=dz2(jk)/ptimestep
718	ENDDO
719
720	DO ig=1,ngrid
721	z1(ig)=zdz2(nsoil)+dz1(nsoil-1)
722	zc(ig,nsoil-1)=zdz2(nsoil)*ptsoil(ig,nsoil)/z1(ig)
723	zd(ig,nsoil-1)=dz1(nsoil-1)/z1(ig)
724	ENDDO
725
726	DO jk=nsoil-1,2,-1
727	DO ig=1,ngrid
728	z1(ig)=1./(zdz2(jk)+dz1(jk-1)+dz1(jk)*(1.-zd(ig,jk)))
729	zc(ig,jk-1)=
730	s (ptsoil(ig,jk)zdz2(jk)+dz1(jk)zc(ig,jk))*z1(ig)
731	zd(ig,jk-1)=dz1(jk-1)*z1(ig)
732	ENDDO
733	ENDDO
734
735	c-----------------------------------------------------------------------
736	c computation of the surface diffusive flux from ground and
737	c calorific capacity of the ground:
738	c ---------------------------------
739
740	DO ig=1,ngrid
741	pfluxgrd(ig)=ptherm_i(ig)dz1(1)
742	s (zc(ig,1)+(zd(ig,1)-1.)*ptsoil(ig,1))
743	pcapcal(ig)=ptherm_i(ig)*
744	s (dz2(1)+ptimestep(1.-zd(ig,1))dz1(1))
745	z1(ig)=lambda*(1.-zd(ig,1))+1.
746	pcapcal(ig)=pcapcal(ig)/z1(ig)
747	pfluxgrd(ig)=pfluxgrd(ig)
748	s +pcapcal(ig)(ptsoil(ig,1)z1(ig)-lambda*zc(ig,1)-ptsrf(ig))
749	s /ptimestep
750	ENDDO
751
752	RETURN
753	END SUBROUTINE SOIL
754	!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
755
756	SUBROUTINE vdif_cd( ngrid,pz0,pg,pz,pu,pv,pts,ph,pcdv,pcdh)
757	IMPLICIT NONE
758	c=======================================================================
759	c
760	c Subject: computation of the surface drag coefficient using the
761	c ------- approch developed by Loui for ECMWF.
762	c
763	c Author: Frederic Hourdin 15 /10 /93
764	c -------
765	c
766	c Arguments:
767	c ----------
768	c
769	c inputs:
770	c ------
771	c ngrid size of the horizontal grid
772	c pg gravity (m s -2)
773	c pz(ngrid) height of the first atmospheric layer
774	c pu(ngrid) u component of the wind in that layer
775	c pv(ngrid) v component of the wind in that layer
776	c pts(ngrid) surfacte temperature
777	c ph(ngrid) potential temperature T*(p/ps)^kappa
778	c
779	c outputs:
780	c --------
781	c pcdv(ngrid) Cd for the wind
782	c pcdh(ngrid) Cd for potential temperature
783	c
784	c=======================================================================
785	c
786	c-----------------------------------------------------------------------
787	c Declarations:
788	c -------------
789
790	c Arguments:
791	c ----------
792
793	INTEGER ngrid,nlay
794	REAL pz0(ngrid)
795	REAL pg,pz(ngrid)
796	REAL pu(ngrid),pv(ngrid)
797	REAL pts(ngrid),ph(ngrid)
798	REAL pcdv(ngrid),pcdh(ngrid)
799
800	c Local:
801	c ------
802
803	INTEGER ig
804
805	REAL zu2,z1,zri,zcd0,zz
806
807	REAL karman,b,c,d,c2b,c3bc,c3b,z0,umin2
808	LOGICAL firstcal
809	DATA karman,b,c,d,umin2/.4,5.,5.,5.,1.e-12/
810	DATA firstcal/.true./
811	SAVE b,c,d,karman,c2b,c3bc,c3b,firstcal,umin2
812
813	c-----------------------------------------------------------------------
814	c couche de surface:
815	c ------------------
816
817	c DO ig=1,ngrid
818	c zu2=pu(ig)pu(ig)+pv(ig)pv(ig)+umin2
819	c pcdv(ig)=pz0(ig)*(1.+sqrt(zu2))
820	c pcdh(ig)=pcdv(ig)
821	c ENDDO
822	c RETURN
823
824	IF (firstcal) THEN
825	c2b=2.*b
826	c3bc=3.bc
827	c3b=3.*b
828	firstcal=.false.
829	ENDIF
830
831	c!!!! WARNING, verifier la formule originale de Louis!
832	DO ig=1,ngrid
833	zu2=pu(ig)pu(ig)+pv(ig)pv(ig)+umin2
834	zri=pgpz(ig)(ph(ig)-pts(ig))/(ph(ig)*zu2)
835	z1=1.+pz(ig)/pz0(ig)
836	zcd0=karman/log(z1)
837	zcd0=zcd0zcd0sqrt(zu2)
838	IF(zri.LT.0.) THEN
839	z1=bzri/(1.+c3bczcd0sqrt(-z1zri))
840	pcdv(ig)=zcd0(1.-2.z1)
841	pcdh(ig)=zcd0(1.-3.z1)
842	ELSE
843	zz=sqrt(1.+d*zri)
844	pcdv(ig)=zcd0/(1.+c2b*zri/zz)
845	pcdh(ig)=zcd0/(1.+c3bzrizz)
846	ENDIF
847	ENDDO
848
849	c-----------------------------------------------------------------------
850
851	RETURN
852	END SUBROUTINE vdif_cd
853	!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
854
855	SUBROUTINE vdif_k(ngrid,nlay,
856	s ptimestep,pg,pzlev,pzlay,pz0,pu,pv,ph,pcdv,pkv,pkh,pq2,pq2l)
857
858	IMPLICIT NONE
859
860	INTEGER ngrid,nlay
861
862	REAL ptimestep
863	REAL pzlay(ngrid,nlay),pzlev(ngrid,nlay+1)
864	REAL pz0(ngrid)
865	REAL pu(ngrid,nlay),pv(ngrid,nlay),ph(ngrid,nlay)
866	REAL pg,pcdv(ngrid)
867	REAL pkv(ngrid,nlay+1),pkh(ngrid,nlay+1)
868	REAL pq2(ngrid,nlay+1),pq2l(ngrid,nlay+1) !!!! SARVESH ADDED to
869
870	INTEGER ig,il
871	REAL zdu,zdv,zri,zdvodz2,zdz,z1,lmix
872
873	REAL karman
874	SAVE karman
875	! DATA lmixmin,emin_turb,karman/100.,1.e-8,.4/
876	!!!!! SARVESH !!!!!!
877	!Error: Host associated variable 'lmixmin' may not be in the DATA statement
878
879	! print*,'LMIXMIN',lmixmin
880	DO ig=1,ngrid
881	pkv(ig,1)=0.
882	pkh(ig,1)=0.
883	pkv(ig,nlay+1)=0.
884	pkh(ig,nlay+1)=0.
885	ENDDO
886	c s ' zdu,zdv,zdz,zdovdz2,ph(ig,il)+ph(ig,il-1)'
887	DO il=2,nlay
888	DO ig=1,ngrid
889	z1=pzlev(ig,il)+pz0(ig)
890	lmix=karmanz1/(1.+karmanz1/lmixmin)
891	c lmix=lmixmin
892	c WARNING test lmix=lmixmin
893	zdu=pu(ig,il)-pu(ig,il-1)
894	zdv=pv(ig,il)-pv(ig,il-1)
895	zdz=pzlay(ig,il)-pzlay(ig,il-1)
896	zdvodz2=(zduzdu+zdvzdv)/(zdz*zdz)
897	IF(zdvodz2.LT.1.e-5) THEN
898	pkv(ig,il)=lmix*sqrt(emin_turb)
899	ELSE
900	zri=2.pg(ph(ig,il)-ph(ig,il-1))
901	s / (zdz* (ph(ig,il)+ph(ig,il-1)) *zdvodz2 )
902	pkv(ig,il)=
903	s lmixsqrt(MAX(lmixlmixzdvodz2(1-zri/.4),emin_turb))
904	ENDIF
905	pkh(ig,il)=pkv(ig,il)
906	c IF(ig.EQ.ngrid/2+1) PRINT*,il,lmix,pkv(ig,il),
907	c s zdu,zdv,zdz,zdvodz2,ph(ig,il)+ph(ig,il-1),
908	c s lmixlmixzdvodz2*(1-zri/.4),emin_turb,zri,ph(ig,il)-ph(ig,il-1),
909	c s ph(ig,il),ph(ig,il-1)
910	ENDDO
911	ENDDO
912
913	RETURN
914	END SUBROUTINE vdif_k
915	!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
916
917	SUBROUTINE multipl(n,x1,x2,y)
918	IMPLICIT NONE
919	c====================================================================
920	c
921	c multiplication de deux vecteurs
922	c
923	c=======================================================================
924	c
925	INTEGER n,i
926	REAL x1(n),x2(n),y(n)
927	c
928	DO 10 i=1,n
929	y(i)=x1(i)*x2(i)
930	10 CONTINUE
931	c
932	RETURN
933	END SUBROUTINE multipl
934	!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
935	END MODULE SURFACE_PROCESS

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

Context Navigation

source: codes/icosagcm/trunk/src/surface_mod.F @ 149

Download in other formats: