Context Navigation

1DUTILS.h_with_writelim @ 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: 79.7 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 fq_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	! fq_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	fq_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	real plev_min
1476
1477	plev_min = 55. ! pas de tendance de vap. d eau au-dessus de 55 hPa
1478
1479	open(21,file=trim(fich_toga),form='formatted')
1480	read(21,'(a)')
1481	do ip = 1, nt_toga
1482	read(21,'(a)')
1483	read(21,'(a)')
1484	read(21,223) iy, im, id, ih, bid, ts_toga(ip), bid,bid,bid,bid
1485	read(21,'(a)')
1486	read(21,'(a)')
1487
1488	do k = 1, nlev_toga
1489	read(21,230) plev_toga(k,ip), t_toga(k,ip), q_toga(k,ip)
1490	: ,u_toga(k,ip), v_toga(k,ip), w_toga(k,ip)
1491	: ,ht_toga(k,ip), vt_toga(k,ip), hq_toga(k,ip), vq_toga(k,ip)
1492
1493	! conversion in SI units:
1494	t_toga(k,ip)=t_toga(k,ip)+273.15 ! K
1495	q_toga(k,ip)=q_toga(k,ip)*0.001 ! kg/kg
1496	w_toga(k,ip)=w_toga(k,ip)*100./3600. ! Pa/s
1497	! no water vapour tendency above 55 hPa
1498	if (plev_toga(k,ip) .lt. plev_min) then
1499	q_toga(k,ip) = 0.
1500	hq_toga(k,ip) = 0.
1501	vq_toga(k,ip) =0.
1502	endif
1503	enddo
1504
1505	ts_toga(ip)=ts_toga(ip)+273.15 ! K
1506	enddo
1507	close(21)
1508
1509	223 format(4i3,6f8.2)
1510	226 format(f7.1,1x,10f8.2)
1511	227 format(f7.1,1x,1p,4e11.3)
1512	230 format(6f9.3,4e11.3)
1513
1514	return
1515	end
1516
1517	!=====================================================================
1518	subroutine read_twpice(fich_twpice,nlevel,ntime
1519	: ,T_srf,plev,T,q,u,v,omega
1520	: ,T_adv_h,T_adv_v,q_adv_h,q_adv_v)
1521
1522	!program reading forcings of the TWP-ICE experiment
1523
1524	! use netcdf
1525
1526	implicit none
1527
1528	#include "netcdf.inc"
1529
1530	integer ntime,nlevel
1531	integer l,k
1532	character*80 :: fich_twpice
1533	real*8 time(ntime)
1534	real*8 lat, lon, alt, phis
1535	real*8 lev(nlevel)
1536	real*8 plev(nlevel,ntime)
1537
1538	real*8 T(nlevel,ntime)
1539	real*8 q(nlevel,ntime),u(nlevel,ntime)
1540	real*8 v(nlevel,ntime)
1541	real*8 omega(nlevel,ntime), div(nlevel,ntime)
1542	real*8 T_adv_h(nlevel,ntime)
1543	real*8 T_adv_v(nlevel,ntime), q_adv_h(nlevel,ntime)
1544	real*8 q_adv_v(nlevel,ntime)
1545	real*8 s(nlevel,ntime), s_adv_h(nlevel,ntime)
1546	real*8 s_adv_v(nlevel,ntime)
1547	real*8 p_srf_aver(ntime), p_srf_center(ntime)
1548	real*8 T_srf(ntime)
1549
1550	integer nid, ierr
1551	integer nbvar3d
1552	parameter(nbvar3d=20)
1553	integer var3didin(nbvar3d)
1554
1555	ierr = NF_OPEN(fich_twpice,NF_NOWRITE,nid)
1556	if (ierr.NE.NF_NOERR) then
1557	write(,) 'ERROR: Pb opening forcings cdf file '
1558	write(,) NF_STRERROR(ierr)
1559	stop ""
1560	endif
1561
1562	ierr=NF_INQ_VARID(nid,"lat",var3didin(1))
1563	if(ierr/=NF_NOERR) then
1564	write(,) NF_STRERROR(ierr)
1565	stop 'lat'
1566	endif
1567
1568	ierr=NF_INQ_VARID(nid,"lon",var3didin(2))
1569	if(ierr/=NF_NOERR) then
1570	write(,) NF_STRERROR(ierr)
1571	stop 'lon'
1572	endif
1573
1574	ierr=NF_INQ_VARID(nid,"alt",var3didin(3))
1575	if(ierr/=NF_NOERR) then
1576	write(,) NF_STRERROR(ierr)
1577	stop 'alt'
1578	endif
1579
1580	ierr=NF_INQ_VARID(nid,"phis",var3didin(4))
1581	if(ierr/=NF_NOERR) then
1582	write(,) NF_STRERROR(ierr)
1583	stop 'phis'
1584	endif
1585
1586	ierr=NF_INQ_VARID(nid,"T",var3didin(5))
1587	if(ierr/=NF_NOERR) then
1588	write(,) NF_STRERROR(ierr)
1589	stop 'T'
1590	endif
1591
1592	ierr=NF_INQ_VARID(nid,"q",var3didin(6))
1593	if(ierr/=NF_NOERR) then
1594	write(,) NF_STRERROR(ierr)
1595	stop 'q'
1596	endif
1597
1598	ierr=NF_INQ_VARID(nid,"u",var3didin(7))
1599	if(ierr/=NF_NOERR) then
1600	write(,) NF_STRERROR(ierr)
1601	stop 'u'
1602	endif
1603
1604	ierr=NF_INQ_VARID(nid,"v",var3didin(8))
1605	if(ierr/=NF_NOERR) then
1606	write(,) NF_STRERROR(ierr)
1607	stop 'v'
1608	endif
1609
1610	ierr=NF_INQ_VARID(nid,"omega",var3didin(9))
1611	if(ierr/=NF_NOERR) then
1612	write(,) NF_STRERROR(ierr)
1613	stop 'omega'
1614	endif
1615
1616	ierr=NF_INQ_VARID(nid,"div",var3didin(10))
1617	if(ierr/=NF_NOERR) then
1618	write(,) NF_STRERROR(ierr)
1619	stop 'div'
1620	endif
1621
1622	ierr=NF_INQ_VARID(nid,"T_adv_h",var3didin(11))
1623	if(ierr/=NF_NOERR) then
1624	write(,) NF_STRERROR(ierr)
1625	stop 'T_adv_h'
1626	endif
1627
1628	ierr=NF_INQ_VARID(nid,"T_adv_v",var3didin(12))
1629	if(ierr/=NF_NOERR) then
1630	write(,) NF_STRERROR(ierr)
1631	stop 'T_adv_v'
1632	endif
1633
1634	ierr=NF_INQ_VARID(nid,"q_adv_h",var3didin(13))
1635	if(ierr/=NF_NOERR) then
1636	write(,) NF_STRERROR(ierr)
1637	stop 'q_adv_h'
1638	endif
1639
1640	ierr=NF_INQ_VARID(nid,"q_adv_v",var3didin(14))
1641	if(ierr/=NF_NOERR) then
1642	write(,) NF_STRERROR(ierr)
1643	stop 'q_adv_v'
1644	endif
1645
1646	ierr=NF_INQ_VARID(nid,"s",var3didin(15))
1647	if(ierr/=NF_NOERR) then
1648	write(,) NF_STRERROR(ierr)
1649	stop 's'
1650	endif
1651
1652	ierr=NF_INQ_VARID(nid,"s_adv_h",var3didin(16))
1653	if(ierr/=NF_NOERR) then
1654	write(,) NF_STRERROR(ierr)
1655	stop 's_adv_h'
1656	endif
1657
1658	ierr=NF_INQ_VARID(nid,"s_adv_v",var3didin(17))
1659	if(ierr/=NF_NOERR) then
1660	write(,) NF_STRERROR(ierr)
1661	stop 's_adv_v'
1662	endif
1663
1664	ierr=NF_INQ_VARID(nid,"p_srf_aver",var3didin(18))
1665	if(ierr/=NF_NOERR) then
1666	write(,) NF_STRERROR(ierr)
1667	stop 'p_srf_aver'
1668	endif
1669
1670	ierr=NF_INQ_VARID(nid,"p_srf_center",var3didin(19))
1671	if(ierr/=NF_NOERR) then
1672	write(,) NF_STRERROR(ierr)
1673	stop 'p_srf_center'
1674	endif
1675
1676	ierr=NF_INQ_VARID(nid,"T_srf",var3didin(20))
1677	if(ierr/=NF_NOERR) then
1678	write(,) NF_STRERROR(ierr)
1679	stop 'T_srf'
1680	endif
1681
1682	!dimensions lecture
1683	call catchaxis(nid,ntime,nlevel,time,lev,ierr)
1684
1685	!pressure
1686	do l=1,ntime
1687	do k=1,nlevel
1688	plev(k,l)=lev(k)
1689	enddo
1690	enddo
1691
1692	#ifdef NC_DOUBLE
1693	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(1),lat)
1694	#else
1695	ierr = NF_GET_VAR_REAL(nid,var3didin(1),lat)
1696	#endif
1697	if(ierr/=NF_NOERR) then
1698	write(,) NF_STRERROR(ierr)
1699	stop "getvarup"
1700	endif
1701	! write(,)'lecture lat ok',lat
1702
1703	#ifdef NC_DOUBLE
1704	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(2),lon)
1705	#else
1706	ierr = NF_GET_VAR_REAL(nid,var3didin(2),lon)
1707	#endif
1708	if(ierr/=NF_NOERR) then
1709	write(,) NF_STRERROR(ierr)
1710	stop "getvarup"
1711	endif
1712	! write(,)'lecture lon ok',lon
1713
1714	#ifdef NC_DOUBLE
1715	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(3),alt)
1716	#else
1717	ierr = NF_GET_VAR_REAL(nid,var3didin(3),alt)
1718	#endif
1719	if(ierr/=NF_NOERR) then
1720	write(,) NF_STRERROR(ierr)
1721	stop "getvarup"
1722	endif
1723	! write(,)'lecture alt ok',alt
1724
1725	#ifdef NC_DOUBLE
1726	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(4),phis)
1727	#else
1728	ierr = NF_GET_VAR_REAL(nid,var3didin(4),phis)
1729	#endif
1730	if(ierr/=NF_NOERR) then
1731	write(,) NF_STRERROR(ierr)
1732	stop "getvarup"
1733	endif
1734	! write(,)'lecture phis ok',phis
1735
1736	#ifdef NC_DOUBLE
1737	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(5),T)
1738	#else
1739	ierr = NF_GET_VAR_REAL(nid,var3didin(5),T)
1740	#endif
1741	if(ierr/=NF_NOERR) then
1742	write(,) NF_STRERROR(ierr)
1743	stop "getvarup"
1744	endif
1745	! write(,)'lecture T ok'
1746
1747	#ifdef NC_DOUBLE
1748	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(6),q)
1749	#else
1750	ierr = NF_GET_VAR_REAL(nid,var3didin(6),q)
1751	#endif
1752	if(ierr/=NF_NOERR) then
1753	write(,) NF_STRERROR(ierr)
1754	stop "getvarup"
1755	endif
1756	! write(,)'lecture q ok'
1757	!q in kg/kg
1758	do l=1,ntime
1759	do k=1,nlevel
1760	q(k,l)=q(k,l)/1000.
1761	enddo
1762	enddo
1763	#ifdef NC_DOUBLE
1764	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(7),u)
1765	#else
1766	ierr = NF_GET_VAR_REAL(nid,var3didin(7),u)
1767	#endif
1768	if(ierr/=NF_NOERR) then
1769	write(,) NF_STRERROR(ierr)
1770	stop "getvarup"
1771	endif
1772	! write(,)'lecture u ok'
1773
1774	#ifdef NC_DOUBLE
1775	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(8),v)
1776	#else
1777	ierr = NF_GET_VAR_REAL(nid,var3didin(8),v)
1778	#endif
1779	if(ierr/=NF_NOERR) then
1780	write(,) NF_STRERROR(ierr)
1781	stop "getvarup"
1782	endif
1783	! write(,)'lecture v ok'
1784
1785	#ifdef NC_DOUBLE
1786	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(9),omega)
1787	#else
1788	ierr = NF_GET_VAR_REAL(nid,var3didin(9),omega)
1789	#endif
1790	if(ierr/=NF_NOERR) then
1791	write(,) NF_STRERROR(ierr)
1792	stop "getvarup"
1793	endif
1794	! write(,)'lecture omega ok'
1795	!omega in mb/hour
1796	do l=1,ntime
1797	do k=1,nlevel
1798	omega(k,l)=omega(k,l)*100./3600.
1799	enddo
1800	enddo
1801
1802	#ifdef NC_DOUBLE
1803	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(10),div)
1804	#else
1805	ierr = NF_GET_VAR_REAL(nid,var3didin(10),div)
1806	#endif
1807	if(ierr/=NF_NOERR) then
1808	write(,) NF_STRERROR(ierr)
1809	stop "getvarup"
1810	endif
1811	! write(,)'lecture div ok'
1812
1813	#ifdef NC_DOUBLE
1814	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(11),T_adv_h)
1815	#else
1816	ierr = NF_GET_VAR_REAL(nid,var3didin(11),T_adv_h)
1817	#endif
1818	if(ierr/=NF_NOERR) then
1819	write(,) NF_STRERROR(ierr)
1820	stop "getvarup"
1821	endif
1822	! write(,)'lecture T_adv_h ok'
1823	!T adv in K/s
1824	do l=1,ntime
1825	do k=1,nlevel
1826	T_adv_h(k,l)=T_adv_h(k,l)/3600.
1827	enddo
1828	enddo
1829
1830
1831	#ifdef NC_DOUBLE
1832	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(12),T_adv_v)
1833	#else
1834	ierr = NF_GET_VAR_REAL(nid,var3didin(12),T_adv_v)
1835	#endif
1836	if(ierr/=NF_NOERR) then
1837	write(,) NF_STRERROR(ierr)
1838	stop "getvarup"
1839	endif
1840	! write(,)'lecture T_adv_v ok'
1841	!T adv in K/s
1842	do l=1,ntime
1843	do k=1,nlevel
1844	T_adv_v(k,l)=T_adv_v(k,l)/3600.
1845	enddo
1846	enddo
1847
1848	#ifdef NC_DOUBLE
1849	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(13),q_adv_h)
1850	#else
1851	ierr = NF_GET_VAR_REAL(nid,var3didin(13),q_adv_h)
1852	#endif
1853	if(ierr/=NF_NOERR) then
1854	write(,) NF_STRERROR(ierr)
1855	stop "getvarup"
1856	endif
1857	! write(,)'lecture q_adv_h ok'
1858	!q adv in kg/kg/s
1859	do l=1,ntime
1860	do k=1,nlevel
1861	q_adv_h(k,l)=q_adv_h(k,l)/1000./3600.
1862	enddo
1863	enddo
1864
1865
1866	#ifdef NC_DOUBLE
1867	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(14),q_adv_v)
1868	#else
1869	ierr = NF_GET_VAR_REAL(nid,var3didin(14),q_adv_v)
1870	#endif
1871	if(ierr/=NF_NOERR) then
1872	write(,) NF_STRERROR(ierr)
1873	stop "getvarup"
1874	endif
1875	! write(,)'lecture q_adv_v ok'
1876	!q adv in kg/kg/s
1877	do l=1,ntime
1878	do k=1,nlevel
1879	q_adv_v(k,l)=q_adv_v(k,l)/1000./3600.
1880	enddo
1881	enddo
1882
1883
1884	#ifdef NC_DOUBLE
1885	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(15),s)
1886	#else
1887	ierr = NF_GET_VAR_REAL(nid,var3didin(15),s)
1888	#endif
1889	if(ierr/=NF_NOERR) then
1890	write(,) NF_STRERROR(ierr)
1891	stop "getvarup"
1892	endif
1893
1894	#ifdef NC_DOUBLE
1895	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(16),s_adv_h)
1896	#else
1897	ierr = NF_GET_VAR_REAL(nid,var3didin(16),s_adv_h)
1898	#endif
1899	if(ierr/=NF_NOERR) then
1900	write(,) NF_STRERROR(ierr)
1901	stop "getvarup"
1902	endif
1903
1904	#ifdef NC_DOUBLE
1905	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(17),s_adv_v)
1906	#else
1907	ierr = NF_GET_VAR_REAL(nid,var3didin(17),s_adv_v)
1908	#endif
1909	if(ierr/=NF_NOERR) then
1910	write(,) NF_STRERROR(ierr)
1911	stop "getvarup"
1912	endif
1913
1914	#ifdef NC_DOUBLE
1915	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(18),p_srf_aver)
1916	#else
1917	ierr = NF_GET_VAR_REAL(nid,var3didin(18),p_srf_aver)
1918	#endif
1919	if(ierr/=NF_NOERR) then
1920	write(,) NF_STRERROR(ierr)
1921	stop "getvarup"
1922	endif
1923
1924	#ifdef NC_DOUBLE
1925	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(19),p_srf_center)
1926	#else
1927	ierr = NF_GET_VAR_REAL(nid,var3didin(19),p_srf_center)
1928	#endif
1929	if(ierr/=NF_NOERR) then
1930	write(,) NF_STRERROR(ierr)
1931	stop "getvarup"
1932	endif
1933
1934	#ifdef NC_DOUBLE
1935	ierr = NF_GET_VAR_DOUBLE(nid,var3didin(20),T_srf)
1936	#else
1937	ierr = NF_GET_VAR_REAL(nid,var3didin(20),T_srf)
1938	#endif
1939	if(ierr/=NF_NOERR) then
1940	write(,) NF_STRERROR(ierr)
1941	stop "getvarup"
1942	endif
1943	! write(,)'lecture T_srf ok', T_srf
1944
1945	return
1946	end subroutine read_twpice
1947	!=====================================================================
1948	subroutine catchaxis(nid,ttm,llm,time,lev,ierr)
1949
1950	! use netcdf
1951
1952	implicit none
1953	#include "netcdf.inc"
1954	integer nid,ttm,llm
1955	real*8 time(ttm)
1956	real*8 lev(llm)
1957	integer ierr
1958
1959	integer i
1960	integer timevar,levvar
1961	integer timelen,levlen
1962	integer timedimin,levdimin
1963
1964	! Control & lecture on dimensions
1965	! ===============================
1966	ierr=NF_INQ_DIMID(nid,"time",timedimin)
1967	ierr=NF_INQ_VARID(nid,"time",timevar)
1968	if (ierr.NE.NF_NOERR) then
1969	write(,) 'ERROR: Field <time> is missing'
1970	stop ""
1971	endif
1972	ierr=NF_INQ_DIMLEN(nid,timedimin,timelen)
1973
1974	ierr=NF_INQ_DIMID(nid,"lev",levdimin)
1975	ierr=NF_INQ_VARID(nid,"lev",levvar)
1976	if (ierr.NE.NF_NOERR) then
1977	write(,) 'ERROR: Field <lev> is lacking'
1978	stop ""
1979	endif
1980	ierr=NF_INQ_DIMLEN(nid,levdimin,levlen)
1981
1982	if((timelen/=ttm).or.(levlen/=llm)) then
1983	write(,) 'ERROR: Not the good lenght for axis'
1984	write(,) 'longitude: ',timelen,ttm+1
1985	write(,) 'latitude: ',levlen,llm
1986	stop ""
1987	endif
1988
1989	!#ifdef NC_DOUBLE
1990	ierr = NF_GET_VAR_DOUBLE(nid,timevar,time)
1991	ierr = NF_GET_VAR_DOUBLE(nid,levvar,lev)
1992	!#else
1993	! ierr = NF_GET_VAR_REAL(nid,timevar,time)
1994	! ierr = NF_GET_VAR_REAL(nid,levvar,lev)
1995	!#endif
1996
1997	return
1998	end
1999
2000	!======================================================================
2001	SUBROUTINE read_rico(fich_rico,nlev_rico,ps_rico,play
2002	: ,ts_rico,t_rico,q_rico,u_rico,v_rico,w_rico
2003	: ,dth_dyn,dqh_dyn)
2004	implicit none
2005
2006	c-------------------------------------------------------------------------
2007	c Read RICO forcing data
2008	c-------------------------------------------------------------------------
2009	#include "dimensions.h"
2010
2011
2012	integer nlev_rico
2013	real ts_rico,ps_rico
2014	real t_rico(llm),q_rico(llm)
2015	real u_rico(llm),v_rico(llm)
2016	real w_rico(llm)
2017	real dth_dyn(llm)
2018	real dqh_dyn(llm)
2019
2020
2021	real play(llm),zlay(llm)
2022
2023
2024	real prico(nlev_rico),zrico(nlev_rico)
2025
2026	character*80 fich_rico
2027
2028	integer k,l
2029
2030
2031	print*,fich_rico
2032	open(21,file=trim(fich_rico),form='formatted')
2033	do k=1,llm
2034	zlay(k)=0.
2035	enddo
2036
2037	read(21,*) ps_rico,ts_rico
2038	prico(1)=ps_rico
2039	zrico(1)=0.0
2040	do l=2,nlev_rico
2041	read(21,*) k,prico(l),zrico(l)
2042	enddo
2043	close(21)
2044
2045	do k=1,llm
2046	do l=1,80
2047	if(prico(l)>play(k)) then
2048	if(play(k)>prico(l+1)) then
2049	zlay(k)=zrico(l)+(play(k)-prico(l)) *
2050	& (zrico(l+1)-zrico(l))/(prico(l+1)-prico(l))
2051	else
2052	zlay(k)=zrico(l)+(play(k)-prico(80))*
2053	& (zrico(81)-zrico(80))/(prico(81)-prico(80))
2054	endif
2055	endif
2056	enddo
2057	print*,k,zlay(k)
2058	! U
2059	if(0 < zlay(k) .and. zlay(k) < 4000) then
2060	u_rico(k)=-9.9 + (-1.9 + 9.9)*zlay(k)/4000
2061	elseif(4000 < zlay(k) .and. zlay(k) < 12000) then
2062	u_rico(k)= -1.9 + (30.0 + 1.9) /
2063	: (12000 - 4000) * (zlay(k) - 4000)
2064	elseif(12000 < zlay(k) .and. zlay(k) < 13000) then
2065	u_rico(k)=30.0
2066	elseif(13000 < zlay(k) .and. zlay(k) < 20000) then
2067	u_rico(k)=30.0 - (30.0) /
2068	: (20000 - 13000) * (zlay(k) - 13000)
2069	else
2070	u_rico(k)=0.0
2071	endif
2072
2073	!Q_v
2074	if(0 < zlay(k) .and. zlay(k) < 740) then
2075	q_rico(k)=16.0 + (13.8 - 16.0) / (740) * zlay(k)
2076	elseif(740 < zlay(k) .and. zlay(k) < 3260) then
2077	q_rico(k)=13.8 + (2.4 - 13.8) /
2078	: (3260 - 740) * (zlay(k) - 740)
2079	elseif(3260 < zlay(k) .and. zlay(k) < 4000) then
2080	q_rico(k)=2.4 + (1.8 - 2.4) /
2081	: (4000 - 3260) * (zlay(k) - 3260)
2082	elseif(4000 < zlay(k) .and. zlay(k) < 9000) then
2083	q_rico(k)=1.8 + (0 - 1.8) /
2084	: (10000 - 4000) * (zlay(k) - 4000)
2085	else
2086	q_rico(k)=0.0
2087	endif
2088
2089	!T
2090	if(0 < zlay(k) .and. zlay(k) < 740) then
2091	t_rico(k)=299.2 + (292.0 - 299.2) / (740) * zlay(k)
2092	elseif(740 < zlay(k) .and. zlay(k) < 4000) then
2093	t_rico(k)=292.0 + (278.0 - 292.0) /
2094	: (4000 - 740) * (zlay(k) - 740)
2095	elseif(4000 < zlay(k) .and. zlay(k) < 15000) then
2096	t_rico(k)=278.0 + (203.0 - 278.0) /
2097	: (15000 - 4000) * (zlay(k) - 4000)
2098	elseif(15000 < zlay(k) .and. zlay(k) < 17500) then
2099	t_rico(k)=203.0 + (194.0 - 203.0) /
2100	: (17500 - 15000)* (zlay(k) - 15000)
2101	elseif(17500 < zlay(k) .and. zlay(k) < 20000) then
2102	t_rico(k)=194.0 + (206.0 - 194.0) /
2103	: (20000 - 17500)* (zlay(k) - 17500)
2104	elseif(20000 < zlay(k) .and. zlay(k) < 60000) then
2105	t_rico(k)=206.0 + (270.0 - 206.0) /
2106	: (60000 - 20000)* (zlay(k) - 20000)
2107	endif
2108
2109	! W
2110	if(0 < zlay(k) .and. zlay(k) < 2260 ) then
2111	w_rico(k)=- (0.005/2260) * zlay(k)
2112	elseif(2260 < zlay(k) .and. zlay(k) < 4000 ) then
2113	w_rico(k)=- 0.005
2114	elseif(4000 < zlay(k) .and. zlay(k) < 5000 ) then
2115	w_rico(k)=- 0.005 + (0.005/ (5000 - 4000)) * (zlay(k) - 4000)
2116	else
2117	w_rico(k)=0.0
2118	endif
2119
2120	! dThrz+dTsw0+dTlw0
2121	if(0 < zlay(k) .and. zlay(k) < 4000) then
2122	dth_dyn(k)=- 2.51 / 86400 + (-2.18 + 2.51 )/
2123	: (864004000) zlay(k)
2124	elseif(4000 < zlay(k) .and. zlay(k) < 5000) then
2125	dth_dyn(k)=- 2.18 / 86400 + ( 2.18 ) /
2126	: (86400(5000 - 4000)) (zlay(k) - 4000)
2127	else
2128	dth_dyn(k)=0.0
2129	endif
2130	! dQhrz
2131	if(0 < zlay(k) .and. zlay(k) < 3000) then
2132	dqh_dyn(k)=-1.0 / 86400 + (0.345 + 1.0)/
2133	: (864003000) (zlay(k))
2134	elseif(3000 < zlay(k) .and. zlay(k) < 4000) then
2135	dqh_dyn(k)=0.345 / 86400
2136	elseif(4000 < zlay(k) .and. zlay(k) < 5000) then
2137	dqh_dyn(k)=0.345 / 86400 +
2138	+ (-0.345)/(86400 * (5000 - 4000)) * (zlay(k)-4000)
2139	else
2140	dqh_dyn(k)=0.0
2141	endif
2142
2143	!? if(play(k)>6e4) then
2144	!? ratqs0(1,k)=ratqsbas*(plev(1)-play(k))/(plev(1)-6e4)
2145	!? elseif((play(k)>3e4).and.(play(k)<6e4)) then
2146	!? ratqs0(1,k)=ratqsbas+(ratqshaut-ratqsbas)&
2147	!? *(6e4-play(k))/(6e4-3e4)
2148	!? else
2149	!? ratqs0(1,k)=ratqshaut
2150	!? endif
2151
2152	enddo
2153
2154	do k=1,llm
2155	q_rico(k)=q_rico(k)/1e3
2156	dqh_dyn(k)=dqh_dyn(k)/1e3
2157	v_rico(k)=-3.8
2158	enddo
2159
2160	return
2161	end
2162
2163	!=====================================================================
2164	c-------------------------------------------------------------------------
2165	SUBROUTINE read_armcu(fich_armcu,nlev_armcu,nt_armcu,
2166	: sens,flat,adv_theta,rad_theta,adv_qt)
2167	implicit none
2168
2169	c-------------------------------------------------------------------------
2170	c Read ARM_CU case forcing data
2171	c-------------------------------------------------------------------------
2172
2173	integer nlev_armcu,nt_armcu
2174	real sens(nt_armcu),flat(nt_armcu)
2175	real adv_theta(nt_armcu),rad_theta(nt_armcu),adv_qt(nt_armcu)
2176	character*80 fich_armcu
2177
2178	integer no,l,k,ip
2179	real riy,rim,rid,rih,bid
2180
2181	integer iy,im,id,ih,in
2182
2183	open(21,file=trim(fich_armcu),form='formatted')
2184	read(21,'(a)')
2185	do ip = 1, nt_armcu
2186	read(21,'(a)')
2187	read(21,'(a)')
2188	read(21,223) iy, im, id, ih, in, sens(ip),flat(ip),
2189	: adv_theta(ip),rad_theta(ip),adv_qt(ip)
2190	print *,'forcages=',iy,im,id,ih,in, sens(ip),flat(ip),
2191	: adv_theta(ip),rad_theta(ip),adv_qt(ip)
2192	enddo
2193	close(21)
2194
2195	223 format(5i3,5f8.3)
2196	226 format(f7.1,1x,10f8.2)
2197	227 format(f7.1,1x,1p,4e11.3)
2198	230 format(6f9.3,4e11.3)
2199
2200	return
2201	end
2202
2203	!=====================================================================
2204	SUBROUTINE interp_toga_vertical(play,nlev_toga,plev_prof
2205	: ,t_prof,q_prof,u_prof,v_prof,w_prof
2206	: ,ht_prof,vt_prof,hq_prof,vq_prof
2207	: ,t_mod,q_mod,u_mod,v_mod,w_mod
2208	: ,ht_mod,vt_mod,hq_mod,vq_mod,mxcalc)
2209
2210	implicit none
2211
2212	#include "dimensions.h"
2213
2214	c-------------------------------------------------------------------------
2215	c Vertical interpolation of TOGA-COARE forcing data onto model levels
2216	c-------------------------------------------------------------------------
2217
2218	integer nlevmax
2219	parameter (nlevmax=41)
2220	integer nlev_toga,mxcalc
2221	! real play(llm), plev_prof(nlevmax)
2222	! real t_prof(nlevmax),q_prof(nlevmax)
2223	! real u_prof(nlevmax),v_prof(nlevmax), w_prof(nlevmax)
2224	! real ht_prof(nlevmax),vt_prof(nlevmax)
2225	! real hq_prof(nlevmax),vq_prof(nlevmax)
2226
2227	real play(llm), plev_prof(nlev_toga)
2228	real t_prof(nlev_toga),q_prof(nlev_toga)
2229	real u_prof(nlev_toga),v_prof(nlev_toga), w_prof(nlev_toga)
2230	real ht_prof(nlev_toga),vt_prof(nlev_toga)
2231	real hq_prof(nlev_toga),vq_prof(nlev_toga)
2232
2233	real t_mod(llm),q_mod(llm)
2234	real u_mod(llm),v_mod(llm), w_mod(llm)
2235	real ht_mod(llm),vt_mod(llm)
2236	real hq_mod(llm),vq_mod(llm)
2237
2238	integer l,k,k1,k2,kp
2239	real aa,frac,frac1,frac2,fact
2240
2241	do l = 1, llm
2242
2243	if (play(l).ge.plev_prof(nlev_toga)) then
2244
2245	mxcalc=l
2246	k1=0
2247	k2=0
2248
2249	if (play(l).le.plev_prof(1)) then
2250
2251	do k = 1, nlev_toga-1
2252	if (play(l).le.plev_prof(k)
2253	: .and. play(l).gt.plev_prof(k+1)) then
2254	k1=k
2255	k2=k+1
2256	endif
2257	enddo
2258
2259	if (k1.eq.0 .or. k2.eq.0) then
2260	write(,) 'PB! k1, k2 = ',k1,k2
2261	write(,) 'l,play(l) = ',l,play(l)/100
2262	do k = 1, nlev_toga-1
2263	write(,) 'k,plev_prof(k) = ',k,plev_prof(k)/100
2264	enddo
2265	endif
2266
2267	frac = (plev_prof(k2)-play(l))/(plev_prof(k2)-plev_prof(k1))
2268	t_mod(l)= t_prof(k2) - frac*(t_prof(k2)-t_prof(k1))
2269	q_mod(l)= q_prof(k2) - frac*(q_prof(k2)-q_prof(k1))
2270	u_mod(l)= u_prof(k2) - frac*(u_prof(k2)-u_prof(k1))
2271	v_mod(l)= v_prof(k2) - frac*(v_prof(k2)-v_prof(k1))
2272	w_mod(l)= w_prof(k2) - frac*(w_prof(k2)-w_prof(k1))
2273	ht_mod(l)= ht_prof(k2) - frac*(ht_prof(k2)-ht_prof(k1))
2274	vt_mod(l)= vt_prof(k2) - frac*(vt_prof(k2)-vt_prof(k1))
2275	hq_mod(l)= hq_prof(k2) - frac*(hq_prof(k2)-hq_prof(k1))
2276	vq_mod(l)= vq_prof(k2) - frac*(vq_prof(k2)-vq_prof(k1))
2277
2278	else !play>plev_prof(1)
2279
2280	k1=1
2281	k2=2
2282	frac1 = (play(l)-plev_prof(k2))/(plev_prof(k1)-plev_prof(k2))
2283	frac2 = (play(l)-plev_prof(k1))/(plev_prof(k1)-plev_prof(k2))
2284	t_mod(l)= frac1t_prof(k1) - frac2t_prof(k2)
2285	q_mod(l)= frac1q_prof(k1) - frac2q_prof(k2)
2286	u_mod(l)= frac1u_prof(k1) - frac2u_prof(k2)
2287	v_mod(l)= frac1v_prof(k1) - frac2v_prof(k2)
2288	w_mod(l)= frac1w_prof(k1) - frac2w_prof(k2)
2289	ht_mod(l)= frac1ht_prof(k1) - frac2ht_prof(k2)
2290	vt_mod(l)= frac1vt_prof(k1) - frac2vt_prof(k2)
2291	hq_mod(l)= frac1hq_prof(k1) - frac2hq_prof(k2)
2292	vq_mod(l)= frac1vq_prof(k1) - frac2vq_prof(k2)
2293
2294	endif ! play.le.plev_prof(1)
2295
2296	else ! above max altitude of forcing file
2297
2298	cjyg
2299	fact=20.*(plev_prof(nlev_toga)-play(l))/plev_prof(nlev_toga) !jyg
2300	fact = max(fact,0.) !jyg
2301	fact = exp(-fact) !jyg
2302	t_mod(l)= t_prof(nlev_toga) !jyg
2303	q_mod(l)= q_prof(nlev_toga)*fact !jyg
2304	u_mod(l)= u_prof(nlev_toga)*fact !jyg
2305	v_mod(l)= v_prof(nlev_toga)*fact !jyg
2306	w_mod(l)= 0.0 !jyg
2307	ht_mod(l)= ht_prof(nlev_toga) !jyg
2308	vt_mod(l)= vt_prof(nlev_toga) !jyg
2309	hq_mod(l)= hq_prof(nlev_toga)*fact !jyg
2310	vq_mod(l)= vq_prof(nlev_toga)*fact !jyg
2311
2312	endif ! play
2313
2314	enddo ! l
2315
2316	! do l = 1,llm
2317	! print *,'t_mod(l),q_mod(l),ht_mod(l),hq_mod(l) ',
2318	! $ l,t_mod(l),q_mod(l),ht_mod(l),hq_mod(l)
2319	! enddo
2320
2321	return
2322	end
2323
2324	!======================================================================
2325	SUBROUTINE interp_toga_time(day,day1,annee_ref
2326	i ,year_ini_toga,day_ini_toga,nt_toga,dt_toga,nlev_toga
2327	i ,ts_toga,plev_toga,t_toga,q_toga,u_toga,v_toga,w_toga
2328	i ,ht_toga,vt_toga,hq_toga,vq_toga
2329	o ,ts_prof,plev_prof,t_prof,q_prof,u_prof,v_prof,w_prof
2330	o ,ht_prof,vt_prof,hq_prof,vq_prof)
2331	implicit none
2332
2333	!---------------------------------------------------------------------------------------
2334	! Time interpolation of a 2D field to the timestep corresponding to day
2335	!
2336	! day: current julian day (e.g. 717538.2)
2337	! day1: first day of the simulation
2338	! nt_toga: total nb of data in the forcing (e.g. 480 for TOGA-COARE)
2339	! dt_toga: total time interval (in sec) between 2 forcing data (e.g. 6h for TOGA-COARE)
2340	!---------------------------------------------------------------------------------------
2341
2342	#include "compar1d.h"
2343
2344	! inputs:
2345	integer annee_ref
2346	integer nt_toga,nlev_toga
2347	integer year_ini_toga
2348	real day, day1,day_ini_toga,dt_toga
2349	real ts_toga(nt_toga)
2350	real plev_toga(nlev_toga,nt_toga),t_toga(nlev_toga,nt_toga)
2351	real q_toga(nlev_toga,nt_toga),u_toga(nlev_toga,nt_toga)
2352	real v_toga(nlev_toga,nt_toga),w_toga(nlev_toga,nt_toga)
2353	real ht_toga(nlev_toga,nt_toga),vt_toga(nlev_toga,nt_toga)
2354	real hq_toga(nlev_toga,nt_toga),vq_toga(nlev_toga,nt_toga)
2355	! outputs:
2356	real ts_prof
2357	real plev_prof(nlev_toga),t_prof(nlev_toga)
2358	real q_prof(nlev_toga),u_prof(nlev_toga)
2359	real v_prof(nlev_toga),w_prof(nlev_toga)
2360	real ht_prof(nlev_toga),vt_prof(nlev_toga)
2361	real hq_prof(nlev_toga),vq_prof(nlev_toga)
2362	! local:
2363	integer it_toga1, it_toga2,k
2364	real timeit,time_toga1,time_toga2,frac
2365
2366
2367	if (forcing_type.eq.2) then
2368	! Check that initial day of the simulation consistent with TOGA-COARE period:
2369	if (annee_ref.ne.1992 .and. annee_ref.ne.1993) then
2370	print*,'Pour TOGA-COARE, annee_ref doit etre 1992 ou 1993'
2371	print*,'Changer annee_ref dans run.def'
2372	stop
2373	endif
2374	if (annee_ref.eq.1992 .and. day1.lt.day_ini_toga) then
2375	print*,'TOGA-COARE a débuté le 1er Nov 1992 (jour julien=306)'
2376	print*,'Changer dayref dans run.def'
2377	stop
2378	endif
2379	if (annee_ref.eq.1993 .and. day1.gt.day_ini_toga+119) then
2380	print*,'TOGA-COARE a fini le 28 Feb 1993 (jour julien=59)'
2381	print*,'Changer dayref ou nday dans run.def'
2382	stop
2383	endif
2384
2385	else if (forcing_type.eq.4) then
2386
2387	! Check that initial day of the simulation consistent with TWP-ICE period:
2388	if (annee_ref.ne.2006) then
2389	print*,'Pour TWP-ICE, annee_ref doit etre 2006'
2390	print*,'Changer annee_ref dans run.def'
2391	stop
2392	endif
2393	if (annee_ref.eq.2006 .and. day1.lt.day_ini_toga) then
2394	print*,'TWP-ICE a debute le 17 Jan 2006 (jour julien=17)'
2395	print*,'Changer dayref dans run.def'
2396	stop
2397	endif
2398	if (annee_ref.eq.2006 .and. day1.gt.day_ini_toga+26) then
2399	print*,'TWP-ICE a fini le 12 Feb 2006 (jour julien=43)'
2400	print*,'Changer dayref ou nday dans run.def'
2401	stop
2402	endif
2403
2404	endif
2405
2406	! Determine timestep relative to the 1st day of TOGA-COARE:
2407	! timeit=(day-day1)*86400.
2408	! if (annee_ref.eq.1992) then
2409	! timeit=(day-day_ini_toga)*86400.
2410	! else
2411	! timeit=(day+61.-1.)*86400. ! 61 days between Nov01 and Dec31 1992
2412	! endif
2413	timeit=(day-day_ini_toga)*86400
2414
2415	! Determine the closest observation times:
2416	it_toga1=INT(timeit/dt_toga)+1
2417	it_toga2=it_toga1 + 1
2418	time_toga1=(it_toga1-1)*dt_toga
2419	time_toga2=(it_toga2-1)*dt_toga
2420
2421	if (it_toga1 .ge. nt_toga) then
2422	write(,) 'PB-stop: day, it_toga1, it_toga2, timeit: '
2423	: ,day,it_toga1,it_toga2,timeit/86400.
2424	stop
2425	endif
2426
2427	! time interpolation:
2428	frac=(time_toga2-timeit)/(time_toga2-time_toga1)
2429	frac=max(frac,0.0)
2430
2431	ts_prof = ts_toga(it_toga2)
2432	: -frac*(ts_toga(it_toga2)-ts_toga(it_toga1))
2433
2434	! print*,
2435	! :'day,annee_ref,day_ini_toga,timeit,it_toga1,it_toga2,SST:',
2436	! :day,annee_ref,day_ini_toga,timeit/86400.,it_toga1,it_toga2,ts_prof
2437
2438	do k=1,nlev_toga
2439	plev_prof(k) = 100.*(plev_toga(k,it_toga2)
2440	: -frac*(plev_toga(k,it_toga2)-plev_toga(k,it_toga1)))
2441	t_prof(k) = t_toga(k,it_toga2)
2442	: -frac*(t_toga(k,it_toga2)-t_toga(k,it_toga1))
2443	q_prof(k) = q_toga(k,it_toga2)
2444	: -frac*(q_toga(k,it_toga2)-q_toga(k,it_toga1))
2445	u_prof(k) = u_toga(k,it_toga2)
2446	: -frac*(u_toga(k,it_toga2)-u_toga(k,it_toga1))
2447	v_prof(k) = v_toga(k,it_toga2)
2448	: -frac*(v_toga(k,it_toga2)-v_toga(k,it_toga1))
2449	w_prof(k) = w_toga(k,it_toga2)
2450	: -frac*(w_toga(k,it_toga2)-w_toga(k,it_toga1))
2451	ht_prof(k) = ht_toga(k,it_toga2)
2452	: -frac*(ht_toga(k,it_toga2)-ht_toga(k,it_toga1))
2453	vt_prof(k) = vt_toga(k,it_toga2)
2454	: -frac*(vt_toga(k,it_toga2)-vt_toga(k,it_toga1))
2455	hq_prof(k) = hq_toga(k,it_toga2)
2456	: -frac*(hq_toga(k,it_toga2)-hq_toga(k,it_toga1))
2457	vq_prof(k) = vq_toga(k,it_toga2)
2458	: -frac*(vq_toga(k,it_toga2)-vq_toga(k,it_toga1))
2459	enddo
2460
2461	return
2462	END
2463
2464	!======================================================================
2465	SUBROUTINE interp_armcu_time(day,day1,annee_ref
2466	i ,year_ini_armcu,day_ini_armcu,nt_armcu,dt_armcu
2467	i ,nlev_armcu,fs_armcu,fl_armcu,at_armcu,rt_armcu
2468	i ,aqt_armcu,fs_prof,fl_prof,at_prof,rt_prof,aqt_prof)
2469	implicit none
2470
2471	!---------------------------------------------------------------------------------------
2472	! Time interpolation of a 2D field to the timestep corresponding to day
2473	!
2474	! day: current julian day (e.g. 717538.2)
2475	! day1: first day of the simulation
2476	! nt_armcu: total nb of data in the forcing (e.g. 31 for armcu)
2477	! dt_armcu: total time interval (in sec) between 2 forcing data (e.g. 1/2h for armcu)
2478	! fs= sensible flux
2479	! fl= latent flux
2480	! at,rt,aqt= advective and radiative tendencies
2481	!---------------------------------------------------------------------------------------
2482
2483	! inputs:
2484	integer annee_ref
2485	integer nt_armcu,nlev_armcu
2486	integer year_ini_armcu
2487	real day, day1,day_ini_armcu,dt_armcu
2488	real fs_armcu(nt_armcu),fl_armcu(nt_armcu),at_armcu(nt_armcu)
2489	real rt_armcu(nt_armcu),aqt_armcu(nt_armcu)
2490	! outputs:
2491	real fs_prof,fl_prof,at_prof,rt_prof,aqt_prof
2492	! local:
2493	integer it_armcu1, it_armcu2,k
2494	real timeit,time_armcu1,time_armcu2,frac
2495
2496	! Check that initial day of the simulation consistent with ARMCU period:
2497	if (annee_ref.ne.1997 ) then
2498	print*,'Pour ARMCU, annee_ref doit etre 1997 '
2499	print*,'Changer annee_ref dans run.def'
2500	stop
2501	endif
2502
2503	timeit=(day-day_ini_armcu)*86400
2504
2505	! Determine the closest observation times:
2506	it_armcu1=INT(timeit/dt_armcu)+1
2507	it_armcu2=it_armcu1 + 1
2508	time_armcu1=(it_armcu1-1)*dt_armcu
2509	time_armcu2=(it_armcu2-1)*dt_armcu
2510	print *,'timeit day day_ini_armcu',timeit,day,day_ini_armcu
2511	print *,'it_armcu1,it_armcu2,time_armcu1,time_armcu2',
2512	. it_armcu1,it_armcu2,time_armcu1,time_armcu2
2513
2514	if (it_armcu1 .ge. nt_armcu) then
2515	write(,) 'PB-stop: day, it_armcu1, it_armcu2, timeit: '
2516	: ,day,it_armcu1,it_armcu2,timeit/86400.
2517	stop
2518	endif
2519
2520	! time interpolation:
2521	frac=(time_armcu2-timeit)/(time_armcu2-time_armcu1)
2522	frac=max(frac,0.0)
2523
2524	fs_prof = fs_armcu(it_armcu2)
2525	: -frac*(fs_armcu(it_armcu2)-fs_armcu(it_armcu1))
2526	fl_prof = fl_armcu(it_armcu2)
2527	: -frac*(fl_armcu(it_armcu2)-fl_armcu(it_armcu1))
2528	at_prof = at_armcu(it_armcu2)
2529	: -frac*(at_armcu(it_armcu2)-at_armcu(it_armcu1))
2530	rt_prof = rt_armcu(it_armcu2)
2531	: -frac*(rt_armcu(it_armcu2)-rt_armcu(it_armcu1))
2532	aqt_prof = aqt_armcu(it_armcu2)
2533	: -frac*(aqt_armcu(it_armcu2)-aqt_armcu(it_armcu1))
2534
2535	print*,
2536	:'day,annee_ref,day_ini_armcu,timeit,it_armcu1,it_armcu2,SST:',
2537	:day,annee_ref,day_ini_armcu,timeit/86400.,it_armcu1,
2538	:it_armcu2,fs_prof,fl_prof,at_prof,rt_prof,aqt_prof
2539
2540	return
2541	END
2542
2543	!=====================================================================
2544	subroutine readprofiles(nlev_max,kmax,height,
2545	. thlprof,qtprof,uprof,
2546	. vprof,e12prof,ugprof,vgprof,
2547	. wfls,dqtdxls,dqtdyls,dqtdtls,
2548	. thlpcar)
2549	implicit none
2550
2551	integer nlev_max,kmax,kmax2
2552	logical :: llesread = .true.
2553
2554	real height(nlev_max),thlprof(nlev_max),qtprof(nlev_max),
2555	. uprof(nlev_max),vprof(nlev_max),e12prof(nlev_max),
2556	. ugprof(nlev_max),vgprof(nlev_max),wfls(nlev_max),
2557	. dqtdxls(nlev_max),dqtdyls(nlev_max),dqtdtls(nlev_max),
2558	. thlpcar(nlev_max)
2559
2560	integer, parameter :: ilesfile=1
2561	integer :: ierr,irad,imax,jtot,k
2562	logical :: lmoist,lcoriol,ltimedep
2563	real :: xsize,ysize
2564	real :: ustin,wsvsurf,timerad
2565	character(80) :: chmess
2566
2567	if(.not.(llesread)) return
2568
2569	open (ilesfile,file='prof.inp.001',status='old',iostat=ierr)
2570	if (ierr /= 0) stop 'ERROR:Prof.inp does not exist'
2571	read (ilesfile,*) kmax
2572	do k=1,kmax
2573	read (ilesfile,*) height(k),thlprof(k),qtprof (k),
2574	. uprof (k),vprof (k),e12prof(k)
2575	enddo
2576	close(ilesfile)
2577
2578	open(ilesfile,file='lscale.inp.001',status='old',iostat=ierr)
2579	if (ierr /= 0) stop 'ERROR:Lscale.inp does not exist'
2580	read (ilesfile,*) kmax2
2581	if (kmax .ne. kmax2) then
2582	print *, 'fichiers prof.inp et lscale.inp incompatibles :'
2583	print *, 'nbre de niveaux : ',kmax,' et ',kmax2
2584	stop 'lecture profiles'
2585	endif
2586	do k=1,kmax
2587	read (ilesfile,*) height(k),ugprof(k),vgprof(k),wfls(k),
2588	. dqtdxls(k),dqtdyls(k),dqtdtls(k),thlpcar(k)
2589	end do
2590	close(ilesfile)
2591
2592	return
2593	end
2594
2595
2596	!======================================================================
2597	subroutine readprofile_armcu(nlev_max,kmax,height,pprof,uprof,
2598	. vprof,thetaprof,tprof,qvprof,rvprof,aprof,bprof)
2599	!======================================================================
2600	implicit none
2601
2602	integer nlev_max,kmax,kmax2
2603	logical :: llesread = .true.
2604
2605	real height(nlev_max),pprof(nlev_max),tprof(nlev_max),
2606	. thetaprof(nlev_max),rvprof(nlev_max),
2607	. qvprof(nlev_max),uprof(nlev_max),vprof(nlev_max),
2608	. aprof(nlev_max+1),bprof(nlev_max+1)
2609
2610	integer, parameter :: ilesfile=1
2611	integer, parameter :: ifile=2
2612	integer :: ierr,irad,imax,jtot,k
2613	logical :: lmoist,lcoriol,ltimedep
2614	real :: xsize,ysize
2615	real :: ustin,wsvsurf,timerad
2616	character(80) :: chmess
2617
2618	if(.not.(llesread)) return
2619
2620	! Read profiles at full levels
2621	IF(nlev_max.EQ.19) THEN
2622	open (ilesfile,file='prof.inp.19',status='old',iostat=ierr)
2623	print *,'On ouvre prof.inp.19'
2624	ELSE
2625	open (ilesfile,file='prof.inp.40',status='old',iostat=ierr)
2626	print *,'On ouvre prof.inp.40'
2627	ENDIF
2628	if (ierr /= 0) stop 'ERROR:Prof.inp does not exist'
2629	read (ilesfile,*) kmax
2630	do k=1,kmax
2631	read (ilesfile,*) height(k) ,pprof(k), uprof(k), vprof(k),
2632	. thetaprof(k) ,tprof(k), qvprof(k),rvprof(k)
2633	enddo
2634	close(ilesfile)
2635
2636	! Vertical coordinates half levels for eta-coordinates (plev = alpha + beta * psurf)
2637	IF(nlev_max.EQ.19) THEN
2638	open (ifile,file='proh.inp.19',status='old',iostat=ierr)
2639	print *,'On ouvre proh.inp.19'
2640	if (ierr /= 0) stop 'ERROR:Proh.inp.19 does not exist'
2641	ELSE
2642	open (ifile,file='proh.inp.40',status='old',iostat=ierr)
2643	print *,'On ouvre proh.inp.40'
2644	if (ierr /= 0) stop 'ERROR:Proh.inp.40 does not exist'
2645	ENDIF
2646	read (ifile,*) kmax
2647	do k=1,kmax
2648	read (ifile,*) jtot,aprof(k),bprof(k)
2649	enddo
2650	close(ifile)
2651
2652	return
2653	end
2654

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

Download in other formats:

Original Format