Context Navigation

1DUTILS.h_with_writelim @ 1608

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

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

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

File size: 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