Context Navigation

1DUTILS.h_with_writelim @ 1702

Last change on this file since 1702 was 1702, checked in by Laurent Fairhead, 12 years ago
To avoid conflict (incompatible and redundant open statements) with conf_gcm FH
File size: 80.5 KB

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

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

Download in other formats:

Original Format