Context Navigation

1DUTILS.h_with_writelim_old @ 1795

Last change on this file since 1795 was 1795, checked in by Ehouarn Millour, 11 years ago

Version testing basee sur la r1794

Testing release based on r1794

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

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

Download in other formats:

Original Format