Context Navigation

1DUTILS.h_no_writelim @ 1691

Last change on this file since 1691 was 1608, checked in by lguez, 13 years ago

In "phy1d/*", replaced obsolete calls to "lnblnk" by calls to "trim".

In "lmdz1d.F", use collector module "ioipsl" rather than specific
module "calendar", which is an internal module of IOIPSL. Removed line
containing only "#" (causes compilation error). Bug fix in call to
"init_phys_lmdz".

File size: 75.6 KB

Line
1	!
2	! $Id: conf_unicol.F 1279 2010-08-04 17:20:56Z lahellec $
3	!
4	c
5	c
6	SUBROUTINE conf_unicol( tapedef )
7	c
8	#ifdef CPP_IOIPSL
9	use IOIPSL
10	#else
11	! if not using IOIPSL, we still need to use (a local version of) getin
12	use ioipsl_getincom
13	#endif
14	IMPLICIT NONE
15	c-----------------------------------------------------------------------
16	c Auteurs : A. Lahellec .
17	c
18	c Arguments :
19	c
20	c tapedef :
21
22	INTEGER tapedef
23	c
24	c Declarations :
25	c --------------
26
27	#include "compar1d.h"
28	#include "flux_arp.h"
29	#include "fcg_gcssold.h"
30	#include "fcg_racmo.h"
31	#include "iniprint.h"
32	c
33	c
34	c local:
35	c ------
36
37	c CHARACTER ch172,ch272,ch372,ch412
38
39	c
40	c -------------------------------------------------------------------
41	c
42	c ......... Initilisation parametres du lmdz1D ..........
43	c
44	c---------------------------------------------------------------------
45	c initialisations:
46	c ----------------
47
48	!Config Key = lunout
49	!Config Desc = unite de fichier pour les impressions
50	!Config Def = 6
51	!Config Help = unite de fichier pour les impressions
52	!Config (defaut sortie standard = 6)
53	lunout=6
54	CALL getin('lunout', lunout)
55	IF (lunout /= 5 .and. lunout /= 6) THEN
56	OPEN(lunout,FILE='lmdz.out')
57	ENDIF
58
59	!Config Key = prt_level
60	!Config Desc = niveau d'impressions de débogage
61	!Config Def = 0
62	!Config Help = Niveau d'impression pour le débogage
63	!Config (0 = minimum d'impression)
64	c prt_level = 0
65	c CALL getin('prt_level',prt_level)
66
67	c-----------------------------------------------------------------------
68	c Parametres de controle du run:
69	c-----------------------------------------------------------------------
70
71	!Config Key = restart
72	!Config Desc = on repart des startphy et start1dyn
73	!Config Def = false
74	!Config Help = les fichiers restart doivent etre renomme en start
75	restart = .FALSE.
76	CALL getin('restart',restart)
77
78	!Config Key = forcing_type
79	!Config Desc = defines the way the SCM is forced:
80	!Config Def = 0
81	!!Config Help = 0 ==> forcing_les = .true.
82	! initial profiles from file prof.inp.001
83	! no forcing by LS convergence ;
84	! surface temperature imposed ;
85	! radiative cooling may be imposed (iflag_radia=0 in physiq.def)
86	! = 1 ==> forcing_radconv = .true.
87	! idem forcing_type = 0, but the imposed radiative cooling
88	! is set to 0 (hence, if iflag_radia=0 in physiq.def,
89	! then there is no radiative cooling at all)
90	! = 2 ==> forcing_toga = .true.
91	! initial profiles from TOGA-COARE IFA files
92	! LS convergence and SST imposed from TOGA-COARE IFA files
93	! = 3 ==> forcing_GCM2SCM = .true.
94	! initial profiles from the GCM output
95	! LS convergence imposed from the GCM output
96	! = 4 ==> forcing_twpi = .true.
97	! initial profiles from TWPICE nc files
98	! LS convergence and SST imposed from TWPICE nc files
99	! = 5 ==> forcing_rico = .true.
100	! initial profiles from RICO idealized
101	! LS convergence imposed from RICO (cst)
102	! = 40 ==> forcing_GCSSold = .true.
103	! initial profile from GCSS file
104	! LS convergence imposed from GCSS file
105	!
106	forcing_type = 0
107	CALL getin('forcing_type',forcing_type)
108	imp_fcg_gcssold = .false.
109	ts_fcg_gcssold = .false.
110	Tp_fcg_gcssold = .false.
111	Tp_ini_gcssold = .false.
112	xTurb_fcg_gcssold = .false.
113	IF (forcing_type .eq.40) THEN
114	CALL getin('imp_fcg',imp_fcg_gcssold)
115	CALL getin('ts_fcg',ts_fcg_gcssold)
116	CALL getin('tp_fcg',Tp_fcg_gcssold)
117	CALL getin('tp_ini',Tp_ini_gcssold)
118	CALL getin('turb_fcg',xTurb_fcg_gcssold)
119	ENDIF
120
121	!Config Key = ok_flux_surf
122	!Config Desc = forcage ou non par les flux de surface
123	!Config Def = false
124	!Config Help = forcage ou non par les flux de surface
125	ok_flux_surf = .FALSE.
126	CALL getin('ok_flux_surf',ok_flux_surf)
127
128	!Config Key = time_ini
129	!Config Desc = meaningless in this case
130	!Config Def = 0.
131	!Config Help =
132	tsurf = 0.
133	CALL getin('time_ini',time_ini)
134
135	!Config Key = rlat et rlon
136	!Config Desc = latitude et longitude
137	!Config Def = 0.0 0.0
138	!Config Help = fixe la position de la colonne
139	xlat = 0.
140	xlon = 0.
141	CALL getin('rlat',xlat)
142	CALL getin('rlon',xlon)
143
144	!Config Key = airephy
145	!Config Desc = Grid cell area
146	!Config Def = 1.e11
147	!Config Help =
148	airefi = 1.e11
149	CALL getin('airephy',airefi)
150
151	!Config Key = nat_surf
152	!Config Desc = surface type
153	!Config Def = 0 (ocean)
154	!Config Help = 0=ocean,1=land,2=glacier,3=banquise
155	nat_surf = 0.
156	CALL getin('nat_surf',nat_surf)
157
158	!Config Key = tsurf
159	!Config Desc = surface temperature
160	!Config Def = 290.
161	!Config Help = not used if type_ts_forcing=1 in lmdz1d.F
162	tsurf = 290.
163	CALL getin('tsurf',tsurf)
164
165	!Config Key = psurf
166	!Config Desc = surface pressure
167	!Config Def = 102400.
168	!Config Help =
169	psurf = 102400.
170	CALL getin('psurf',psurf)
171
172	!Config Key = zsurf
173	!Config Desc = surface altitude
174	!Config Def = 0.
175	!Config Help =
176	zsurf = 0.
177	CALL getin('zsurf',zsurf)
178
179	!Config Key = rugos
180	!Config Desc = coefficient de frottement
181	!Config Def = 0.0001
182	!Config Help = calcul du Cdrag
183	rugos = 0.0001
184	CALL getin('rugos',rugos)
185
186	!Config Key = wtsurf et wqsurf
187	!Config Desc = ???
188	!Config Def = 0.0 0.0
189	!Config Help =
190	wtsurf = 0.0
191	wqsurf = 0.0
192	CALL getin('wtsurf',wtsurf)
193	CALL getin('wqsurf',wqsurf)
194
195	!Config Key = albedo
196	!Config Desc = albedo
197	!Config Def = 0.09
198	!Config Help =
199	albedo = 0.09
200	CALL getin('albedo',albedo)
201
202	!Config Key = agesno
203	!Config Desc = age de la neige
204	!Config Def = 30.0
205	!Config Help =
206	xagesno = 30.0
207	CALL getin('agesno',xagesno)
208
209	!Config Key = restart_runoff
210	!Config Desc = age de la neige
211	!Config Def = 30.0
212	!Config Help =
213	restart_runoff = 0.0
214	CALL getin('restart_runoff',restart_runoff)
215
216	!Config Key = qsolinp
217	!Config Desc = initial bucket water content (kg/m2) when land (5std)
218	!Config Def = 30.0
219	!Config Help =
220	qsolinp = 1.
221	CALL getin('qsolinp',qsolinp)
222
223	!Config Key = zpicinp
224	!Config Desc = denivellation orographie
225	!Config Def = 300.
226	!Config Help = input brise
227	zpicinp = 300.
228	CALL getin('zpicinp',zpicinp)
229
230
231	write(lunout,*)' +++++++++++++++++++++++++++++++++++++++'
232	write(lunout,*)' Configuration des parametres du gcm1D: '
233	write(lunout,*)' +++++++++++++++++++++++++++++++++++++++'
234	write(lunout,*)' restart = ', restart
235	write(lunout,*)' forcing_type = ', forcing_type
236	write(lunout,*)' time_ini = ', time_ini
237	write(lunout,*)' rlat = ', xlat
238	write(lunout,*)' rlon = ', xlon
239	write(lunout,*)' airephy = ', airefi
240	write(lunout,*)' nat_surf = ', nat_surf
241	write(lunout,*)' tsurf = ', tsurf
242	write(lunout,*)' psurf = ', psurf
243	write(lunout,*)' zsurf = ', zsurf
244	write(lunout,*)' rugos = ', rugos
245	write(lunout,*)' wtsurf = ', wtsurf
246	write(lunout,*)' wqsurf = ', wqsurf
247	write(lunout,*)' albedo = ', albedo
248	write(lunout,*)' xagesno = ', xagesno
249	write(lunout,*)' restart_runoff = ', restart_runoff
250	write(lunout,*)' qsolinp = ', qsolinp
251	write(lunout,*)' zpicinp = ', zpicinp
252	IF (forcing_type .eq.40) THEN
253	write(lunout,*) '--- Forcing type GCSS Old --- with:'
254	write(lunout,*)'imp_fcg',imp_fcg_gcssold
255	write(lunout,*)'ts_fcg',ts_fcg_gcssold
256	write(lunout,*)'tp_fcg',Tp_fcg_gcssold
257	write(lunout,*)'tp_ini',Tp_ini_gcssold
258	write(lunout,*)'xturb_fcg',xTurb_fcg_gcssold
259	ENDIF
260
261	write(lunout,*)' +++++++++++++++++++++++++++++++++++++++'
262	write(lunout,*)
263	c
264	RETURN
265	END
266	!
267	! $Id: dyn1deta0.F 1279 2010/07/30 A Lahellec$
268	!
269	c
270	SUBROUTINE dyn1deta0(fichnom,plev,play,phi,phis,presnivs,
271	& ucov,vcov,temp,q,omega2)
272	USE dimphy
273	USE mod_grid_phy_lmdz
274	USE mod_phys_lmdz_para
275	USE iophy
276	USE phys_state_var_mod
277	USE iostart
278	USE write_field_phy
279	USE infotrac
280	use control_mod
281
282	IMPLICIT NONE
283	c=======================================================
284	c Ecriture du fichier de redemarrage sous format NetCDF
285	c=======================================================
286	c Declarations:
287	c -------------
288	#include "dimensions.h"
289	#include "comconst.h"
290	#include "temps.h"
291	!!#include "control.h"
292	#include "logic.h"
293	#include "netcdf.inc"
294
295	c Arguments:
296	c ----------
297	CHARACTER() fichnom
298	cAl1 plev tronque pour .nc mais plev(klev+1):=0
299	real :: plev(klon,klev),play (klon,klev),phi(klon,klev)
300	real :: presnivs(klon,klev)
301	real :: ucov(klon,klev),vcov(klon,klev),temp(klon,klev)
302	real :: q(klon,klev,nqtot),omega2(klon,klev)
303	c real :: ug(klev),vg(klev),fcoriolis
304	real :: phis(klon)
305
306	c Variables locales pour NetCDF:
307	c ------------------------------
308	INTEGER nid, nvarid
309	INTEGER idim_s
310	INTEGER ierr, ierr_file
311	INTEGER iq
312	INTEGER length
313	PARAMETER (length = 100)
314	REAL tab_cntrl(length) ! tableau des parametres du run
315	character*4 nmq(nqtot)
316	character*12 modname
317	character*80 abort_message
318	LOGICAL found
319	c
320	INTEGER nb
321
322	modname = 'dyn1deta0 : '
323	nmq(1)="vap"
324	nmq(2)="cond"
325	do iq=3,nqtot
326	write(nmq(iq),'("tra",i1)') iq-2
327	enddo
328	print*,'in dyn1deta0 ',fichnom,klon,klev,nqtot
329	CALL open_startphy(fichnom)
330	print*,'after open startphy ',fichnom,nmq
331
332	c
333	c Lecture des parametres de controle:
334	c
335	CALL get_var("controle",tab_cntrl)
336
337
338	im = tab_cntrl(1)
339	jm = tab_cntrl(2)
340	lllm = tab_cntrl(3)
341	day_ref = tab_cntrl(4)
342	annee_ref = tab_cntrl(5)
343	c rad = tab_cntrl(6)
344	c omeg = tab_cntrl(7)
345	c g = tab_cntrl(8)
346	c cpp = tab_cntrl(9)
347	c kappa = tab_cntrl(10)
348	c daysec = tab_cntrl(11)
349	c dtvr = tab_cntrl(12)
350	c etot0 = tab_cntrl(13)
351	c ptot0 = tab_cntrl(14)
352	c ztot0 = tab_cntrl(15)
353	c stot0 = tab_cntrl(16)
354	c ang0 = tab_cntrl(17)
355	c pa = tab_cntrl(18)
356	c preff = tab_cntrl(19)
357	c
358	c clon = tab_cntrl(20)
359	c clat = tab_cntrl(21)
360	c grossismx = tab_cntrl(22)
361	c grossismy = tab_cntrl(23)
362	c
363	IF ( tab_cntrl(24).EQ.1. ) THEN
364	fxyhypb = . TRUE .
365	c dzoomx = tab_cntrl(25)
366	c dzoomy = tab_cntrl(26)
367	c taux = tab_cntrl(28)
368	c tauy = tab_cntrl(29)
369	ELSE
370	fxyhypb = . FALSE .
371	ysinus = . FALSE .
372	IF( tab_cntrl(27).EQ.1. ) ysinus = . TRUE.
373	ENDIF
374
375	day_ini = tab_cntrl(30)
376	itau_dyn = tab_cntrl(31)
377	c .................................................................
378	c
379	c
380	c PRINT*,'rad,omeg,g,cpp,kappa',rad,omeg,g,cpp,kappa
381	cAl1
382	Print*,'day_ref,annee_ref,day_ini,itau_dyn',
383	& day_ref,annee_ref,day_ini,itau_dyn
384
385	c Lecture des champs
386	c
387	plev(1,klev+1)=0.
388	CALL get_field("plev",plev,found)
389	IF (.NOT. found) PRINT*, modname//'Le champ <Plev> est absent'
390	CALL get_field("play",play,found)
391	IF (.NOT. found) PRINT*, modname//'Le champ <Play> est absent'
392	CALL get_field("phi",phi,found)
393	IF (.NOT. found) PRINT*, modname//'Le champ <Phi> est absent'
394	CALL get_field("phis",phis,found)
395	IF (.NOT. found) PRINT*, modname//'Le champ <Phis> est absent'
396	CALL get_field("presnivs",presnivs,found)
397	IF (.NOT. found) PRINT*, modname//'Le champ <Presnivs> est absent'
398	CALL get_field("ucov",ucov,found)
399	IF (.NOT. found) PRINT*, modname//'Le champ <ucov> est absent'
400	CALL get_field("vcov",vcov,found)
401	IF (.NOT. found) PRINT*, modname//'Le champ <vcov> est absent'
402	CALL get_field("temp",temp,found)
403	IF (.NOT. found) PRINT*, modname//'Le champ <temp> est absent'
404	CALL get_field("omega2",omega2,found)
405	IF (.NOT. found) PRINT*, modname//'Le champ <omega2> est absent'
406
407	Do iq=1,nqtot
408	CALL get_field("q"//nmq(iq),q(:,:,iq),found)
409	IF (.NOT. found)
410	& PRINT*, modname//'Le champ <q'//nmq//'> est absent'
411	EndDo
412
413	CALL close_startphy
414	print*,' close startphy'
415	. ,fichnom,play(1,1),play(1,klev),temp(1,klev)
416	c
417	RETURN
418	END
419	!
420	! $Id: dyn1dredem.F 1279 2010/07/29 A Lahellec$
421	!
422	c
423	SUBROUTINE dyn1dredem(fichnom,plev,play,phi,phis,presnivs,
424	& ucov,vcov,temp,q,omega2)
425	USE dimphy
426	USE mod_grid_phy_lmdz
427	USE mod_phys_lmdz_para
428	USE phys_state_var_mod
429	USE iostart
430	USE infotrac
431	use control_mod
432
433	IMPLICIT NONE
434	c=======================================================
435	c Ecriture du fichier de redemarrage sous format NetCDF
436	c=======================================================
437	c Declarations:
438	c -------------
439	#include "dimensions.h"
440	#include "comconst.h"
441	#include "temps.h"
442	!!#include "control.h"
443	#include "logic.h"
444	#include "netcdf.inc"
445
446	c Arguments:
447	c ----------
448	CHARACTER() fichnom
449	REAL time
450	cAl1 plev tronque pour .nc mais plev(klev+1):=0
451	real :: plev(klon,klev),play (klon,klev),phi(klon,klev)
452	real :: presnivs(klon,klev)
453	real :: ucov(klon,klev),vcov(klon,klev),temp(klon,klev)
454	real :: q(klon,klev,nqtot)
455	real :: omega2(klon,klev),rho(klon,klev+1)
456	c real :: ug(klev),vg(klev),fcoriolis
457	real :: phis(klon)
458
459	c Variables locales pour NetCDF:
460	c ------------------------------
461	INTEGER nid, nvarid
462	INTEGER idim_s
463	INTEGER ierr, ierr_file
464	INTEGER iq,l
465	INTEGER length
466	PARAMETER (length = 100)
467	REAL tab_cntrl(length) ! tableau des parametres du run
468	character*4 nmq(nqtot)
469	character*20 modname
470	character*80 abort_message
471	c
472	INTEGER nb
473	SAVE nb
474	DATA nb / 0 /
475
476	REAL zan0,zjulian,hours
477	INTEGER yyears0,jjour0, mmois0
478	character*30 unites
479
480	cDbg
481	CALL open_restartphy(fichnom)
482	print*,'redm1 ',fichnom,klon,klev,nqtot
483	nmq(1)="vap"
484	nmq(2)="cond"
485	nmq(3)="tra1"
486	nmq(4)="tra2"
487
488	modname = 'dyn1dredem'
489	ierr = NF_OPEN(fichnom, NF_WRITE, nid)
490	IF (ierr .NE. NF_NOERR) THEN
491	PRINT*, "Pb. d ouverture "//fichnom
492	CALL abort
493	ENDIF
494
495	DO l=1,length
496	tab_cntrl(l) = 0.
497	ENDDO
498	tab_cntrl(1) = FLOAT(iim)
499	tab_cntrl(2) = FLOAT(jjm)
500	tab_cntrl(3) = FLOAT(llm)
501	tab_cntrl(4) = FLOAT(day_ref)
502	tab_cntrl(5) = FLOAT(annee_ref)
503	tab_cntrl(6) = rad
504	tab_cntrl(7) = omeg
505	tab_cntrl(8) = g
506	tab_cntrl(9) = cpp
507	tab_cntrl(10) = kappa
508	tab_cntrl(11) = daysec
509	tab_cntrl(12) = dtvr
510	c tab_cntrl(13) = etot0
511	c tab_cntrl(14) = ptot0
512	c tab_cntrl(15) = ztot0
513	c tab_cntrl(16) = stot0
514	c tab_cntrl(17) = ang0
515	c tab_cntrl(18) = pa
516	c tab_cntrl(19) = preff
517	c
518	c ..... parametres pour le zoom ......
519
520	c tab_cntrl(20) = clon
521	c tab_cntrl(21) = clat
522	c tab_cntrl(22) = grossismx
523	c tab_cntrl(23) = grossismy
524	c
525	IF ( fxyhypb ) THEN
526	tab_cntrl(24) = 1.
527	c tab_cntrl(25) = dzoomx
528	c tab_cntrl(26) = dzoomy
529	tab_cntrl(27) = 0.
530	c tab_cntrl(28) = taux
531	c tab_cntrl(29) = tauy
532	ELSE
533	tab_cntrl(24) = 0.
534	c tab_cntrl(25) = dzoomx
535	c tab_cntrl(26) = dzoomy
536	tab_cntrl(27) = 0.
537	tab_cntrl(28) = 0.
538	tab_cntrl(29) = 0.
539	IF( ysinus ) tab_cntrl(27) = 1.
540	ENDIF
541	CAl1 iday_end -> day_end
542	tab_cntrl(30) = FLOAT(day_end)
543	tab_cntrl(31) = FLOAT(itau_dyn + itaufin)
544	c
545	CALL put_var("controle","Param. de controle Dyn1D",tab_cntrl)
546	c
547
548	c Ecriture/extension de la coordonnee temps
549
550	nb = nb + 1
551
552	c Ecriture des champs
553	c
554	CALL put_field("plev","p interfaces sauf la nulle",plev)
555	CALL put_field("play","",play)
556	CALL put_field("phi","geopotentielle",phi)
557	CALL put_field("phis","geopotentiell de surface",phis)
558	CALL put_field("presnivs","",presnivs)
559	CALL put_field("ucov","",ucov)
560	CALL put_field("vcov","",vcov)
561	CALL put_field("temp","",temp)
562	CALL put_field("omega2","",omega2)
563
564	Do iq=1,nqtot
565	CALL put_field("q"//nmq(iq),"eau vap ou condens et traceurs",
566	. q(:,:,iq))
567	EndDo
568	CALL close_restartphy
569
570	c
571	RETURN
572	END
573	SUBROUTINE gr_fi_dyn(nfield,ngrid,im,jm,pfi,pdyn)
574	IMPLICIT NONE
575	!=======================================================================
576	! passage d'un champ de la grille scalaire a la grille physique
577	!=======================================================================
578
579	!-----------------------------------------------------------------------
580	! declarations:
581	! -------------
582
583	INTEGER im,jm,ngrid,nfield
584	REAL pdyn(im,jm,nfield)
585	REAL pfi(ngrid,nfield)
586
587	INTEGER i,j,ifield,ig
588
589	!-----------------------------------------------------------------------
590	! calcul:
591	! -------
592
593	DO ifield=1,nfield
594	! traitement des poles
595	DO i=1,im
596	pdyn(i,1,ifield)=pfi(1,ifield)
597	pdyn(i,jm,ifield)=pfi(ngrid,ifield)
598	ENDDO
599
600	! traitement des point normaux
601	DO j=2,jm-1
602	ig=2+(j-2)*(im-1)
603	CALL SCOPY(im-1,pfi(ig,ifield),1,pdyn(1,j,ifield),1)
604	pdyn(im,j,ifield)=pdyn(1,j,ifield)
605	ENDDO
606	ENDDO
607
608	RETURN
609	END
610
611
612
613	SUBROUTINE abort_gcm(modname, message, ierr)
614
615	USE IOIPSL
616	!
617	! Stops the simulation cleanly, closing files and printing various
618	! comments
619	!
620	! Input: modname = name of calling program
621	! message = stuff to print
622	! ierr = severity of situation ( = 0 normal )
623
624	character*20 modname
625	integer ierr
626	character*80 message
627
628	write(,) 'in abort_gcm'
629	call histclo
630	! call histclo(2)
631	! call histclo(3)
632	! call histclo(4)
633	! call histclo(5)
634	write(,) 'out of histclo'
635	write(,) 'Stopping in ', modname
636	write(,) 'Reason = ',message
637	call getin_dump
638	!
639	if (ierr .eq. 0) then
640	write(,) 'Everything is cool'
641	else
642	write(,) 'Houston, we have a problem ', ierr
643	endif
644	STOP
645	END
646	REAL FUNCTION fq_sat(kelvin, millibar)
647	!
648	IMPLICIT none
649	!======================================================================
650	! Autheur(s): Z.X. Li (LMD/CNRS)
651	! Objet: calculer la vapeur d'eau saturante (formule Centre Euro.)
652	!======================================================================
653	! Arguments:
654	! kelvin---input-R: temperature en Kelvin
655	! millibar--input-R: pression en mb
656	!
657	! fq_sat----output-R: vapeur d'eau saturante en kg/kg
658	!======================================================================
659	!
660	REAL kelvin, millibar
661	!
662	REAL r2es
663	PARAMETER (r2es=611.14 *18.0153/28.9644)
664	!
665	REAL r3les, r3ies, r3es
666	PARAMETER (R3LES=17.269)
667	PARAMETER (R3IES=21.875)
668	!
669	REAL r4les, r4ies, r4es
670	PARAMETER (R4LES=35.86)
671	PARAMETER (R4IES=7.66)
672	!
673	REAL rtt
674	PARAMETER (rtt=273.16)
675	!
676	REAL retv
677	PARAMETER (retv=28.9644/18.0153 - 1.0)
678	!
679	REAL zqsat
680	REAL temp, pres
681	! ------------------------------------------------------------------
682	!
683	!
684	temp = kelvin
685	pres = millibar * 100.0
686	! write(,)'kelvin,millibar=',kelvin,millibar
687	! write(,)'temp,pres=',temp,pres
688	!
689	IF (temp .LE. rtt) THEN
690	r3es = r3ies
691	r4es = r4ies
692	ELSE
693	r3es = r3les
694	r4es = r4les
695	ENDIF
696	!
697	zqsat=r2es/pres * EXP ( r3es*(temp-rtt) / (temp-r4es) )
698	zqsat=MIN(0.5,ZQSAT)
699	zqsat=zqsat/(1.-retv *zqsat)
700	!
701	fq_sat = zqsat
702	!
703	RETURN
704	END
705	subroutine scopy(n,sx,incx,sy,incy)
706	!
707	IMPLICIT NONE
708	!
709	integer n,incx,incy,ix,iy,i
710	real sx((n-1)incx+1),sy((n-1)incy+1)
711	!
712	iy=1
713	ix=1
714	do 10 i=1,n
715	sy(iy)=sx(ix)
716	ix=ix+incx
717	iy=iy+incy
718	10 continue
719	!
720	return
721	end
722	subroutine wrgradsfi(if,nl,field,name,titlevar)
723	implicit none
724
725	! Declarations
726
727	#include "gradsdef.h"
728
729	! arguments
730	integer if,nl
731	real field(imxjmxlmx)
732	character*10 name,file
733	character*10 titlevar
734
735	! local
736
737	integer im,jm,lm,i,j,l,iv,iii,iji,iif,ijf
738
739	logical writectl
740
741
742	writectl=.false.
743
744	! print*,if,iid(if),jid(if),ifd(if),jfd(if)
745	iii=iid(if)
746	iji=jid(if)
747	iif=ifd(if)
748	ijf=jfd(if)
749	im=iif-iii+1
750	jm=ijf-iji+1
751	lm=lmd(if)
752
753
754	! print*,'im,jm,lm,name,firsttime(if)'
755	! print*,im,jm,lm,name,firsttime(if)
756
757	if(firsttime(if)) then
758	if(name.eq.var(1,if)) then
759	firsttime(if)=.false.
760	ivar(if)=1
761	print*,'fin de l initialiation de l ecriture du fichier'
762	print*,file
763	print*,'fichier no: ',if
764	print*,'unit ',unit(if)
765	print*,'nvar ',nvar(if)
766	print*,'vars ',(var(iv,if),iv=1,nvar(if))
767	else
768	ivar(if)=ivar(if)+1
769	nvar(if)=ivar(if)
770	var(ivar(if),if)=name
771	tvar(ivar(if),if)=trim(titlevar)
772	nld(ivar(if),if)=nl
773	print*,'initialisation ecriture de ',var(ivar(if),if)
774	print*,'if ivar(if) nld ',if,ivar(if),nld(ivar(if),if)
775	endif
776	writectl=.true.
777	itime(if)=1
778	else
779	ivar(if)=mod(ivar(if),nvar(if))+1
780	if (ivar(if).eq.nvar(if)) then
781	writectl=.true.
782	itime(if)=itime(if)+1
783	endif
784
785	if(var(ivar(if),if).ne.name) then
786	print*,'Il faut stoker la meme succession de champs a chaque'
787	print*,'pas de temps'
788	print*,'fichier no: ',if
789	print*,'unit ',unit(if)
790	print*,'nvar ',nvar(if)
791	print*,'vars ',(var(iv,if),iv=1,nvar(if))
792
793	stop
794	endif
795	endif
796
797	! print*,'ivar(if),nvar(if),var(ivar(if),if),writectl'
798	! print*,ivar(if),nvar(if),var(ivar(if),if),writectl
799	do l=1,nl
800	irec(if)=irec(if)+1
801	! print*,'Ecrit rec=',irec(if),iii,iif,iji,ijf,
802	! s (l-1)imd(if)jmd(if)+(iji-1)*imd(if)+iii
803	! s ,(l-1)imd(if)jmd(if)+(ijf-1)*imd(if)+iif
804	write(unit(if)+1,rec=irec(if))
805	s ((field((l-1)imd(if)jmd(if)+(j-1)*imd(if)+i)
806	s ,i=iii,iif),j=iji,ijf)
807	enddo
808	if (writectl) then
809
810	file=fichier(if)
811	! WARNING! on reecrase le fichier .ctl a chaque ecriture
812	open(unit(if),file=trim(file)//'.ctl',
813	& form='formatted',status='unknown')
814	write(unit(if),'(a5,1x,a40)')
815	& 'DSET ','^'//trim(file)//'.dat'
816
817	write(unit(if),'(a12)') 'UNDEF 1.0E30'
818	write(unit(if),'(a5,1x,a40)') 'TITLE ',title(if)
819	call formcoord(unit(if),im,xd(iii,if),1.,.false.,'XDEF')
820	call formcoord(unit(if),jm,yd(iji,if),1.,.true.,'YDEF')
821	call formcoord(unit(if),lm,zd(1,if),1.,.false.,'ZDEF')
822	write(unit(if),'(a4,i10,a30)')
823	& 'TDEF ',itime(if),' LINEAR 07AUG1998 30MN '
824	write(unit(if),'(a4,2x,i5)') 'VARS',nvar(if)
825	do iv=1,nvar(if)
826	! print*,'if,var(iv,if),nld(iv,if),nld(iv,if)-1/nld(iv,if)'
827	! print*,if,var(iv,if),nld(iv,if),nld(iv,if)-1/nld(iv,if)
828	write(unit(if),1000) var(iv,if),nld(iv,if)-1/nld(iv,if)
829	& ,99,tvar(iv,if)
830	enddo
831	write(unit(if),'(a7)') 'ENDVARS'
832	!
833	1000 format(a5,3x,i4,i3,1x,a39)
834
835	close(unit(if))
836
837	endif ! writectl
838
839	return
840
841	END
842
843	subroutine inigrads(if,im
844	s ,x,fx,xmin,xmax,jm,y,ymin,ymax,fy,lm,z,fz
845	s ,dt,file,titlel)
846
847
848	implicit none
849
850	integer if,im,jm,lm,i,j,l
851	real x(im),y(jm),z(lm),fx,fy,fz,dt
852	real xmin,xmax,ymin,ymax
853	integer nf
854
855	character file10,titlel40
856
857	#include "gradsdef.h"
858
859	data unit/24,32,34,36,38,40,42,44,46,48/
860	data nf/0/
861
862	if (if.le.nf) stop'verifier les appels a inigrads'
863
864	print*,'Entree dans inigrads'
865
866	nf=if
867	title(if)=titlel
868	ivar(if)=0
869
870	fichier(if)=trim(file)
871
872	firsttime(if)=.true.
873	dtime(if)=dt
874
875	iid(if)=1
876	ifd(if)=im
877	imd(if)=im
878	do i=1,im
879	xd(i,if)=x(i)*fx
880	if(xd(i,if).lt.xmin) iid(if)=i+1
881	if(xd(i,if).le.xmax) ifd(if)=i
882	enddo
883	print*,'On stoke du point ',iid(if),' a ',ifd(if),' en x'
884
885	jid(if)=1
886	jfd(if)=jm
887	jmd(if)=jm
888	do j=1,jm
889	yd(j,if)=y(j)*fy
890	if(yd(j,if).gt.ymax) jid(if)=j+1
891	if(yd(j,if).ge.ymin) jfd(if)=j
892	enddo
893	print*,'On stoke du point ',jid(if),' a ',jfd(if),' en y'
894
895	print*,'Open de dat'
896	print*,'file=',file
897	print*,'fichier(if)=',fichier(if)
898
899	print,4(ifd(if)-iid(if))*(jfd(if)-jid(if))
900	print*,trim(file)//'.dat'
901
902	OPEN (unit(if)+1,FILE=trim(file)//'.dat',
903	s FORM='UNFORMATTED',
904	s ACCESS='DIRECT'
905	s ,RECL=4(ifd(if)-iid(if)+1)(jfd(if)-jid(if)+1))
906
907	print*,'Open de dat ok'
908
909	lmd(if)=lm
910	do l=1,lm
911	zd(l,if)=z(l)*fz
912	enddo
913
914	irec(if)=0
915	!CR
916	! print*,if,imd(if),jmd(if),lmd(if)
917	! print*,'if,imd(if),jmd(if),lmd(if)'
918
919	return
920	end
921	SUBROUTINE gr_dyn_fi(nfield,im,jm,ngrid,pdyn,pfi)
922	IMPLICIT NONE
923	!=======================================================================
924	! passage d'un champ de la grille scalaire a la grille physique
925	!=======================================================================
926
927	!-----------------------------------------------------------------------
928	! declarations:
929	! -------------
930
931	INTEGER im,jm,ngrid,nfield
932	REAL pdyn(im,jm,nfield)
933	REAL pfi(ngrid,nfield)
934
935	INTEGER j,ifield,ig
936
937	!-----------------------------------------------------------------------
938	! calcul:
939	! -------
940
941	IF(ngrid.NE.2+(jm-2)*(im-1).AND.ngrid.NE.1)
942	s STOP 'probleme de dim'
943	! traitement des poles
944	CALL SCOPY(nfield,pdyn,im*jm,pfi,ngrid)
945	CALL SCOPY(nfield,pdyn(1,jm,1),im*jm,pfi(ngrid,1),ngrid)
946
947	! traitement des point normaux
948	DO ifield=1,nfield
949	DO j=2,jm-1
950	ig=2+(j-2)*(im-1)
951	CALL SCOPY(im-1,pdyn(1,j,ifield),1,pfi(ig,ifield),1)
952	ENDDO
953	ENDDO
954
955	RETURN
956	END
957
958	SUBROUTINE disvert(pa,preff,ap,bp,dpres,presnivs,nivsigs,nivsig)
959
960	! Auteur : P. Le Van .
961	!
962	IMPLICIT NONE
963
964	#include "dimensions.h"
965	#include "paramet.h"
966	!
967	!=======================================================================
968	!
969	!
970	! s = sigma ** kappa : coordonnee verticale
971	! dsig(l) : epaisseur de la couche l ds la coord. s
972	! sig(l) : sigma a l'interface des couches l et l-1
973	! ds(l) : distance entre les couches l et l-1 en coord.s
974	!
975	!=======================================================================
976	!
977	REAL pa,preff
978	REAL ap(llmp1),bp(llmp1),dpres(llm),nivsigs(llm),nivsig(llmp1)
979	REAL presnivs(llm)
980	!
981	! declarations:
982	! -------------
983	!
984	REAL sig(llm+1),dsig(llm)
985	!
986	INTEGER l
987	REAL snorm
988	REAL alpha,beta,gama,delta,deltaz,h
989	INTEGER np,ierr
990	REAL pi,x
991
992	!-----------------------------------------------------------------------
993	!
994	pi=2.*ASIN(1.)
995
996	OPEN(99,file='sigma.def',status='old',form='formatted',
997	s iostat=ierr)
998
999	!-----------------------------------------------------------------------
1000	! cas 1 on lit les options dans sigma.def:
1001	! ----------------------------------------
1002
1003	IF (ierr.eq.0) THEN
1004
1005	print*,'WARNING!!! on lit les options dans sigma.def'
1006	READ(99,*) deltaz
1007	READ(99,*) h
1008	READ(99,*) beta
1009	READ(99,*) gama
1010	READ(99,*) delta
1011	READ(99,*) np
1012	CLOSE(99)
1013	alpha=deltaz/(llm*h)
1014	!
1015
1016	DO 1 l = 1, llm
1017	dsig(l) = (alpha+(1.-alpha)exp(-beta(llm-l)))*
1018	$ ( (tanh(gamal)/tanh(gamallm))**np +
1019	$ (1.-l/FLOAT(llm))*delta )
1020	1 CONTINUE
1021
1022	sig(1)=1.
1023	DO 101 l=1,llm-1
1024	sig(l+1)=sig(l)*(1.-dsig(l))/(1.+dsig(l))
1025	101 CONTINUE
1026	sig(llm+1)=0.
1027
1028	DO 2 l = 1, llm
1029	dsig(l) = sig(l)-sig(l+1)
1030	2 CONTINUE
1031	!
1032
1033	ELSE
1034	!-----------------------------------------------------------------------
1035	! cas 2 ancienne discretisation (LMD5...):
1036	! ----------------------------------------
1037
1038	PRINT*,'WARNING!!! Ancienne discretisation verticale'
1039
1040	h=7.
1041	snorm = 0.
1042	DO l = 1, llm
1043	x = 2.asin(1.) (FLOAT(l)-0.5) / float(llm+1)
1044	dsig(l) = 1.0 + 7.0 * SIN(x)**2
1045	snorm = snorm + dsig(l)
1046	ENDDO
1047	snorm = 1./snorm
1048	DO l = 1, llm
1049	dsig(l) = dsig(l)*snorm
1050	ENDDO
1051	sig(llm+1) = 0.
1052	DO l = llm, 1, -1
1053	sig(l) = sig(l+1) + dsig(l)
1054	ENDDO
1055
1056	ENDIF
1057
1058
1059	DO l=1,llm
1060	nivsigs(l) = FLOAT(l)
1061	ENDDO
1062
1063	DO l=1,llmp1
1064	nivsig(l)= FLOAT(l)
1065	ENDDO
1066
1067	!
1068	! .... Calculs de ap(l) et de bp(l) ....
1069	! .........................................
1070	!
1071	!
1072	! ..... pa et preff sont lus sur les fichiers start par lectba .....
1073	!
1074
1075	bp(llmp1) = 0.
1076
1077	DO l = 1, llm
1078	!c
1079	!cc ap(l) = 0.
1080	!cc bp(l) = sig(l)
1081
1082	bp(l) = EXP( 1. -1./( sig(l)*sig(l)) )
1083	ap(l) = pa * ( sig(l) - bp(l) )
1084	!
1085	ENDDO
1086	ap(llmp1) = pa * ( sig(llmp1) - bp(llmp1) )
1087
1088	PRINT *,' BP '
1089	PRINT *, bp
1090	PRINT *,' AP '
1091	PRINT *, ap
1092
1093	DO l = 1, llm
1094	dpres(l) = bp(l) - bp(l+1)
1095	presnivs(l) = 0.5 ( ap(l)+bp(l)preff + ap(l+1)+bp(l+1)*preff )
1096	ENDDO
1097
1098	PRINT *,' PRESNIVS '
1099	PRINT *,presnivs
1100
1101	RETURN
1102	END
1103
1104	!======================================================================
1105	SUBROUTINE read_tsurf1d(knon,knindex,sst_out)
1106
1107	! This subroutine specifies the surface temperature to be used in 1D simulations
1108
1109	USE dimphy, ONLY : klon
1110
1111	INTEGER, INTENT(IN) :: knon ! nomber of points on compressed grid
1112	INTEGER, DIMENSION(klon), INTENT(IN) :: knindex ! grid point number for compressed grid
1113	REAL, DIMENSION(klon), INTENT(OUT) :: sst_out ! tsurf used to force the single-column model
1114
1115	INTEGER :: i
1116	! COMMON defined in lmdz1d.F:
1117	real ts_cur
1118	common /sst_forcing/ts_cur
1119
1120	DO i = 1, knon
1121	sst_out(i) = ts_cur
1122	ENDDO
1123
1124	END SUBROUTINE read_tsurf1d
1125
1126	!===============================================================
1127	subroutine advect_vert(llm,w,dt,q,plev)
1128	!===============================================================
1129	! Schema amont pour l'advection verticale en 1D
1130	! w est la vitesse verticale dp/dt en Pa/s
1131	! Traitement en volumes finis
1132	! d / dt ( zm q ) = delta_z ( omega q )
1133	! d / dt ( zm ) = delta_z ( omega )
1134	! avec zm = delta_z ( p )
1135	! si * designe la valeur au pas de temps t+dt
1136	! zm(l) q(l) - zm(l) q(l) = w(l+1) q(l+1) - w(l) q(l)
1137	! zm*(l) -zm(l) = w(l+1) - w(l)
1138	! avec w=omega * dt
1139	!---------------------------------------------------------------
1140	implicit none
1141	! arguments
1142	integer llm
1143	real w(llm+1),q(llm),plev(llm+1),dt
1144
1145	! local
1146	integer l
1147	real zwq(llm+1),zm(llm+1),zw(llm+1)
1148	real qold
1149
1150	!---------------------------------------------------------------
1151
1152	do l=1,llm
1153	zw(l)=dt*w(l)
1154	zm(l)=plev(l)-plev(l+1)
1155	zwq(l)=q(l)*zw(l)
1156	enddo
1157	zwq(llm+1)=0.
1158	zw(llm+1)=0.
1159
1160	do l=1,llm
1161	qold=q(l)
1162	q(l)=(q(l)*zm(l)+zwq(l+1)-zwq(l))/(zm(l)+zw(l+1)-zw(l))
1163	print*,'ADV Q ',zm(l),zw(l),zwq(l),qold,q(l)
1164	enddo
1165
1166
1167	return
1168	end
1169
1170	!===============================================================
1171
1172
1173	SUBROUTINE advect_va(llm,omega,d_t_va,d_q_va,d_u_va,d_v_va,
1174	! q,temp,u,v,
1175	! play,plev)
1176	!itlmd
1177	!----------------------------------------------------------------------
1178	! Calcul de l'advection verticale (ascendance et subsidence) de
1179	! température et d'humidité. Hypothèse : ce qui rentre de l'extérieur
1180	! a les mêmes caractéristiques que l'air de la colonne 1D (WTG) ou
1181	! sans WTG rajouter une advection horizontale
1182	!----------------------------------------------------------------------
1183	implicit none
1184	#include "YOMCST.h"
1185	! argument
1186	integer llm
1187	real omega(llm+1),d_t_va(llm), d_q_va(llm,3)
1188	real d_u_va(llm), d_v_va(llm)
1189	real q(llm,3),temp(llm)
1190	real u(llm),v(llm)
1191	real play(llm),plev(llm+1)
1192	! interne
1193	integer l
1194	real alpha,omgdown,omgup
1195
1196	do l= 1,llm
1197	if(l.eq.1) then
1198	!si omgup pour la couche 1, alors tendance nulle
1199	omgdown=max(omega(2),0.0)
1200	alpha = rkappatemp(l)(1.+q(l,1)rv/rd)/(play(l)
1201	& (1.+q(l,1)))
1202	d_t_va(l)= alpha*(omgdown)-
1203	& omgdown*(temp(l)-temp(l+1))
1204	& /(play(l)-play(l+1))
1205
1206	d_q_va(l,:)= -omgdown*(q(l,:)-q(l+1,:))
1207	& /(play(l)-play(l+1))
1208
1209	d_u_va(l)= -omgdown*(u(l)-u(l+1))
1210	& /(play(l)-play(l+1))
1211	d_v_va(l)= -omgdown*(v(l)-v(l+1))
1212	& /(play(l)-play(l+1))
1213
1214
1215	elseif(l.eq.llm) then
1216	omgup=min(omega(l),0.0)
1217	alpha = rkappatemp(l)(1.+q(l,1)rv/rd)/(play(l)
1218	& (1.+q(l,1)))
1219	d_t_va(l)= alpha*(omgup)-
1220	!bug? & omgup*(temp(l-1)-temp(l))/(play(l-1)-plev(l))
1221	& omgup*(temp(l-1)-temp(l))/(play(l-1)-play(l))
1222	d_q_va(l,:)= -omgup*(q(l-1,:)-q(l,:))/(play(l-1)-play(l))
1223	d_u_va(l)= -omgup*(u(l-1)-u(l))/(play(l-1)-play(l))
1224	d_v_va(l)= -omgup*(v(l-1)-v(l))/(play(l-1)-play(l))
1225
1226	else
1227	omgup=min(omega(l),0.0)
1228	omgdown=max(omega(l+1),0.0)
1229	alpha = rkappatemp(l)(1.+q(l,1)rv/rd)/(play(l)
1230	& (1.+q(l,1)))
1231	d_t_va(l)= alpha*(omgup+omgdown)-
1232	& omgdown*(temp(l)-temp(l+1))
1233	& /(play(l)-play(l+1))-
1234	!bug? & omgup*(temp(l-1)-temp(l))/(play(l-1)-plev(l))
1235	& omgup*(temp(l-1)-temp(l))/(play(l-1)-play(l))
1236	! print*, ' ??? '
1237
1238	d_q_va(l,:)= -omgdown*(q(l,:)-q(l+1,:))
1239	& /(play(l)-play(l+1))-
1240	& omgup*(q(l-1,:)-q(l,:))/(play(l-1)-play(l))
1241	d_u_va(l)= -omgdown*(u(l)-u(l+1))
1242	& /(play(l)-play(l+1))-
1243	& omgup*(u(l-1)-u(l))/(play(l-1)-play(l))
1244	d_v_va(l)= -omgdown*(v(l)-v(l+1))
1245	& /(play(l)-play(l+1))-
1246	& omgup*(v(l-1)-v(l))/(play(l-1)-play(l))
1247
1248	endif
1249
1250	enddo
1251	!fin itlmd
1252	return
1253	end
1254	! SUBROUTINE lstendH(llm,omega,d_t_va,d_q_va,d_u_va,d_v_va,
1255	SUBROUTINE lstendH(llm,nqtot,omega,d_t_va,d_q_va,
1256	! q,temp,u,v,play)
1257	!itlmd
1258	!----------------------------------------------------------------------
1259	! Calcul de l'advection verticale (ascendance et subsidence) de
1260	! température et d'humidité. Hypothèse : ce qui rentre de l'extérieur
1261	! a les mêmes caractéristiques que l'air de la colonne 1D (WTG) ou
1262	! sans WTG rajouter une advection horizontale
1263	!----------------------------------------------------------------------
1264	implicit none
1265	#include "YOMCST.h"
1266	! argument
1267	integer llm,nqtot
1268	real omega(llm+1),d_t_va(llm), d_q_va(llm,nqtot)
1269	! real d_u_va(llm), d_v_va(llm)
1270	real q(llm,nqtot),temp(llm)
1271	real u(llm),v(llm)
1272	real play(llm)
1273	real cor(llm)
1274	! real dph(llm),dudp(llm),dvdp(llm),dqdp(llm),dtdp(llm)
1275	real dph(llm),dqdp(llm),dtdp(llm)
1276	! interne
1277	integer l,k
1278	real alpha,omdn,omup
1279
1280	! dudp=0.
1281	! dvdp=0.
1282	dqdp=0.
1283	dtdp=0.
1284	! d_u_va=0.
1285	! d_v_va=0.
1286
1287	cor(:) = rkappatemp(1.+q(:,1)rv/rd)/(play(1.+q(:,1)))
1288
1289
1290	do k=2,llm-1
1291
1292	dph (k-1) = (play(k )- play(k-1 ))
1293	! dudp (k-1) = (u (k )- u (k-1 ))/dph(k-1)
1294	! dvdp (k-1) = (v (k )- v (k-1 ))/dph(k-1)
1295	dqdp (k-1) = (q (k,1)- q (k-1,1))/dph(k-1)
1296	dtdp (k-1) = (temp(k )- temp(k-1 ))/dph(k-1)
1297
1298	enddo
1299
1300	! dudp ( llm ) = dudp ( llm-1 )
1301	! dvdp ( llm ) = dvdp ( llm-1 )
1302	dqdp ( llm ) = dqdp ( llm-1 )
1303	dtdp ( llm ) = dtdp ( llm-1 )
1304
1305	do k=2,llm-1
1306	omdn=max(0.0,omega(k+1))
1307	omup=min(0.0,omega( k ))
1308
1309	! d_u_va(k) = -omdndudp(k)-omupdudp(k-1)
1310	! d_v_va(k) = -omdndvdp(k)-omupdvdp(k-1)
1311	d_q_va(k,1)= -omdndqdp(k)-omupdqdp(k-1)
1312	d_t_va(k) = -omdndtdp(k)-omupdtdp(k-1)
1313	: +(omup+omdn)*cor(k)
1314	enddo
1315
1316	omdn=max(0.0,omega( 2 ))
1317	omup=min(0.0,omega(llm))
1318	! d_u_va( 1 ) = -omdn*dudp( 1 )
1319	! d_u_va(llm) = -omup*dudp(llm)
1320	! d_v_va( 1 ) = -omdn*dvdp( 1 )
1321	! d_v_va(llm) = -omup*dvdp(llm)
1322	d_q_va( 1 ,1) = -omdn*dqdp( 1 )
1323	d_q_va(llm,1) = -omup*dqdp(llm)
1324	d_t_va( 1 ) = -omdndtdp( 1 )+omdncor( 1 )
1325	d_t_va(llm) = -omupdtdp(llm)!+omupcor(llm)
1326
1327	! if(abs(rlat(1))>10.) then
1328	! Calculate the tendency due agestrophic motions
1329	! du_age = fcoriolis*(v-vg)
1330	! dv_age = fcoriolis*(ug-u)
1331	! endif
1332
1333	! call writefield_phy('d_t_va',d_t_va,llm)
1334
1335	return
1336	end
1337
1338	!======================================================================
1339	SUBROUTINE read_togacoare(fich_toga,nlev_toga,nt_toga
1340	: ,ts_toga,plev_toga,t_toga,q_toga,u_toga,v_toga,w_toga
1341	: ,ht_toga,vt_toga,hq_toga,vq_toga)
1342	implicit none
1343
1344	c-------------------------------------------------------------------------
1345	c Read TOGA-COARE forcing data
1346	c-------------------------------------------------------------------------
1347
1348	integer nlev_toga,nt_toga
1349	real ts_toga(nt_toga),plev_toga(nlev_toga,nt_toga)
1350	real t_toga(nlev_toga,nt_toga),q_toga(nlev_toga,nt_toga)
1351	real u_toga(nlev_toga,nt_toga),v_toga(nlev_toga,nt_toga)
1352	real w_toga(nlev_toga,nt_toga)
1353	real ht_toga(nlev_toga,nt_toga),vt_toga(nlev_toga,nt_toga)
1354	real hq_toga(nlev_toga,nt_toga),vq_toga(nlev_toga,nt_toga)
1355	character*80 fich_toga
1356
1357	integer no,l,k,ip
1358	real riy,rim,rid,rih,bid
1359
1360	integer iy,im,id,ih
1361
1362	real plev_min
1363
1364	plev_min = 55. ! pas de tendance de vap. d eau au-dessus de 55 hPa
1365
1366	open(21,file=trim(fich_toga),form='formatted')
1367	read(21,'(a)')
1368	do ip = 1, nt_toga
1369	read(21,'(a)')
1370	read(21,'(a)')
1371	read(21,223) iy, im, id, ih, bid, ts_toga(ip), bid,bid,bid,bid
1372	read(21,'(a)')
1373	read(21,'(a)')
1374
1375	do k = 1, nlev_toga
1376	read(21,230) plev_toga(k,ip), t_toga(k,ip), q_toga(k,ip)
1377	: ,u_toga(k,ip), v_toga(k,ip), w_toga(k,ip)
1378	: ,ht_toga(k,ip), vt_toga(k,ip), hq_toga(k,ip), vq_toga(k,ip)
1379
1380	! conversion in SI units:
1381	t_toga(k,ip)=t_toga(k,ip)+273.15 ! K
1382	q_toga(k,ip)=q_toga(k,ip)*0.001 ! kg/kg
1383	w_toga(k,ip)=w_toga(k,ip)*100./3600. ! Pa/s
1384	! no water vapour tendency above 55 hPa
1385	if (plev_toga(k,ip) .lt. plev_min) then
1386	q_toga(k,ip) = 0.
1387	hq_toga(k,ip) = 0.
1388	vq_toga(k,ip) =0.
1389	endif
1390	enddo
1391
1392	ts_toga(ip)=ts_toga(ip)+273.15 ! K
1393	enddo
1394	close(21)
1395
1396	223 format(4i3,6f8.2)
1397	226 format(f7.1,1x,10f8.2)
1398	227 format(f7.1,1x,1p,4e11.3)
1399	230 format(6f9.3,4e11.3)
1400
1401	return
1402	end
1403
1404	!=====================================================================
1405	subroutine read_twpice(fich_twpice,nlevel,ntime
1406	: ,T_srf,plev,T,q,u,v,omega
1407	: ,T_adv_h,T_adv_v,q_adv_h,q_adv_v)
1408
1409	!program reading forcings of the TWP-ICE experiment
1410
1411	! use netcdf
1412
1413	implicit none
1414
1415	#include "netcdf.inc"
1416
1417	integer ntime,nlevel
1418	integer l,k
1419	character*80 :: fich_twpice
1420	real*8 time(ntime)
1421	real*8 lat, lon, alt, phis
1422	real*8 lev(nlevel)
1423	real*8 plev(nlevel,ntime)
1424
1425	real*8 T(nlevel,ntime)
1426	real*8 q(nlevel,ntime),u(nlevel,ntime)
1427	real*8 v(nlevel,ntime)
1428	real*8 omega(nlevel,ntime), div(nlevel,ntime)
1429	real*8 T_adv_h(nlevel,ntime)
1430	real*8 T_adv_v(nlevel,ntime), q_adv_h(nlevel,ntime)
1431	real*8 q_adv_v(nlevel,ntime)
1432	real*8 s(nlevel,ntime), s_adv_h(nlevel,ntime)
1433	real*8 s_adv_v(nlevel,ntime)
1434	real*8 p_srf_aver(ntime), p_srf_center(ntime)
1435	real*8 T_srf(ntime)
1436
1437	integer nid, ierr
1438	integer nbvar3d
1439	parameter(nbvar3d=20)
1440	integer var3didin(nbvar3d)
1441
1442	ierr = NF_OPEN(fich_twpice,NF_NOWRITE,nid)
1443	if (ierr.NE.NF_NOERR) then
1444	write(,) 'ERROR: Pb opening forcings cdf file '
1445	write(,) NF_STRERROR(ierr)
1446	stop ""
1447	endif
1448
1449	ierr=NF_INQ_VARID(nid,"lat",var3didin(1))
1450	if(ierr/=NF_NOERR) then
1451	write(,) NF_STRERROR(ierr)
1452	stop 'lat'
1453	endif
1454
1455	ierr=NF_INQ_VARID(nid,"lon",var3didin(2))
1456	if(ierr/=NF_NOERR) then
1457	write(,) NF_STRERROR(ierr)
1458	stop 'lon'
1459	endif
1460
1461	ierr=NF_INQ_VARID(nid,"alt",var3didin(3))
1462	if(ierr/=NF_NOERR) then
1463	write(,) NF_STRERROR(ierr)
1464	stop 'alt'
1465	endif
1466
1467	ierr=NF_INQ_VARID(nid,"phis",var3didin(4))
1468	if(ierr/=NF_NOERR) then
1469	write(,) NF_STRERROR(ierr)
1470	stop 'phis'
1471	endif
1472
1473	ierr=NF_INQ_VARID(nid,"T",var3didin(5))
1474	if(ierr/=NF_NOERR) then
1475	write(,) NF_STRERROR(ierr)
1476	stop 'T'
1477	endif
1478
1479	ierr=NF_INQ_VARID(nid,"q",var3didin(6))
1480	if(ierr/=NF_NOERR) then
1481	write(,) NF_STRERROR(ierr)
1482	stop 'q'
1483	endif
1484
1485	ierr=NF_INQ_VARID(nid,"u",var3didin(7))
1486	if(ierr/=NF_NOERR) then
1487	write(,) NF_STRERROR(ierr)
1488	stop 'u'
1489	endif
1490
1491	ierr=NF_INQ_VARID(nid,"v",var3didin(8))
1492	if(ierr/=NF_NOERR) then
1493	write(,) NF_STRERROR(ierr)
1494	stop 'v'
1495	endif
1496
1497	ierr=NF_INQ_VARID(nid,"omega",var3didin(9))
1498	if(ierr/=NF_NOERR) then
1499	write(,) NF_STRERROR(ierr)
1500	stop 'omega'
1501	endif
1502
1503	ierr=NF_INQ_VARID(nid,"div",var3didin(10))
1504	if(ierr/=NF_NOERR) then
1505	write(,) NF_STRERROR(ierr)
1506	stop 'div'
1507	endif
1508
1509	ierr=NF_INQ_VARID(nid,"T_adv_h",var3didin(11))
1510	if(ierr/=NF_NOERR) then
1511	write(,) NF_STRERROR(ierr)
1512	stop 'T_adv_h'
1513	endif
1514
1515	ierr=NF_INQ_VARID(nid,"T_adv_v",var3didin(12))
1516	if(ierr/=NF_NOERR) then
1517	write(,) NF_STRERROR(ierr)
1518	stop 'T_adv_v'
1519	endif
1520
1521	ierr=NF_INQ_VARID(nid,"q_adv_h",var3didin(13))
1522	if(ierr/=NF_NOERR) then
1523	write(,) NF_STRERROR(ierr)
1524	stop 'q_adv_h'
1525	endif
1526
1527	ierr=NF_INQ_VARID(nid,"q_adv_v",var3didin(14))
1528	if(ierr/=NF_NOERR) then
1529	write(,) NF_STRERROR(ierr)
1530	stop 'q_adv_v'
1531	endif
1532
1533	ierr=NF_INQ_VARID(nid,"s",var3didin(15))
1534	if(ierr/=NF_NOERR) then
1535	write(,) NF_STRERROR(ierr)
1536	stop 's'
1537	endif
1538
1539	ierr=NF_INQ_VARID(nid,"s_adv_h",var3didin(16))
1540	if(ierr/=NF_NOERR) then
1541	write(,) NF_STRERROR(ierr)
1542	stop 's_adv_h'
1543	endif
1544
1545	ierr=NF_INQ_VARID(nid,"s_adv_v",var3didin(17))
1546	if(ierr/=NF_NOERR) then
1547	write(,) NF_STRERROR(ierr)
1548	stop 's_adv_v'
1549	endif
1550
1551	ierr=NF_INQ_VARID(nid,"p_srf_aver",var3didin(18))
1552	if(ierr/=NF_NOERR) then
1553	write(,) NF_STRERROR(ierr)
1554	stop 'p_srf_aver'
1555	endif
1556
1557	ierr=NF_INQ_VARID(nid,"p_srf_center",var3didin(19))
1558	if(ierr/=NF_NOERR) then
1559	write(,) NF_STRERROR(ierr)
1560	stop 'p_srf_center'
1561	endif
1562
1563	ierr=NF_INQ_VARID(nid,"T_srf",var3didin(20))
1564	if(ierr/=NF_NOERR) then
1565	write(,) NF_STRERROR(ierr)
1566	stop 'T_srf'
1567	endif
1568
1569	!dimensions lecture
1570	call catchaxis(nid,ntime,nlevel,time,lev,ierr)
1571
1572	!pressure
1573	do l=1,ntime
1574	do k=1,nlevel
1575	plev(k,l)=lev(k)
1576	enddo
1577	enddo
1578
1579	#ifdef NC_DOUBLE
1580	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(1),lat)
1581	#else
1582	ierr = NF_GET_VAR_REAL(nid,var3didin(1),lat)
1583	#endif
1584	if(ierr/=NF_NOERR) then
1585	write(,) NF_STRERROR(ierr)
1586	stop "getvarup"
1587	endif
1588	! write(,)'lecture lat ok',lat
1589
1590	#ifdef NC_DOUBLE
1591	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(2),lon)
1592	#else
1593	ierr = NF_GET_VAR_REAL(nid,var3didin(2),lon)
1594	#endif
1595	if(ierr/=NF_NOERR) then
1596	write(,) NF_STRERROR(ierr)
1597	stop "getvarup"
1598	endif
1599	! write(,)'lecture lon ok',lon
1600
1601	#ifdef NC_DOUBLE
1602	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(3),alt)
1603	#else
1604	ierr = NF_GET_VAR_REAL(nid,var3didin(3),alt)
1605	#endif
1606	if(ierr/=NF_NOERR) then
1607	write(,) NF_STRERROR(ierr)
1608	stop "getvarup"
1609	endif
1610	! write(,)'lecture alt ok',alt
1611
1612	#ifdef NC_DOUBLE
1613	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(4),phis)
1614	#else
1615	ierr = NF_GET_VAR_REAL(nid,var3didin(4),phis)
1616	#endif
1617	if(ierr/=NF_NOERR) then
1618	write(,) NF_STRERROR(ierr)
1619	stop "getvarup"
1620	endif
1621	! write(,)'lecture phis ok',phis
1622
1623	#ifdef NC_DOUBLE
1624	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(5),T)
1625	#else
1626	ierr = NF_GET_VAR_REAL(nid,var3didin(5),T)
1627	#endif
1628	if(ierr/=NF_NOERR) then
1629	write(,) NF_STRERROR(ierr)
1630	stop "getvarup"
1631	endif
1632	! write(,)'lecture T ok'
1633
1634	#ifdef NC_DOUBLE
1635	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(6),q)
1636	#else
1637	ierr = NF_GET_VAR_REAL(nid,var3didin(6),q)
1638	#endif
1639	if(ierr/=NF_NOERR) then
1640	write(,) NF_STRERROR(ierr)
1641	stop "getvarup"
1642	endif
1643	! write(,)'lecture q ok'
1644	!q in kg/kg
1645	do l=1,ntime
1646	do k=1,nlevel
1647	q(k,l)=q(k,l)/1000.
1648	enddo
1649	enddo
1650	#ifdef NC_DOUBLE
1651	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(7),u)
1652	#else
1653	ierr = NF_GET_VAR_REAL(nid,var3didin(7),u)
1654	#endif
1655	if(ierr/=NF_NOERR) then
1656	write(,) NF_STRERROR(ierr)
1657	stop "getvarup"
1658	endif
1659	! write(,)'lecture u ok'
1660
1661	#ifdef NC_DOUBLE
1662	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(8),v)
1663	#else
1664	ierr = NF_GET_VAR_REAL(nid,var3didin(8),v)
1665	#endif
1666	if(ierr/=NF_NOERR) then
1667	write(,) NF_STRERROR(ierr)
1668	stop "getvarup"
1669	endif
1670	! write(,)'lecture v ok'
1671
1672	#ifdef NC_DOUBLE
1673	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(9),omega)
1674	#else
1675	ierr = NF_GET_VAR_REAL(nid,var3didin(9),omega)
1676	#endif
1677	if(ierr/=NF_NOERR) then
1678	write(,) NF_STRERROR(ierr)
1679	stop "getvarup"
1680	endif
1681	! write(,)'lecture omega ok'
1682	!omega in mb/hour
1683	do l=1,ntime
1684	do k=1,nlevel
1685	omega(k,l)=omega(k,l)*100./3600.
1686	enddo
1687	enddo
1688
1689	#ifdef NC_DOUBLE
1690	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(10),div)
1691	#else
1692	ierr = NF_GET_VAR_REAL(nid,var3didin(10),div)
1693	#endif
1694	if(ierr/=NF_NOERR) then
1695	write(,) NF_STRERROR(ierr)
1696	stop "getvarup"
1697	endif
1698	! write(,)'lecture div ok'
1699
1700	#ifdef NC_DOUBLE
1701	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(11),T_adv_h)
1702	#else
1703	ierr = NF_GET_VAR_REAL(nid,var3didin(11),T_adv_h)
1704	#endif
1705	if(ierr/=NF_NOERR) then
1706	write(,) NF_STRERROR(ierr)
1707	stop "getvarup"
1708	endif
1709	! write(,)'lecture T_adv_h ok'
1710	!T adv in K/s
1711	do l=1,ntime
1712	do k=1,nlevel
1713	T_adv_h(k,l)=T_adv_h(k,l)/3600.
1714	enddo
1715	enddo
1716
1717
1718	#ifdef NC_DOUBLE
1719	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(12),T_adv_v)
1720	#else
1721	ierr = NF_GET_VAR_REAL(nid,var3didin(12),T_adv_v)
1722	#endif
1723	if(ierr/=NF_NOERR) then
1724	write(,) NF_STRERROR(ierr)
1725	stop "getvarup"
1726	endif
1727	! write(,)'lecture T_adv_v ok'
1728	!T adv in K/s
1729	do l=1,ntime
1730	do k=1,nlevel
1731	T_adv_v(k,l)=T_adv_v(k,l)/3600.
1732	enddo
1733	enddo
1734
1735	#ifdef NC_DOUBLE
1736	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(13),q_adv_h)
1737	#else
1738	ierr = NF_GET_VAR_REAL(nid,var3didin(13),q_adv_h)
1739	#endif
1740	if(ierr/=NF_NOERR) then
1741	write(,) NF_STRERROR(ierr)
1742	stop "getvarup"
1743	endif
1744	! write(,)'lecture q_adv_h ok'
1745	!q adv in kg/kg/s
1746	do l=1,ntime
1747	do k=1,nlevel
1748	q_adv_h(k,l)=q_adv_h(k,l)/1000./3600.
1749	enddo
1750	enddo
1751
1752
1753	#ifdef NC_DOUBLE
1754	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(14),q_adv_v)
1755	#else
1756	ierr = NF_GET_VAR_REAL(nid,var3didin(14),q_adv_v)
1757	#endif
1758	if(ierr/=NF_NOERR) then
1759	write(,) NF_STRERROR(ierr)
1760	stop "getvarup"
1761	endif
1762	! write(,)'lecture q_adv_v ok'
1763	!q adv in kg/kg/s
1764	do l=1,ntime
1765	do k=1,nlevel
1766	q_adv_v(k,l)=q_adv_v(k,l)/1000./3600.
1767	enddo
1768	enddo
1769
1770
1771	#ifdef NC_DOUBLE
1772	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(15),s)
1773	#else
1774	ierr = NF_GET_VAR_REAL(nid,var3didin(15),s)
1775	#endif
1776	if(ierr/=NF_NOERR) then
1777	write(,) NF_STRERROR(ierr)
1778	stop "getvarup"
1779	endif
1780
1781	#ifdef NC_DOUBLE
1782	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(16),s_adv_h)
1783	#else
1784	ierr = NF_GET_VAR_REAL(nid,var3didin(16),s_adv_h)
1785	#endif
1786	if(ierr/=NF_NOERR) then
1787	write(,) NF_STRERROR(ierr)
1788	stop "getvarup"
1789	endif
1790
1791	#ifdef NC_DOUBLE
1792	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(17),s_adv_v)
1793	#else
1794	ierr = NF_GET_VAR_REAL(nid,var3didin(17),s_adv_v)
1795	#endif
1796	if(ierr/=NF_NOERR) then
1797	write(,) NF_STRERROR(ierr)
1798	stop "getvarup"
1799	endif
1800
1801	#ifdef NC_DOUBLE
1802	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(18),p_srf_aver)
1803	#else
1804	ierr = NF_GET_VAR_REAL(nid,var3didin(18),p_srf_aver)
1805	#endif
1806	if(ierr/=NF_NOERR) then
1807	write(,) NF_STRERROR(ierr)
1808	stop "getvarup"
1809	endif
1810
1811	#ifdef NC_DOUBLE
1812	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(19),p_srf_center)
1813	#else
1814	ierr = NF_GET_VAR_REAL(nid,var3didin(19),p_srf_center)
1815	#endif
1816	if(ierr/=NF_NOERR) then
1817	write(,) NF_STRERROR(ierr)
1818	stop "getvarup"
1819	endif
1820
1821	#ifdef NC_DOUBLE
1822	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(20),T_srf)
1823	#else
1824	ierr = NF_GET_VAR_REAL(nid,var3didin(20),T_srf)
1825	#endif
1826	if(ierr/=NF_NOERR) then
1827	write(,) NF_STRERROR(ierr)
1828	stop "getvarup"
1829	endif
1830	! write(,)'lecture T_srf ok', T_srf
1831
1832	return
1833	end subroutine read_twpice
1834	!=====================================================================
1835	subroutine catchaxis(nid,ttm,llm,time,lev,ierr)
1836
1837	! use netcdf
1838
1839	implicit none
1840	#include "netcdf.inc"
1841	integer nid,ttm,llm
1842	real*8 time(ttm)
1843	real*8 lev(llm)
1844	integer ierr
1845
1846	integer i
1847	integer timevar,levvar
1848	integer timelen,levlen
1849	integer timedimin,levdimin
1850
1851	! Control & lecture on dimensions
1852	! ===============================
1853	ierr=NF_INQ_DIMID(nid,"time",timedimin)
1854	ierr=NF_INQ_VARID(nid,"time",timevar)
1855	if (ierr.NE.NF_NOERR) then
1856	write(,) 'ERROR: Field <time> is missing'
1857	stop ""
1858	endif
1859	ierr=NF_INQ_DIMLEN(nid,timedimin,timelen)
1860
1861	ierr=NF_INQ_DIMID(nid,"lev",levdimin)
1862	ierr=NF_INQ_VARID(nid,"lev",levvar)
1863	if (ierr.NE.NF_NOERR) then
1864	write(,) 'ERROR: Field <lev> is lacking'
1865	stop ""
1866	endif
1867	ierr=NF_INQ_DIMLEN(nid,levdimin,levlen)
1868
1869	if((timelen/=ttm).or.(levlen/=llm)) then
1870	write(,) 'ERROR: Not the good lenght for axis'
1871	write(,) 'longitude: ',timelen,ttm+1
1872	write(,) 'latitude: ',levlen,llm
1873	stop ""
1874	endif
1875
1876	!#ifdef NC_DOUBLE
1877	ierr = NF_GET_VAR_DOUBLE(nid,timevar,time)
1878	ierr = NF_GET_VAR_DOUBLE(nid,levvar,lev)
1879	!#else
1880	! ierr = NF_GET_VAR_REAL(nid,timevar,time)
1881	! ierr = NF_GET_VAR_REAL(nid,levvar,lev)
1882	!#endif
1883
1884	return
1885	end
1886
1887	!======================================================================
1888	SUBROUTINE read_rico(fich_rico,nlev_rico,ps_rico,play
1889	: ,ts_rico,t_rico,q_rico,u_rico,v_rico,w_rico
1890	: ,dth_dyn,dqh_dyn)
1891	implicit none
1892
1893	c-------------------------------------------------------------------------
1894	c Read RICO forcing data
1895	c-------------------------------------------------------------------------
1896	#include "dimensions.h"
1897
1898
1899	integer nlev_rico
1900	real ts_rico,ps_rico
1901	real t_rico(llm),q_rico(llm)
1902	real u_rico(llm),v_rico(llm)
1903	real w_rico(llm)
1904	real dth_dyn(llm)
1905	real dqh_dyn(llm)
1906
1907
1908	real play(llm),zlay(llm)
1909
1910
1911	real prico(nlev_rico),zrico(nlev_rico)
1912
1913	character*80 fich_rico
1914
1915	integer k,l
1916
1917
1918	print*,fich_rico
1919	open(21,file=trim(fich_rico),form='formatted')
1920	do k=1,llm
1921	zlay(k)=0.
1922	enddo
1923
1924	read(21,*) ps_rico,ts_rico
1925	prico(1)=ps_rico
1926	zrico(1)=0.0
1927	do l=2,nlev_rico
1928	read(21,*) k,prico(l),zrico(l)
1929	enddo
1930	close(21)
1931
1932	do k=1,llm
1933	do l=1,80
1934	if(prico(l)>play(k)) then
1935	if(play(k)>prico(l+1)) then
1936	zlay(k)=zrico(l)+(play(k)-prico(l)) *
1937	& (zrico(l+1)-zrico(l))/(prico(l+1)-prico(l))
1938	else
1939	zlay(k)=zrico(l)+(play(k)-prico(80))*
1940	& (zrico(81)-zrico(80))/(prico(81)-prico(80))
1941	endif
1942	endif
1943	enddo
1944	print*,k,zlay(k)
1945	! U
1946	if(0 < zlay(k) .and. zlay(k) < 4000) then
1947	u_rico(k)=-9.9 + (-1.9 + 9.9)*zlay(k)/4000
1948	elseif(4000 < zlay(k) .and. zlay(k) < 12000) then
1949	u_rico(k)= -1.9 + (30.0 + 1.9) /
1950	: (12000 - 4000) * (zlay(k) - 4000)
1951	elseif(12000 < zlay(k) .and. zlay(k) < 13000) then
1952	u_rico(k)=30.0
1953	elseif(13000 < zlay(k) .and. zlay(k) < 20000) then
1954	u_rico(k)=30.0 - (30.0) /
1955	: (20000 - 13000) * (zlay(k) - 13000)
1956	else
1957	u_rico(k)=0.0
1958	endif
1959
1960	!Q_v
1961	if(0 < zlay(k) .and. zlay(k) < 740) then
1962	q_rico(k)=16.0 + (13.8 - 16.0) / (740) * zlay(k)
1963	elseif(740 < zlay(k) .and. zlay(k) < 3260) then
1964	q_rico(k)=13.8 + (2.4 - 13.8) /
1965	: (3260 - 740) * (zlay(k) - 740)
1966	elseif(3260 < zlay(k) .and. zlay(k) < 4000) then
1967	q_rico(k)=2.4 + (1.8 - 2.4) /
1968	: (4000 - 3260) * (zlay(k) - 3260)
1969	elseif(4000 < zlay(k) .and. zlay(k) < 9000) then
1970	q_rico(k)=1.8 + (0 - 1.8) /
1971	: (10000 - 4000) * (zlay(k) - 4000)
1972	else
1973	q_rico(k)=0.0
1974	endif
1975
1976	!T
1977	if(0 < zlay(k) .and. zlay(k) < 740) then
1978	t_rico(k)=299.2 + (292.0 - 299.2) / (740) * zlay(k)
1979	elseif(740 < zlay(k) .and. zlay(k) < 4000) then
1980	t_rico(k)=292.0 + (278.0 - 292.0) /
1981	: (4000 - 740) * (zlay(k) - 740)
1982	elseif(4000 < zlay(k) .and. zlay(k) < 15000) then
1983	t_rico(k)=278.0 + (203.0 - 278.0) /
1984	: (15000 - 4000) * (zlay(k) - 4000)
1985	elseif(15000 < zlay(k) .and. zlay(k) < 17500) then
1986	t_rico(k)=203.0 + (194.0 - 203.0) /
1987	: (17500 - 15000)* (zlay(k) - 15000)
1988	elseif(17500 < zlay(k) .and. zlay(k) < 20000) then
1989	t_rico(k)=194.0 + (206.0 - 194.0) /
1990	: (20000 - 17500)* (zlay(k) - 17500)
1991	elseif(20000 < zlay(k) .and. zlay(k) < 60000) then
1992	t_rico(k)=206.0 + (270.0 - 206.0) /
1993	: (60000 - 20000)* (zlay(k) - 20000)
1994	endif
1995
1996	! W
1997	if(0 < zlay(k) .and. zlay(k) < 2260 ) then
1998	w_rico(k)=- (0.005/2260) * zlay(k)
1999	elseif(2260 < zlay(k) .and. zlay(k) < 4000 ) then
2000	w_rico(k)=- 0.005
2001	elseif(4000 < zlay(k) .and. zlay(k) < 5000 ) then
2002	w_rico(k)=- 0.005 + (0.005/ (5000 - 4000)) * (zlay(k) - 4000)
2003	else
2004	w_rico(k)=0.0
2005	endif
2006
2007	! dThrz+dTsw0+dTlw0
2008	if(0 < zlay(k) .and. zlay(k) < 4000) then
2009	dth_dyn(k)=- 2.51 / 86400 + (-2.18 + 2.51 )/
2010	: (864004000) zlay(k)
2011	elseif(4000 < zlay(k) .and. zlay(k) < 5000) then
2012	dth_dyn(k)=- 2.18 / 86400 + ( 2.18 ) /
2013	: (86400(5000 - 4000)) (zlay(k) - 4000)
2014	else
2015	dth_dyn(k)=0.0
2016	endif
2017	! dQhrz
2018	if(0 < zlay(k) .and. zlay(k) < 3000) then
2019	dqh_dyn(k)=-1.0 / 86400 + (0.345 + 1.0)/
2020	: (864003000) (zlay(k))
2021	elseif(3000 < zlay(k) .and. zlay(k) < 4000) then
2022	dqh_dyn(k)=0.345 / 86400
2023	elseif(4000 < zlay(k) .and. zlay(k) < 5000) then
2024	dqh_dyn(k)=0.345 / 86400 +
2025	+ (-0.345)/(86400 * (5000 - 4000)) * (zlay(k)-4000)
2026	else
2027	dqh_dyn(k)=0.0
2028	endif
2029
2030	!? if(play(k)>6e4) then
2031	!? ratqs0(1,k)=ratqsbas*(plev(1)-play(k))/(plev(1)-6e4)
2032	!? elseif((play(k)>3e4).and.(play(k)<6e4)) then
2033	!? ratqs0(1,k)=ratqsbas+(ratqshaut-ratqsbas)&
2034	!? *(6e4-play(k))/(6e4-3e4)
2035	!? else
2036	!? ratqs0(1,k)=ratqshaut
2037	!? endif
2038
2039	enddo
2040
2041	do k=1,llm
2042	q_rico(k)=q_rico(k)/1e3
2043	dqh_dyn(k)=dqh_dyn(k)/1e3
2044	v_rico(k)=-3.8
2045	enddo
2046
2047	return
2048	end
2049
2050	!=====================================================================
2051	c-------------------------------------------------------------------------
2052	SUBROUTINE read_armcu(fich_armcu,nlev_armcu,nt_armcu,
2053	: sens,flat,adv_theta,rad_theta,adv_qt)
2054	implicit none
2055
2056	c-------------------------------------------------------------------------
2057	c Read ARM_CU case forcing data
2058	c-------------------------------------------------------------------------
2059
2060	integer nlev_armcu,nt_armcu
2061	real sens(nt_armcu),flat(nt_armcu)
2062	real adv_theta(nt_armcu),rad_theta(nt_armcu),adv_qt(nt_armcu)
2063	character*80 fich_armcu
2064
2065	integer no,l,k,ip
2066	real riy,rim,rid,rih,bid
2067
2068	integer iy,im,id,ih,in
2069
2070	open(21,file=trim(fich_armcu),form='formatted')
2071	read(21,'(a)')
2072	do ip = 1, nt_armcu
2073	read(21,'(a)')
2074	read(21,'(a)')
2075	read(21,223) iy, im, id, ih, in, sens(ip),flat(ip),
2076	: adv_theta(ip),rad_theta(ip),adv_qt(ip)
2077	print *,'forcages=',iy,im,id,ih,in, sens(ip),flat(ip),
2078	: adv_theta(ip),rad_theta(ip),adv_qt(ip)
2079	enddo
2080	close(21)
2081
2082	223 format(5i3,5f8.3)
2083	226 format(f7.1,1x,10f8.2)
2084	227 format(f7.1,1x,1p,4e11.3)
2085	230 format(6f9.3,4e11.3)
2086
2087	return
2088	end
2089
2090	!=====================================================================
2091	SUBROUTINE interp_toga_vertical(play,nlev_toga,plev_prof
2092	: ,t_prof,q_prof,u_prof,v_prof,w_prof
2093	: ,ht_prof,vt_prof,hq_prof,vq_prof
2094	: ,t_mod,q_mod,u_mod,v_mod,w_mod
2095	: ,ht_mod,vt_mod,hq_mod,vq_mod,mxcalc)
2096
2097	implicit none
2098
2099	#include "dimensions.h"
2100
2101	c-------------------------------------------------------------------------
2102	c Vertical interpolation of TOGA-COARE forcing data onto model levels
2103	c-------------------------------------------------------------------------
2104
2105	integer nlevmax
2106	parameter (nlevmax=41)
2107	integer nlev_toga,mxcalc
2108	! real play(llm), plev_prof(nlevmax)
2109	! real t_prof(nlevmax),q_prof(nlevmax)
2110	! real u_prof(nlevmax),v_prof(nlevmax), w_prof(nlevmax)
2111	! real ht_prof(nlevmax),vt_prof(nlevmax)
2112	! real hq_prof(nlevmax),vq_prof(nlevmax)
2113
2114	real play(llm), plev_prof(nlev_toga)
2115	real t_prof(nlev_toga),q_prof(nlev_toga)
2116	real u_prof(nlev_toga),v_prof(nlev_toga), w_prof(nlev_toga)
2117	real ht_prof(nlev_toga),vt_prof(nlev_toga)
2118	real hq_prof(nlev_toga),vq_prof(nlev_toga)
2119
2120	real t_mod(llm),q_mod(llm)
2121	real u_mod(llm),v_mod(llm), w_mod(llm)
2122	real ht_mod(llm),vt_mod(llm)
2123	real hq_mod(llm),vq_mod(llm)
2124
2125	integer l,k,k1,k2,kp
2126	real aa,frac,frac1,frac2,fact
2127
2128	do l = 1, llm
2129
2130	if (play(l).ge.plev_prof(nlev_toga)) then
2131
2132	mxcalc=l
2133	k1=0
2134	k2=0
2135
2136	if (play(l).le.plev_prof(1)) then
2137
2138	do k = 1, nlev_toga-1
2139	if (play(l).le.plev_prof(k)
2140	: .and. play(l).gt.plev_prof(k+1)) then
2141	k1=k
2142	k2=k+1
2143	endif
2144	enddo
2145
2146	if (k1.eq.0 .or. k2.eq.0) then
2147	write(,) 'PB! k1, k2 = ',k1,k2
2148	write(,) 'l,play(l) = ',l,play(l)/100
2149	do k = 1, nlev_toga-1
2150	write(,) 'k,plev_prof(k) = ',k,plev_prof(k)/100
2151	enddo
2152	endif
2153
2154	frac = (plev_prof(k2)-play(l))/(plev_prof(k2)-plev_prof(k1))
2155	t_mod(l)= t_prof(k2) - frac*(t_prof(k2)-t_prof(k1))
2156	q_mod(l)= q_prof(k2) - frac*(q_prof(k2)-q_prof(k1))
2157	u_mod(l)= u_prof(k2) - frac*(u_prof(k2)-u_prof(k1))
2158	v_mod(l)= v_prof(k2) - frac*(v_prof(k2)-v_prof(k1))
2159	w_mod(l)= w_prof(k2) - frac*(w_prof(k2)-w_prof(k1))
2160	ht_mod(l)= ht_prof(k2) - frac*(ht_prof(k2)-ht_prof(k1))
2161	vt_mod(l)= vt_prof(k2) - frac*(vt_prof(k2)-vt_prof(k1))
2162	hq_mod(l)= hq_prof(k2) - frac*(hq_prof(k2)-hq_prof(k1))
2163	vq_mod(l)= vq_prof(k2) - frac*(vq_prof(k2)-vq_prof(k1))
2164
2165	else !play>plev_prof(1)
2166
2167	k1=1
2168	k2=2
2169	frac1 = (play(l)-plev_prof(k2))/(plev_prof(k1)-plev_prof(k2))
2170	frac2 = (play(l)-plev_prof(k1))/(plev_prof(k1)-plev_prof(k2))
2171	t_mod(l)= frac1t_prof(k1) - frac2t_prof(k2)
2172	q_mod(l)= frac1q_prof(k1) - frac2q_prof(k2)
2173	u_mod(l)= frac1u_prof(k1) - frac2u_prof(k2)
2174	v_mod(l)= frac1v_prof(k1) - frac2v_prof(k2)
2175	w_mod(l)= frac1w_prof(k1) - frac2w_prof(k2)
2176	ht_mod(l)= frac1ht_prof(k1) - frac2ht_prof(k2)
2177	vt_mod(l)= frac1vt_prof(k1) - frac2vt_prof(k2)
2178	hq_mod(l)= frac1hq_prof(k1) - frac2hq_prof(k2)
2179	vq_mod(l)= frac1vq_prof(k1) - frac2vq_prof(k2)
2180
2181	endif ! play.le.plev_prof(1)
2182
2183	else ! above max altitude of forcing file
2184
2185	cjyg
2186	fact=20.*(plev_prof(nlev_toga)-play(l))/plev_prof(nlev_toga) !jyg
2187	fact = max(fact,0.) !jyg
2188	fact = exp(-fact) !jyg
2189	t_mod(l)= t_prof(nlev_toga) !jyg
2190	q_mod(l)= q_prof(nlev_toga)*fact !jyg
2191	u_mod(l)= u_prof(nlev_toga)*fact !jyg
2192	v_mod(l)= v_prof(nlev_toga)*fact !jyg
2193	w_mod(l)= 0.0 !jyg
2194	ht_mod(l)= ht_prof(nlev_toga) !jyg
2195	vt_mod(l)= vt_prof(nlev_toga) !jyg
2196	hq_mod(l)= hq_prof(nlev_toga)*fact !jyg
2197	vq_mod(l)= vq_prof(nlev_toga)*fact !jyg
2198
2199	endif ! play
2200
2201	enddo ! l
2202
2203	! do l = 1,llm
2204	! print *,'t_mod(l),q_mod(l),ht_mod(l),hq_mod(l) ',
2205	! $ l,t_mod(l),q_mod(l),ht_mod(l),hq_mod(l)
2206	! enddo
2207
2208	return
2209	end
2210
2211	!======================================================================
2212	SUBROUTINE interp_toga_time(day,day1,annee_ref
2213	i ,year_ini_toga,day_ini_toga,nt_toga,dt_toga,nlev_toga
2214	i ,ts_toga,plev_toga,t_toga,q_toga,u_toga,v_toga,w_toga
2215	i ,ht_toga,vt_toga,hq_toga,vq_toga
2216	o ,ts_prof,plev_prof,t_prof,q_prof,u_prof,v_prof,w_prof
2217	o ,ht_prof,vt_prof,hq_prof,vq_prof)
2218	implicit none
2219
2220	!---------------------------------------------------------------------------------------
2221	! Time interpolation of a 2D field to the timestep corresponding to day
2222	!
2223	! day: current julian day (e.g. 717538.2)
2224	! day1: first day of the simulation
2225	! nt_toga: total nb of data in the forcing (e.g. 480 for TOGA-COARE)
2226	! dt_toga: total time interval (in sec) between 2 forcing data (e.g. 6h for TOGA-COARE)
2227	!---------------------------------------------------------------------------------------
2228
2229	#include "compar1d.h"
2230
2231	! inputs:
2232	integer annee_ref
2233	integer nt_toga,nlev_toga
2234	integer year_ini_toga
2235	real day, day1,day_ini_toga,dt_toga
2236	real ts_toga(nt_toga)
2237	real plev_toga(nlev_toga,nt_toga),t_toga(nlev_toga,nt_toga)
2238	real q_toga(nlev_toga,nt_toga),u_toga(nlev_toga,nt_toga)
2239	real v_toga(nlev_toga,nt_toga),w_toga(nlev_toga,nt_toga)
2240	real ht_toga(nlev_toga,nt_toga),vt_toga(nlev_toga,nt_toga)
2241	real hq_toga(nlev_toga,nt_toga),vq_toga(nlev_toga,nt_toga)
2242	! outputs:
2243	real ts_prof
2244	real plev_prof(nlev_toga),t_prof(nlev_toga)
2245	real q_prof(nlev_toga),u_prof(nlev_toga)
2246	real v_prof(nlev_toga),w_prof(nlev_toga)
2247	real ht_prof(nlev_toga),vt_prof(nlev_toga)
2248	real hq_prof(nlev_toga),vq_prof(nlev_toga)
2249	! local:
2250	integer it_toga1, it_toga2,k
2251	real timeit,time_toga1,time_toga2,frac
2252
2253
2254	if (forcing_type.eq.2) then
2255	! Check that initial day of the simulation consistent with TOGA-COARE period:
2256	if (annee_ref.ne.1992 .and. annee_ref.ne.1993) then
2257	print*,'Pour TOGA-COARE, annee_ref doit etre 1992 ou 1993'
2258	print*,'Changer annee_ref dans run.def'
2259	stop
2260	endif
2261	if (annee_ref.eq.1992 .and. day1.lt.day_ini_toga) then
2262	print*,'TOGA-COARE a débuté le 1er Nov 1992 (jour julien=306)'
2263	print*,'Changer dayref dans run.def'
2264	stop
2265	endif
2266	if (annee_ref.eq.1993 .and. day1.gt.day_ini_toga+119) then
2267	print*,'TOGA-COARE a fini le 28 Feb 1993 (jour julien=59)'
2268	print*,'Changer dayref ou nday dans run.def'
2269	stop
2270	endif
2271
2272	else if (forcing_type.eq.4) then
2273
2274	! Check that initial day of the simulation consistent with TWP-ICE period:
2275	if (annee_ref.ne.2006) then
2276	print*,'Pour TWP-ICE, annee_ref doit etre 2006'
2277	print*,'Changer annee_ref dans run.def'
2278	stop
2279	endif
2280	if (annee_ref.eq.2006 .and. day1.lt.day_ini_toga) then
2281	print*,'TWP-ICE a debute le 17 Jan 2006 (jour julien=17)'
2282	print*,'Changer dayref dans run.def'
2283	stop
2284	endif
2285	if (annee_ref.eq.2006 .and. day1.gt.day_ini_toga+26) then
2286	print*,'TWP-ICE a fini le 12 Feb 2006 (jour julien=43)'
2287	print*,'Changer dayref ou nday dans run.def'
2288	stop
2289	endif
2290
2291	endif
2292
2293	! Determine timestep relative to the 1st day of TOGA-COARE:
2294	! timeit=(day-day1)*86400.
2295	! if (annee_ref.eq.1992) then
2296	! timeit=(day-day_ini_toga)*86400.
2297	! else
2298	! timeit=(day+61.-1.)*86400. ! 61 days between Nov01 and Dec31 1992
2299	! endif
2300	timeit=(day-day_ini_toga)*86400
2301
2302	! Determine the closest observation times:
2303	it_toga1=INT(timeit/dt_toga)+1
2304	it_toga2=it_toga1 + 1
2305	time_toga1=(it_toga1-1)*dt_toga
2306	time_toga2=(it_toga2-1)*dt_toga
2307
2308	if (it_toga1 .ge. nt_toga) then
2309	write(,) 'PB-stop: day, it_toga1, it_toga2, timeit: '
2310	: ,day,it_toga1,it_toga2,timeit/86400.
2311	stop
2312	endif
2313
2314	! time interpolation:
2315	frac=(time_toga2-timeit)/(time_toga2-time_toga1)
2316	frac=max(frac,0.0)
2317
2318	ts_prof = ts_toga(it_toga2)
2319	: -frac*(ts_toga(it_toga2)-ts_toga(it_toga1))
2320
2321	! print*,
2322	! :'day,annee_ref,day_ini_toga,timeit,it_toga1,it_toga2,SST:',
2323	! :day,annee_ref,day_ini_toga,timeit/86400.,it_toga1,it_toga2,ts_prof
2324
2325	do k=1,nlev_toga
2326	plev_prof(k) = 100.*(plev_toga(k,it_toga2)
2327	: -frac*(plev_toga(k,it_toga2)-plev_toga(k,it_toga1)))
2328	t_prof(k) = t_toga(k,it_toga2)
2329	: -frac*(t_toga(k,it_toga2)-t_toga(k,it_toga1))
2330	q_prof(k) = q_toga(k,it_toga2)
2331	: -frac*(q_toga(k,it_toga2)-q_toga(k,it_toga1))
2332	u_prof(k) = u_toga(k,it_toga2)
2333	: -frac*(u_toga(k,it_toga2)-u_toga(k,it_toga1))
2334	v_prof(k) = v_toga(k,it_toga2)
2335	: -frac*(v_toga(k,it_toga2)-v_toga(k,it_toga1))
2336	w_prof(k) = w_toga(k,it_toga2)
2337	: -frac*(w_toga(k,it_toga2)-w_toga(k,it_toga1))
2338	ht_prof(k) = ht_toga(k,it_toga2)
2339	: -frac*(ht_toga(k,it_toga2)-ht_toga(k,it_toga1))
2340	vt_prof(k) = vt_toga(k,it_toga2)
2341	: -frac*(vt_toga(k,it_toga2)-vt_toga(k,it_toga1))
2342	hq_prof(k) = hq_toga(k,it_toga2)
2343	: -frac*(hq_toga(k,it_toga2)-hq_toga(k,it_toga1))
2344	vq_prof(k) = vq_toga(k,it_toga2)
2345	: -frac*(vq_toga(k,it_toga2)-vq_toga(k,it_toga1))
2346	enddo
2347
2348	return
2349	END
2350
2351	!======================================================================
2352	SUBROUTINE interp_armcu_time(day,day1,annee_ref
2353	i ,year_ini_armcu,day_ini_armcu,nt_armcu,dt_armcu
2354	i ,nlev_armcu,fs_armcu,fl_armcu,at_armcu,rt_armcu
2355	i ,aqt_armcu,fs_prof,fl_prof,at_prof,rt_prof,aqt_prof)
2356	implicit none
2357
2358	!---------------------------------------------------------------------------------------
2359	! Time interpolation of a 2D field to the timestep corresponding to day
2360	!
2361	! day: current julian day (e.g. 717538.2)
2362	! day1: first day of the simulation
2363	! nt_armcu: total nb of data in the forcing (e.g. 31 for armcu)
2364	! dt_armcu: total time interval (in sec) between 2 forcing data (e.g. 1/2h for armcu)
2365	! fs= sensible flux
2366	! fl= latent flux
2367	! at,rt,aqt= advective and radiative tendencies
2368	!---------------------------------------------------------------------------------------
2369
2370	! inputs:
2371	integer annee_ref
2372	integer nt_armcu,nlev_armcu
2373	integer year_ini_armcu
2374	real day, day1,day_ini_armcu,dt_armcu
2375	real fs_armcu(nt_armcu),fl_armcu(nt_armcu),at_armcu(nt_armcu)
2376	real rt_armcu(nt_armcu),aqt_armcu(nt_armcu)
2377	! outputs:
2378	real fs_prof,fl_prof,at_prof,rt_prof,aqt_prof
2379	! local:
2380	integer it_armcu1, it_armcu2,k
2381	real timeit,time_armcu1,time_armcu2,frac
2382
2383	! Check that initial day of the simulation consistent with ARMCU period:
2384	if (annee_ref.ne.1997 ) then
2385	print*,'Pour ARMCU, annee_ref doit etre 1997 '
2386	print*,'Changer annee_ref dans run.def'
2387	stop
2388	endif
2389
2390	timeit=(day-day_ini_armcu)*86400
2391
2392	! Determine the closest observation times:
2393	it_armcu1=INT(timeit/dt_armcu)+1
2394	it_armcu2=it_armcu1 + 1
2395	time_armcu1=(it_armcu1-1)*dt_armcu
2396	time_armcu2=(it_armcu2-1)*dt_armcu
2397	print *,'timeit day day_ini_armcu',timeit,day,day_ini_armcu
2398	print *,'it_armcu1,it_armcu2,time_armcu1,time_armcu2',
2399	. it_armcu1,it_armcu2,time_armcu1,time_armcu2
2400
2401	if (it_armcu1 .ge. nt_armcu) then
2402	write(,) 'PB-stop: day, it_armcu1, it_armcu2, timeit: '
2403	: ,day,it_armcu1,it_armcu2,timeit/86400.
2404	stop
2405	endif
2406
2407	! time interpolation:
2408	frac=(time_armcu2-timeit)/(time_armcu2-time_armcu1)
2409	frac=max(frac,0.0)
2410
2411	fs_prof = fs_armcu(it_armcu2)
2412	: -frac*(fs_armcu(it_armcu2)-fs_armcu(it_armcu1))
2413	fl_prof = fl_armcu(it_armcu2)
2414	: -frac*(fl_armcu(it_armcu2)-fl_armcu(it_armcu1))
2415	at_prof = at_armcu(it_armcu2)
2416	: -frac*(at_armcu(it_armcu2)-at_armcu(it_armcu1))
2417	rt_prof = rt_armcu(it_armcu2)
2418	: -frac*(rt_armcu(it_armcu2)-rt_armcu(it_armcu1))
2419	aqt_prof = aqt_armcu(it_armcu2)
2420	: -frac*(aqt_armcu(it_armcu2)-aqt_armcu(it_armcu1))
2421
2422	print*,
2423	:'day,annee_ref,day_ini_armcu,timeit,it_armcu1,it_armcu2,SST:',
2424	:day,annee_ref,day_ini_armcu,timeit/86400.,it_armcu1,
2425	:it_armcu2,fs_prof,fl_prof,at_prof,rt_prof,aqt_prof
2426
2427	return
2428	END
2429
2430	!=====================================================================
2431	subroutine readprofiles(nlev_max,kmax,height,
2432	. thlprof,qtprof,uprof,
2433	. vprof,e12prof,ugprof,vgprof,
2434	. wfls,dqtdxls,dqtdyls,dqtdtls,
2435	. thlpcar)
2436	implicit none
2437
2438	integer nlev_max,kmax,kmax2
2439	logical :: llesread = .true.
2440
2441	real height(nlev_max),thlprof(nlev_max),qtprof(nlev_max),
2442	. uprof(nlev_max),vprof(nlev_max),e12prof(nlev_max),
2443	. ugprof(nlev_max),vgprof(nlev_max),wfls(nlev_max),
2444	. dqtdxls(nlev_max),dqtdyls(nlev_max),dqtdtls(nlev_max),
2445	. thlpcar(nlev_max)
2446
2447	integer, parameter :: ilesfile=1
2448	integer :: ierr,irad,imax,jtot,k
2449	logical :: lmoist,lcoriol,ltimedep
2450	real :: xsize,ysize
2451	real :: ustin,wsvsurf,timerad
2452	character(80) :: chmess
2453
2454	if(.not.(llesread)) return
2455
2456	open (ilesfile,file='prof.inp.001',status='old',iostat=ierr)
2457	if (ierr /= 0) stop 'ERROR:Prof.inp does not exist'
2458	read (ilesfile,*) kmax
2459	do k=1,kmax
2460	read (ilesfile,*) height(k),thlprof(k),qtprof (k),
2461	. uprof (k),vprof (k),e12prof(k)
2462	enddo
2463	close(ilesfile)
2464
2465	open(ilesfile,file='lscale.inp.001',status='old',iostat=ierr)
2466	if (ierr /= 0) stop 'ERROR:Lscale.inp does not exist'
2467	read (ilesfile,*) kmax2
2468	if (kmax .ne. kmax2) then
2469	print *, 'fichiers prof.inp et lscale.inp incompatibles :'
2470	print *, 'nbre de niveaux : ',kmax,' et ',kmax2
2471	stop 'lecture profiles'
2472	endif
2473	do k=1,kmax
2474	read (ilesfile,*) height(k),ugprof(k),vgprof(k),wfls(k),
2475	. dqtdxls(k),dqtdyls(k),dqtdtls(k),thlpcar(k)
2476	end do
2477	close(ilesfile)
2478
2479	return
2480	end
2481
2482
2483	!======================================================================
2484	subroutine readprofile_armcu(nlev_max,kmax,height,pprof,uprof,
2485	. vprof,thetaprof,tprof,qvprof,rvprof,aprof,bprof)
2486	!======================================================================
2487	implicit none
2488
2489	integer nlev_max,kmax,kmax2
2490	logical :: llesread = .true.
2491
2492	real height(nlev_max),pprof(nlev_max),tprof(nlev_max),
2493	. thetaprof(nlev_max),rvprof(nlev_max),
2494	. qvprof(nlev_max),uprof(nlev_max),vprof(nlev_max),
2495	. aprof(nlev_max+1),bprof(nlev_max+1)
2496
2497	integer, parameter :: ilesfile=1
2498	integer, parameter :: ifile=2
2499	integer :: ierr,irad,imax,jtot,k
2500	logical :: lmoist,lcoriol,ltimedep
2501	real :: xsize,ysize
2502	real :: ustin,wsvsurf,timerad
2503	character(80) :: chmess
2504
2505	if(.not.(llesread)) return
2506
2507	! Read profiles at full levels
2508	IF(nlev_max.EQ.19) THEN
2509	open (ilesfile,file='prof.inp.19',status='old',iostat=ierr)
2510	print *,'On ouvre prof.inp.19'
2511	ELSE
2512	open (ilesfile,file='prof.inp.40',status='old',iostat=ierr)
2513	print *,'On ouvre prof.inp.40'
2514	ENDIF
2515	if (ierr /= 0) stop 'ERROR:Prof.inp does not exist'
2516	read (ilesfile,*) kmax
2517	do k=1,kmax
2518	read (ilesfile,*) height(k) ,pprof(k), uprof(k), vprof(k),
2519	. thetaprof(k) ,tprof(k), qvprof(k),rvprof(k)
2520	enddo
2521	close(ilesfile)
2522
2523	! Vertical coordinates half levels for eta-coordinates (plev = alpha + beta * psurf)
2524	IF(nlev_max.EQ.19) THEN
2525	open (ifile,file='proh.inp.19',status='old',iostat=ierr)
2526	print *,'On ouvre proh.inp.19'
2527	if (ierr /= 0) stop 'ERROR:Proh.inp.19 does not exist'
2528	ELSE
2529	open (ifile,file='proh.inp.40',status='old',iostat=ierr)
2530	print *,'On ouvre proh.inp.40'
2531	if (ierr /= 0) stop 'ERROR:Proh.inp.40 does not exist'
2532	ENDIF
2533	read (ifile,*) kmax
2534	do k=1,kmax
2535	read (ifile,*) jtot,aprof(k),bprof(k)
2536	enddo
2537	close(ifile)
2538
2539	return
2540	end
2541	!===============================================================
2542	function ismin(n,sx,incx)
2543
2544	implicit none
2545	integer n,i,incx,ismin,ix
2546	real sx((n-1)*incx+1),sxmin
2547
2548	ix=1
2549	ismin=1
2550	sxmin=sx(1)
2551	do i=1,n-1
2552	ix=ix+incx
2553	if(sx(ix).lt.sxmin) then
2554	sxmin=sx(ix)
2555	ismin=i+1
2556	endif
2557	enddo
2558
2559	return
2560	end
2561
2562	!===============================================================
2563	function ismax(n,sx,incx)
2564
2565	implicit none
2566	integer n,i,incx,ismax,ix
2567	real sx((n-1)*incx+1),sxmax
2568
2569	ix=1
2570	ismax=1
2571	sxmax=sx(1)
2572	do i=1,n-1
2573	ix=ix+incx
2574	if(sx(ix).gt.sxmax) then
2575	sxmax=sx(ix)
2576	ismax=i+1
2577	endif
2578	enddo
2579
2580	return
2581	end
2582

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

Download in other formats:

Original Format