Context Navigation

1DUTILS.h_with_writelim_old @ 1679

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

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

Download in other formats:

Original Format