Context Navigation

1DUTILS.h_with_writelim @ 3034

Last change on this file since 3034 was 1763, checked in by Laurent Fairhead, 11 years ago

In 1DUTILS.h disvert renamed disvert0.
There are now 2 routines to get a vertical grid: disvert0 or disvert (../dyn3d/disvert.F90).
The new flag ok_old_disvert allows the use of either disvert0 (ok_old_disvert=y) or disvert (ok_old_disvert=n)

In lmdz1d.F: rlat_rad and rlon_rad for coordinates in radians, cfdt ans sfdt concern the geostrophic wind

File size: 80.7 KB

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

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

Download in other formats:

Original Format