source: LMDZ5/branches/testing/libf/phylmd/1DUTILS.h @ 2056

Last change on this file since 2056 was 2056, checked in by Laurent Fairhead, 10 years ago

Merged trunk changes r1997:2055 into testing branch

File size: 111.6 KB
Line 
1#include "../dyn3d/conf_gcm.F"
2#include "../dyn3d_common/q_sat.F"
3
4!
5! $Id: conf_unicol.F 1279 2010-08-04 17:20:56Z lahellec $
6!
7!
8!
9      SUBROUTINE conf_unicol
10!
11#ifdef CPP_IOIPSL
12      use IOIPSL
13#else
14! if not using IOIPSL, we still need to use (a local version of) getin
15      use ioipsl_getincom
16#endif
17      IMPLICIT NONE
18!-----------------------------------------------------------------------
19!     Auteurs :   A. Lahellec  .
20!
21!   Declarations :
22!   --------------
23
24#include "compar1d.h"
25#include "flux_arp.h"
26#include "tsoilnudge.h"
27#include "fcg_gcssold.h"
28#include "fcg_racmo.h"
29#include "iniprint.h"
30!
31!
32!   local:
33!   ------
34
35!      CHARACTER ch1*72,ch2*72,ch3*72,ch4*12
36     
37!
38!  -------------------------------------------------------------------
39!
40!      .........    Initilisation parametres du lmdz1D      ..........
41!
42!---------------------------------------------------------------------
43!   initialisations:
44!   ----------------
45
46!Config  Key  = lunout
47!Config  Desc = unite de fichier pour les impressions
48!Config  Def  = 6
49!Config  Help = unite de fichier pour les impressions
50!Config         (defaut sortie standard = 6)
51      lunout=6
52!      CALL getin('lunout', lunout)
53      IF (lunout /= 5 .and. lunout /= 6) THEN
54        OPEN(lunout,FILE='lmdz.out')
55      ENDIF
56
57!Config  Key  = prt_level
58!Config  Desc = niveau d'impressions de débogage
59!Config  Def  = 0
60!Config  Help = Niveau d'impression pour le débogage
61!Config         (0 = minimum d'impression)
62!      prt_level = 0
63!      CALL getin('prt_level',prt_level)
64
65!-----------------------------------------------------------------------
66!  Parametres de controle du run:
67!-----------------------------------------------------------------------
68
69!Config  Key  = restart
70!Config  Desc = on repart des startphy et start1dyn
71!Config  Def  = false
72!Config  Help = les fichiers restart doivent etre renomme en start
73       restart =.false.
74       CALL getin('restart',restart)
75
76!Config  Key  = forcing_type
77!Config  Desc = defines the way the SCM is forced:
78!Config  Def  = 0
79!!Config  Help = 0 ==> forcing_les = .true.
80!             initial profiles from file prof.inp.001
81!             no forcing by LS convergence ; 
82!             surface temperature imposed ;
83!             radiative cooling may be imposed (iflag_radia=0 in physiq.def)
84!         = 1 ==> forcing_radconv = .true.
85!             idem forcing_type = 0, but the imposed radiative cooling
86!             is set to 0 (hence, if iflag_radia=0 in physiq.def, 
87!             then there is no radiative cooling at all)
88!         = 2 ==> forcing_toga = .true.
89!             initial profiles from TOGA-COARE IFA files
90!             LS convergence and SST imposed from TOGA-COARE IFA files
91!         = 3 ==> forcing_GCM2SCM = .true.
92!             initial profiles from the GCM output
93!             LS convergence imposed from the GCM output
94!         = 4 ==> forcing_twpi = .true.
95!             initial profiles from TWPICE nc files
96!             LS convergence and SST imposed from TWPICE nc files
97!         = 5 ==> forcing_rico = .true.
98!             initial profiles from RICO idealized
99!             LS convergence imposed from  RICO (cst) 
100!         = 6 ==> forcing_amma = .true.
101!         = 40 ==> forcing_GCSSold = .true.
102!             initial profile from GCSS file
103!             LS convergence imposed from GCSS file
104!         = 50 ==> forcing_fire = .true.
105!         = 59 ==> forcing_sandu = .true.
106!             initial profiles from sanduref file: see prof.inp.001
107!             SST varying with time and divergence constante: see ifa_sanduref.txt file
108!             Radiation has to be computed interactively
109!         = 60 ==> forcing_astex = .true.
110!             initial profiles from file: see prof.inp.001 
111!             SST,divergence,ug,vg,ufa,vfa varying with time : see ifa_astex.txt file
112!             Radiation has to be computed interactively
113!         = 61 ==> forcing_armcu = .true.
114!             initial profiles from file: see prof.inp.001 
115!             sensible and latent heat flux imposed: see ifa_arm_cu_1.txt
116!             large scale advective forcing & radiative tendencies applied below 1000m: see ifa_arm_cu_2.txt
117!             use geostrophic wind ug=10m/s vg=0m/s. Duration of the case 53100s
118!             Radiation to be switched off
119!
120       forcing_type = 0
121       CALL getin('forcing_type',forcing_type)
122         imp_fcg_gcssold   = .false.
123         ts_fcg_gcssold    = .false. 
124         Tp_fcg_gcssold    = .false. 
125         Tp_ini_gcssold    = .false. 
126         xTurb_fcg_gcssold = .false. 
127        IF (forcing_type .eq.40) THEN
128          CALL getin('imp_fcg',imp_fcg_gcssold)
129          CALL getin('ts_fcg',ts_fcg_gcssold)
130          CALL getin('tp_fcg',Tp_fcg_gcssold)
131          CALL getin('tp_ini',Tp_ini_gcssold)
132          CALL getin('turb_fcg',xTurb_fcg_gcssold)
133        ENDIF
134
135!Config  Key  = ok_flux_surf
136!Config  Desc = forcage ou non par les flux de surface
137!Config  Def  = false
138!Config  Help = forcage ou non par les flux de surface
139       ok_flux_surf =.false.
140       CALL getin('ok_flux_surf',ok_flux_surf)
141
142!Config  Key  = ok_old_disvert
143!Config  Desc = utilisation de l ancien programme disvert0 (dans 1DUTILS.h)
144!Config  Def  = false
145!Config  Help = utilisation de l ancien programme disvert0 (dans 1DUTILS.h)
146       ok_old_disvert = .FALSE.
147       CALL getin('ok_old_disvert',ok_old_disvert)
148
149!Config  Key  = time_ini
150!Config  Desc = meaningless in this  case
151!Config  Def  = 0.
152!Config  Help = 
153       tsurf = 0.
154       CALL getin('time_ini',time_ini)
155
156!Config  Key  = rlat et rlon
157!Config  Desc = latitude et longitude
158!Config  Def  = 0.0  0.0
159!Config  Help = fixe la position de la colonne
160       xlat = 0.
161       xlon = 0.
162       CALL getin('rlat',xlat)
163       CALL getin('rlon',xlon)
164
165!Config  Key  = airephy
166!Config  Desc = Grid cell area
167!Config  Def  = 1.e11
168!Config  Help = 
169       airefi = 1.e11
170       CALL getin('airephy',airefi)
171
172!Config  Key  = nat_surf
173!Config  Desc = surface type
174!Config  Def  = 0 (ocean)
175!Config  Help = 0=ocean,1=land,2=glacier,3=banquise
176       nat_surf = 0.
177       CALL getin('nat_surf',nat_surf)
178
179!Config  Key  = tsurf
180!Config  Desc = surface temperature
181!Config  Def  = 290.
182!Config  Help = not used if type_ts_forcing=1 in lmdz1d.F
183       tsurf = 290.
184       CALL getin('tsurf',tsurf)
185
186!Config  Key  = psurf
187!Config  Desc = surface pressure
188!Config  Def  = 102400.
189!Config  Help = 
190       psurf = 102400.
191       CALL getin('psurf',psurf)
192
193!Config  Key  = zsurf
194!Config  Desc = surface altitude
195!Config  Def  = 0.
196!Config  Help = 
197       zsurf = 0.
198       CALL getin('zsurf',zsurf)
199
200!Config  Key  = rugos
201!Config  Desc = coefficient de frottement
202!Config  Def  = 0.0001
203!Config  Help = calcul du Cdrag
204       rugos = 0.0001
205       CALL getin('rugos',rugos)
206
207!Config  Key  = wtsurf et wqsurf
208!Config  Desc = ???
209!Config  Def  = 0.0 0.0
210!Config  Help = 
211       wtsurf = 0.0
212       wqsurf = 0.0
213       CALL getin('wtsurf',wtsurf)
214       CALL getin('wqsurf',wqsurf)
215
216!Config  Key  = albedo
217!Config  Desc = albedo
218!Config  Def  = 0.09
219!Config  Help = 
220       albedo = 0.09
221       CALL getin('albedo',albedo)
222
223!Config  Key  = agesno
224!Config  Desc = age de la neige
225!Config  Def  = 30.0
226!Config  Help = 
227       xagesno = 30.0
228       CALL getin('agesno',xagesno)
229
230!Config  Key  = restart_runoff
231!Config  Desc = age de la neige
232!Config  Def  = 30.0
233!Config  Help = 
234       restart_runoff = 0.0
235       CALL getin('restart_runoff',restart_runoff)
236
237!Config  Key  = qsolinp
238!Config  Desc = initial bucket water content (kg/m2) when land (5std)
239!Config  Def  = 30.0
240!Config  Help = 
241       qsolinp = 1.
242       CALL getin('qsolinp',qsolinp)
243
244!Config  Key  = zpicinp
245!Config  Desc = denivellation orographie
246!Config  Def  = 300.
247!Config  Help =  input brise
248       zpicinp = 300.
249       CALL getin('zpicinp',zpicinp)
250!Config key = nudge_tsoil
251!Config  Desc = activation of soil temperature nudging
252!Config  Def  = .FALSE.
253!Config  Help = ...
254
255       nudge_tsoil=.FALSE.
256       CALL getin('nudge_tsoil',nudge_tsoil)
257
258!Config key = isoil_nudge
259!Config  Desc = level number where soil temperature is nudged
260!Config  Def  = 3
261!Config  Help = ...
262
263       isoil_nudge=3
264       CALL getin('isoil_nudge',isoil_nudge)
265
266!Config key = Tsoil_nudge
267!Config  Desc = target temperature for tsoil(isoil_nudge)
268!Config  Def  = 300.
269!Config  Help = ...
270
271       Tsoil_nudge=300.
272       CALL getin('Tsoil_nudge',Tsoil_nudge)
273
274!Config key = tau_soil_nudge
275!Config  Desc = nudging relaxation time for tsoil
276!Config  Def  = 3600.
277!Config  Help = ...
278
279       tau_soil_nudge=3600.
280       CALL getin('tau_soil_nudge',tau_soil_nudge)
281
282
283
284
285      write(lunout,*)' +++++++++++++++++++++++++++++++++++++++'
286      write(lunout,*)' Configuration des parametres du gcm1D: '
287      write(lunout,*)' +++++++++++++++++++++++++++++++++++++++'
288      write(lunout,*)' restart = ', restart
289      write(lunout,*)' forcing_type = ', forcing_type
290      write(lunout,*)' time_ini = ', time_ini
291      write(lunout,*)' rlat = ', xlat
292      write(lunout,*)' rlon = ', xlon
293      write(lunout,*)' airephy = ', airefi
294      write(lunout,*)' nat_surf = ', nat_surf
295      write(lunout,*)' tsurf = ', tsurf
296      write(lunout,*)' psurf = ', psurf
297      write(lunout,*)' zsurf = ', zsurf
298      write(lunout,*)' rugos = ', rugos
299      write(lunout,*)' wtsurf = ', wtsurf
300      write(lunout,*)' wqsurf = ', wqsurf
301      write(lunout,*)' albedo = ', albedo
302      write(lunout,*)' xagesno = ', xagesno
303      write(lunout,*)' restart_runoff = ', restart_runoff
304      write(lunout,*)' qsolinp = ', qsolinp
305      write(lunout,*)' zpicinp = ', zpicinp
306      write(lunout,*)' nudge_tsoil = ', nudge_tsoil
307      write(lunout,*)' isoil_nudge = ', isoil_nudge
308      write(lunout,*)' Tsoil_nudge = ', Tsoil_nudge
309      write(lunout,*)' tau_soil_nudge = ', tau_soil_nudge
310      IF (forcing_type .eq.40) THEN
311        write(lunout,*) '--- Forcing type GCSS Old --- with:'
312        write(lunout,*)'imp_fcg',imp_fcg_gcssold
313        write(lunout,*)'ts_fcg',ts_fcg_gcssold
314        write(lunout,*)'tp_fcg',Tp_fcg_gcssold
315        write(lunout,*)'tp_ini',Tp_ini_gcssold
316        write(lunout,*)'xturb_fcg',xTurb_fcg_gcssold
317      ENDIF
318
319      write(lunout,*)' +++++++++++++++++++++++++++++++++++++++'
320      write(lunout,*)
321!
322      RETURN
323      END
324!
325! $Id: dyn1deta0.F 1279 2010/07/30 A Lahellec$
326!
327!
328      SUBROUTINE dyn1deta0(fichnom,plev,play,phi,phis,presnivs,                 &
329     &                          ucov,vcov,temp,q,omega2)
330      USE dimphy
331      USE mod_grid_phy_lmdz
332      USE mod_phys_lmdz_para
333      USE iophy
334      USE phys_state_var_mod
335      USE iostart
336      USE write_field_phy
337      USE infotrac
338      use control_mod
339
340      IMPLICIT NONE
341!=======================================================
342! Ecriture du fichier de redemarrage sous format NetCDF
343!=======================================================
344!   Declarations:
345!   -------------
346#include "dimensions.h"
347#include "comconst.h"
348#include "temps.h"
349!!#include "control.h"
350#include "logic.h"
351#include "netcdf.inc"
352
353!   Arguments:
354!   ----------
355      CHARACTER*(*) fichnom
356!Al1 plev tronque pour .nc mais plev(klev+1):=0
357      real :: plev(klon,klev),play (klon,klev),phi(klon,klev)
358      real :: presnivs(klon,klev)
359      real :: ucov(klon,klev),vcov(klon,klev),temp(klon,klev)
360      real :: q(klon,klev,nqtot),omega2(klon,klev)
361!      real :: ug(klev),vg(klev),fcoriolis
362      real :: phis(klon) 
363
364!   Variables locales pour NetCDF:
365!   ------------------------------
366      INTEGER iq
367      INTEGER length
368      PARAMETER (length = 100)
369      REAL tab_cntrl(length) ! tableau des parametres du run
370      character*4 nmq(nqtot)
371      character*12 modname
372      character*80 abort_message
373      LOGICAL found
374
375      modname = 'dyn1deta0 : '
376      nmq(1)="vap"
377      nmq(2)="cond"
378      do iq=3,nqtot
379        write(nmq(iq),'("tra",i1)') iq-2
380      enddo
381      print*,'in dyn1deta0 ',fichnom,klon,klev,nqtot
382      CALL open_startphy(fichnom)
383      print*,'after open startphy ',fichnom,nmq
384
385!
386! Lecture des parametres de controle:
387!
388      CALL get_var("controle",tab_cntrl)
389       
390
391      im         = tab_cntrl(1)
392      jm         = tab_cntrl(2)
393      lllm       = tab_cntrl(3)
394      day_ref    = tab_cntrl(4)
395      annee_ref  = tab_cntrl(5)
396!      rad        = tab_cntrl(6)
397!      omeg       = tab_cntrl(7)
398!      g          = tab_cntrl(8)
399!      cpp        = tab_cntrl(9)
400!      kappa      = tab_cntrl(10)
401!      daysec     = tab_cntrl(11)
402!      dtvr       = tab_cntrl(12)
403!      etot0      = tab_cntrl(13)
404!      ptot0      = tab_cntrl(14)
405!      ztot0      = tab_cntrl(15)
406!      stot0      = tab_cntrl(16)
407!      ang0       = tab_cntrl(17)
408!      pa         = tab_cntrl(18)
409!      preff      = tab_cntrl(19)
410!
411!      clon       = tab_cntrl(20)
412!      clat       = tab_cntrl(21)
413!      grossismx  = tab_cntrl(22)
414!      grossismy  = tab_cntrl(23)
415!
416      IF ( tab_cntrl(24).EQ.1. )  THEN
417        fxyhypb  =.true.
418!        dzoomx   = tab_cntrl(25)
419!        dzoomy   = tab_cntrl(26)
420!        taux     = tab_cntrl(28)
421!        tauy     = tab_cntrl(29)
422      ELSE
423        fxyhypb = .false.
424        ysinus  = .false.
425        IF( tab_cntrl(27).EQ.1. ) ysinus =.true. 
426      ENDIF
427
428      day_ini = tab_cntrl(30)
429      itau_dyn = tab_cntrl(31)
430!   .................................................................
431!
432!
433!      PRINT*,'rad,omeg,g,cpp,kappa',rad,omeg,g,cpp,kappa
434!Al1
435       Print*,'day_ref,annee_ref,day_ini,itau_dyn',                         &
436     &              day_ref,annee_ref,day_ini,itau_dyn
437
438!  Lecture des champs
439!
440      plev(1,klev+1)=0.
441      CALL get_field("plev",plev,found)
442      IF (.NOT. found) PRINT*, modname//'Le champ <Plev> est absent'
443      CALL get_field("play",play,found)
444      IF (.NOT. found) PRINT*, modname//'Le champ <Play> est absent'
445      CALL get_field("phi",phi,found)
446      IF (.NOT. found) PRINT*, modname//'Le champ <Phi> est absent'
447      CALL get_field("phis",phis,found)
448      IF (.NOT. found) PRINT*, modname//'Le champ <Phis> est absent'
449      CALL get_field("presnivs",presnivs,found)
450      IF (.NOT. found) PRINT*, modname//'Le champ <Presnivs> est absent'
451      CALL get_field("ucov",ucov,found)
452      IF (.NOT. found) PRINT*, modname//'Le champ <ucov> est absent'
453      CALL get_field("vcov",vcov,found)
454      IF (.NOT. found) PRINT*, modname//'Le champ <vcov> est absent'
455      CALL get_field("temp",temp,found)
456      IF (.NOT. found) PRINT*, modname//'Le champ <temp> est absent'
457      CALL get_field("omega2",omega2,found)
458      IF (.NOT. found) PRINT*, modname//'Le champ <omega2> est absent'
459
460      Do iq=1,nqtot
461        CALL get_field("q"//nmq(iq),q(:,:,iq),found)
462        IF (.NOT.found)PRINT*, modname//'Le champ <q'//nmq//'> est absent'
463      EndDo
464
465      CALL close_startphy
466      print*,' close startphy',fichnom,play(1,1),play(1,klev),temp(1,klev)
467!
468      RETURN
469      END
470!
471! $Id: dyn1dredem.F 1279 2010/07/29 A Lahellec$
472!
473!
474      SUBROUTINE dyn1dredem(fichnom,plev,play,phi,phis,presnivs,           &
475     &                          ucov,vcov,temp,q,omega2)
476      USE dimphy
477      USE mod_grid_phy_lmdz
478      USE mod_phys_lmdz_para
479      USE phys_state_var_mod
480      USE iostart
481      USE infotrac
482      use control_mod
483
484      IMPLICIT NONE
485!=======================================================
486! Ecriture du fichier de redemarrage sous format NetCDF
487!=======================================================
488!   Declarations:
489!   -------------
490#include "dimensions.h"
491#include "comconst.h"
492#include "temps.h"
493!!#include "control.h"
494#include "logic.h"
495#include "netcdf.inc"
496
497!   Arguments:
498!   ----------
499      CHARACTER*(*) fichnom
500!Al1 plev tronque pour .nc mais plev(klev+1):=0
501      real :: plev(klon,klev),play (klon,klev),phi(klon,klev)
502      real :: presnivs(klon,klev)
503      real :: ucov(klon,klev),vcov(klon,klev),temp(klon,klev)
504      real :: q(klon,klev,nqtot)
505      real :: omega2(klon,klev),rho(klon,klev+1)
506!      real :: ug(klev),vg(klev),fcoriolis
507      real :: phis(klon) 
508
509!   Variables locales pour NetCDF:
510!   ------------------------------
511      INTEGER nid
512      INTEGER ierr
513      INTEGER iq,l
514      INTEGER length
515      PARAMETER (length = 100)
516      REAL tab_cntrl(length) ! tableau des parametres du run
517      character*4 nmq(nqtot)
518      character*20 modname
519      character*80 abort_message
520!
521      INTEGER nb
522      SAVE nb
523      DATA nb / 0 /
524
525      CALL open_restartphy(fichnom)
526      print*,'redm1 ',fichnom,klon,klev,nqtot
527      nmq(1)="vap"
528      nmq(2)="cond"
529      nmq(3)="tra1"
530      nmq(4)="tra2"
531
532      modname = 'dyn1dredem'
533      ierr = NF_OPEN(fichnom, NF_WRITE, nid)
534      IF (ierr .NE. NF_NOERR) THEN
535         abort_message="Pb. d ouverture "//fichnom
536         CALL abort_gcm('Modele 1D',abort_message,1)
537      ENDIF
538
539      DO l=1,length
540       tab_cntrl(l) = 0.
541      ENDDO
542       tab_cntrl(1)  = FLOAT(iim)
543       tab_cntrl(2)  = FLOAT(jjm)
544       tab_cntrl(3)  = FLOAT(llm)
545       tab_cntrl(4)  = FLOAT(day_ref)
546       tab_cntrl(5)  = FLOAT(annee_ref)
547       tab_cntrl(6)  = rad
548       tab_cntrl(7)  = omeg
549       tab_cntrl(8)  = g
550       tab_cntrl(9)  = cpp
551       tab_cntrl(10) = kappa
552       tab_cntrl(11) = daysec
553       tab_cntrl(12) = dtvr
554!       tab_cntrl(13) = etot0
555!       tab_cntrl(14) = ptot0
556!       tab_cntrl(15) = ztot0
557!       tab_cntrl(16) = stot0
558!       tab_cntrl(17) = ang0
559!       tab_cntrl(18) = pa
560!       tab_cntrl(19) = preff
561!
562!    .....    parametres  pour le zoom      ......   
563
564!       tab_cntrl(20)  = clon
565!       tab_cntrl(21)  = clat
566!       tab_cntrl(22)  = grossismx
567!       tab_cntrl(23)  = grossismy
568!
569      IF ( fxyhypb )   THEN
570       tab_cntrl(24) = 1.
571!       tab_cntrl(25) = dzoomx
572!       tab_cntrl(26) = dzoomy
573       tab_cntrl(27) = 0.
574!       tab_cntrl(28) = taux
575!       tab_cntrl(29) = tauy
576      ELSE
577       tab_cntrl(24) = 0.
578!       tab_cntrl(25) = dzoomx
579!       tab_cntrl(26) = dzoomy
580       tab_cntrl(27) = 0.
581       tab_cntrl(28) = 0.
582       tab_cntrl(29) = 0.
583       IF( ysinus )  tab_cntrl(27) = 1.
584      ENDIF
585!Al1 iday_end -> day_end
586       tab_cntrl(30) = FLOAT(day_end)
587       tab_cntrl(31) = FLOAT(itau_dyn + itaufin)
588!
589      CALL put_var("controle","Param. de controle Dyn1D",tab_cntrl)
590!
591
592!  Ecriture/extension de la coordonnee temps
593
594      nb = nb + 1
595
596!  Ecriture des champs
597!
598      CALL put_field("plev","p interfaces sauf la nulle",plev)
599      CALL put_field("play","",play)
600      CALL put_field("phi","geopotentielle",phi)
601      CALL put_field("phis","geopotentiell de surface",phis)
602      CALL put_field("presnivs","",presnivs)
603      CALL put_field("ucov","",ucov)
604      CALL put_field("vcov","",vcov)
605      CALL put_field("temp","",temp)
606      CALL put_field("omega2","",omega2)
607
608      Do iq=1,nqtot
609        CALL put_field("q"//nmq(iq),"eau vap ou condens et traceurs",           &
610     &                                                      q(:,:,iq))
611      EndDo
612      CALL close_restartphy
613
614!
615      RETURN
616      END
617      SUBROUTINE gr_fi_dyn(nfield,ngrid,im,jm,pfi,pdyn)
618      IMPLICIT NONE
619!=======================================================================
620!   passage d'un champ de la grille scalaire a la grille physique
621!=======================================================================
622 
623!-----------------------------------------------------------------------
624!   declarations:
625!   -------------
626 
627      INTEGER im,jm,ngrid,nfield
628      REAL pdyn(im,jm,nfield)
629      REAL pfi(ngrid,nfield)
630 
631      INTEGER i,j,ifield,ig
632 
633!-----------------------------------------------------------------------
634!   calcul:
635!   -------
636 
637      DO ifield=1,nfield
638!   traitement des poles
639         DO i=1,im
640            pdyn(i,1,ifield)=pfi(1,ifield)
641            pdyn(i,jm,ifield)=pfi(ngrid,ifield)
642         ENDDO
643 
644!   traitement des point normaux
645         DO j=2,jm-1
646            ig=2+(j-2)*(im-1)
647            CALL SCOPY(im-1,pfi(ig,ifield),1,pdyn(1,j,ifield),1)
648            pdyn(im,j,ifield)=pdyn(1,j,ifield)
649         ENDDO
650      ENDDO
651 
652      RETURN
653      END
654 
655 
656
657      SUBROUTINE abort_gcm(modname, message, ierr)
658 
659      USE IOIPSL
660!
661! Stops the simulation cleanly, closing files and printing various
662! comments
663!
664!  Input: modname = name of calling program
665!         message = stuff to print
666!         ierr    = severity of situation ( = 0 normal )
667 
668      character(len=*) modname
669      integer ierr
670      character(len=*) message
671 
672      write(*,*) 'in abort_gcm'
673      call histclo
674!     call histclo(2)
675!     call histclo(3)
676!     call histclo(4)
677!     call histclo(5)
678      write(*,*) 'out of histclo'
679      write(*,*) 'Stopping in ', modname
680      write(*,*) 'Reason = ',message
681      call getin_dump
682!
683      if (ierr .eq. 0) then
684        write(*,*) 'Everything is cool'
685      else
686        write(*,*) 'Houston, we have a problem ', ierr
687      endif
688      STOP
689      END
690      REAL FUNCTION fq_sat(kelvin, millibar)
691!
692      IMPLICIT none
693!======================================================================
694! Autheur(s): Z.X. Li (LMD/CNRS)
695! Objet: calculer la vapeur d'eau saturante (formule Centre Euro.)
696!======================================================================
697! Arguments:
698! kelvin---input-R: temperature en Kelvin
699! millibar--input-R: pression en mb
700!
701! fq_sat----output-R: vapeur d'eau saturante en kg/kg
702!======================================================================
703!
704      REAL kelvin, millibar
705!
706      REAL r2es
707      PARAMETER (r2es=611.14 *18.0153/28.9644)
708!
709      REAL r3les, r3ies, r3es
710      PARAMETER (R3LES=17.269)
711      PARAMETER (R3IES=21.875)
712!
713      REAL r4les, r4ies, r4es
714      PARAMETER (R4LES=35.86)
715      PARAMETER (R4IES=7.66)
716!
717      REAL rtt
718      PARAMETER (rtt=273.16)
719!
720      REAL retv
721      PARAMETER (retv=28.9644/18.0153 - 1.0)
722!
723      REAL zqsat
724      REAL temp, pres
725!     ------------------------------------------------------------------
726!
727!
728      temp = kelvin
729      pres = millibar * 100.0
730!      write(*,*)'kelvin,millibar=',kelvin,millibar
731!      write(*,*)'temp,pres=',temp,pres
732!
733      IF (temp .LE. rtt) THEN
734         r3es = r3ies
735         r4es = r4ies
736      ELSE
737         r3es = r3les
738         r4es = r4les
739      ENDIF
740!
741      zqsat=r2es/pres * EXP ( r3es*(temp-rtt) / (temp-r4es) )
742      zqsat=MIN(0.5,ZQSAT)
743      zqsat=zqsat/(1.-retv  *zqsat)
744!
745      fq_sat = zqsat
746!
747      RETURN
748      END
749 
750      SUBROUTINE gr_dyn_fi(nfield,im,jm,ngrid,pdyn,pfi)
751      IMPLICIT NONE
752!=======================================================================
753!   passage d'un champ de la grille scalaire a la grille physique
754!=======================================================================
755 
756!-----------------------------------------------------------------------
757!   declarations:
758!   -------------
759 
760      INTEGER im,jm,ngrid,nfield
761      REAL pdyn(im,jm,nfield)
762      REAL pfi(ngrid,nfield)
763 
764      INTEGER j,ifield,ig
765 
766!-----------------------------------------------------------------------
767!   calcul:
768!   -------
769 
770      IF(ngrid.NE.2+(jm-2)*(im-1).AND.ngrid.NE.1)                          &
771     &    STOP 'probleme de dim'
772!   traitement des poles
773      CALL SCOPY(nfield,pdyn,im*jm,pfi,ngrid)
774      CALL SCOPY(nfield,pdyn(1,jm,1),im*jm,pfi(ngrid,1),ngrid)
775 
776!   traitement des point normaux
777      DO ifield=1,nfield
778         DO j=2,jm-1
779            ig=2+(j-2)*(im-1)
780            CALL SCOPY(im-1,pdyn(1,j,ifield),1,pfi(ig,ifield),1)
781         ENDDO
782      ENDDO
783 
784      RETURN
785      END
786 
787      SUBROUTINE disvert0(pa,preff,ap,bp,dpres,presnivs,nivsigs,nivsig)
788 
789!    Ancienne version disvert dont on a modifie nom pour utiliser
790!    le disvert de dyn3d (qui permet d'utiliser grille avec ab,bp imposes)
791!    (MPL 18092012)
792!
793!    Auteur :  P. Le Van .
794!
795      IMPLICIT NONE
796 
797#include "dimensions.h"
798#include "paramet.h"
799!
800!=======================================================================
801!
802!
803!    s = sigma ** kappa   :  coordonnee  verticale
804!    dsig(l)            : epaisseur de la couche l ds la coord.  s
805!    sig(l)             : sigma a l'interface des couches l et l-1
806!    ds(l)              : distance entre les couches l et l-1 en coord.s
807!
808!=======================================================================
809!
810      REAL pa,preff
811      REAL ap(llmp1),bp(llmp1),dpres(llm),nivsigs(llm),nivsig(llmp1)
812      REAL presnivs(llm)
813!
814!   declarations:
815!   -------------
816!
817      REAL sig(llm+1),dsig(llm)
818!
819      INTEGER l
820      REAL snorm
821      REAL alpha,beta,gama,delta,deltaz,h
822      INTEGER np,ierr
823      REAL pi,x
824      REAL :: scaleheight=8.
825 
826!-----------------------------------------------------------------------
827!
828      pi=2.*ASIN(1.)
829 
830      OPEN(99,file='sigma.def',status='old',form='formatted',                   &
831     &   iostat=ierr)
832 
833!-----------------------------------------------------------------------
834!   cas 1 on lit les options dans sigma.def:
835!   ----------------------------------------
836 
837      IF (ierr.eq.0) THEN
838 
839      print*,'WARNING!!! on lit les options dans sigma.def'
840      READ(99,*) deltaz
841      READ(99,*) h
842      READ(99,*) beta
843      READ(99,*) gama
844      READ(99,*) delta
845      READ(99,*) np
846      CLOSE(99)
847      alpha=deltaz/(llm*h)
848!
849 
850       DO 1  l = 1, llm
851       dsig(l) = (alpha+(1.-alpha)*exp(-beta*(llm-l)))*                    &
852     &          ( (tanh(gama*l)/tanh(gama*llm))**np +                      &
853     &            (1.-l/FLOAT(llm))*delta )
854   1   CONTINUE
855 
856       sig(1)=1.
857       DO 101 l=1,llm-1
858          sig(l+1)=sig(l)*(1.-dsig(l))/(1.+dsig(l))
859101    CONTINUE
860       sig(llm+1)=0.
861 
862       DO 2  l = 1, llm
863       dsig(l) = sig(l)-sig(l+1)
864   2   CONTINUE
865!
866 
867      ELSE
868!-----------------------------------------------------------------------
869!   cas 2 ancienne discretisation (LMD5...):
870!   ----------------------------------------
871 
872      PRINT*,'WARNING!!! Ancienne discretisation verticale'
873 
874      h=7.
875      snorm  = 0.
876      DO l = 1, llm
877         x = 2.*asin(1.) * (FLOAT(l)-0.5) / float(llm+1)
878         dsig(l) = 1.0 + 7.0 * SIN(x)**2
879         snorm = snorm + dsig(l)
880      ENDDO
881      snorm = 1./snorm
882      DO l = 1, llm
883         dsig(l) = dsig(l)*snorm
884      ENDDO
885      sig(llm+1) = 0.
886      DO l = llm, 1, -1
887         sig(l) = sig(l+1) + dsig(l)
888      ENDDO
889 
890      ENDIF
891 
892 
893      DO l=1,llm
894        nivsigs(l) = FLOAT(l)
895      ENDDO
896 
897      DO l=1,llmp1
898        nivsig(l)= FLOAT(l)
899      ENDDO
900 
901!
902!    ....  Calculs  de ap(l) et de bp(l)  ....
903!    .........................................
904!
905!
906!   .....  pa et preff sont lus  sur les fichiers start par lectba  .....
907!
908 
909      bp(llmp1) =   0.
910 
911      DO l = 1, llm
912!c
913!cc    ap(l) = 0.
914!cc    bp(l) = sig(l)
915 
916      bp(l) = EXP( 1. -1./( sig(l)*sig(l)) )
917      ap(l) = pa * ( sig(l) - bp(l) )
918!
919      ENDDO
920      ap(llmp1) = pa * ( sig(llmp1) - bp(llmp1) )
921 
922      PRINT *,' BP '
923      PRINT *,  bp
924      PRINT *,' AP '
925      PRINT *,  ap
926 
927   print*,'scaleheight=',scaleheight
928  DO l = 1, llm
929     dpres(l) = bp(l) - bp(l+1)
930     presnivs(l) = 0.5 *( ap(l)+bp(l)*preff + ap(l+1)+bp(l+1)*preff )
931     write(*, *)'PRESNIVS(', l, ')=', presnivs(l), ' Z ~ ', &
932          log(preff/presnivs(l))*scaleheight &
933          , ' DZ ~ ', scaleheight*log((ap(l)+bp(l)*preff)/ &
934          max(ap(l+1)+bp(l+1)*preff, 1.e-10))
935  ENDDO
936
937 
938      PRINT *,' PRESNIVS '
939      PRINT *,presnivs
940 
941      RETURN
942      END
943
944!======================================================================
945       SUBROUTINE read_tsurf1d(knon,sst_out)
946
947! This subroutine specifies the surface temperature to be used in 1D simulations
948
949      USE dimphy, ONLY : klon
950
951      INTEGER, INTENT(IN)                  :: knon     ! nomber of points on compressed grid
952      REAL, DIMENSION(klon), INTENT(OUT)   :: sst_out  ! tsurf used to force the single-column model
953
954       INTEGER :: i
955! COMMON defined in lmdz1d.F:
956       real ts_cur
957       common /sst_forcing/ts_cur
958
959       DO i = 1, knon
960        sst_out(i) = ts_cur
961       ENDDO
962
963      END SUBROUTINE read_tsurf1d
964
965!===============================================================
966      subroutine advect_vert(llm,w,dt,q,plev)
967!===============================================================
968!   Schema amont pour l'advection verticale en 1D
969!   w est la vitesse verticale dp/dt en Pa/s
970!   Traitement en volumes finis
971!   d / dt ( zm q ) = delta_z ( omega q )
972!   d / dt ( zm ) = delta_z ( omega )
973!   avec zm = delta_z ( p )
974!   si * designe la valeur au pas de temps t+dt
975!   zm*(l) q*(l) - zm(l) q(l) = w(l+1) q(l+1) - w(l) q(l)
976!   zm*(l) -zm(l) = w(l+1) - w(l)
977!   avec w=omega * dt
978!---------------------------------------------------------------
979      implicit none
980! arguments
981      integer llm
982      real w(llm+1),q(llm),plev(llm+1),dt
983
984! local
985      integer l
986      real zwq(llm+1),zm(llm+1),zw(llm+1)
987      real qold
988
989!---------------------------------------------------------------
990
991      do l=1,llm
992         zw(l)=dt*w(l)
993         zm(l)=plev(l)-plev(l+1)
994         zwq(l)=q(l)*zw(l)
995      enddo
996      zwq(llm+1)=0.
997      zw(llm+1)=0.
998 
999      do l=1,llm
1000         qold=q(l)
1001         q(l)=(q(l)*zm(l)+zwq(l+1)-zwq(l))/(zm(l)+zw(l+1)-zw(l))
1002         print*,'ADV Q ',zm(l),zw(l),zwq(l),qold,q(l)
1003      enddo
1004
1005 
1006      return
1007      end
1008
1009!===============================================================
1010
1011
1012       SUBROUTINE advect_va(llm,omega,d_t_va,d_q_va,d_u_va,d_v_va,              &
1013     &                q,temp,u,v,play)
1014!itlmd
1015!----------------------------------------------------------------------
1016!   Calcul de l'advection verticale (ascendance et subsidence) de
1017!   température et d'humidité. Hypothèse : ce qui rentre de l'extérieur
1018!   a les mêmes caractéristiques que l'air de la colonne 1D (WTG) ou
1019!   sans WTG rajouter une advection horizontale
1020!---------------------------------------------------------------------- 
1021        implicit none
1022#include "YOMCST.h"
1023!        argument
1024        integer llm
1025        real  omega(llm+1),d_t_va(llm), d_q_va(llm,3)
1026        real  d_u_va(llm), d_v_va(llm)
1027        real  q(llm,3),temp(llm)
1028        real  u(llm),v(llm)
1029        real  play(llm)
1030! interne
1031        integer l
1032        real alpha,omgdown,omgup
1033
1034      do l= 1,llm
1035       if(l.eq.1) then
1036!si omgup pour la couche 1, alors tendance nulle
1037        omgdown=max(omega(2),0.0)
1038        alpha = rkappa*temp(l)*(1.+q(l,1)*rv/rd)/(play(l)*(1.+q(l,1)))
1039        d_t_va(l)= alpha*(omgdown)-omgdown*(temp(l)-temp(l+1))             &
1040     &       /(play(l)-play(l+1))
1041
1042        d_q_va(l,:)= -omgdown*(q(l,:)-q(l+1,:))/(play(l)-play(l+1))             
1043
1044        d_u_va(l)= -omgdown*(u(l)-u(l+1))/(play(l)-play(l+1))             
1045        d_v_va(l)= -omgdown*(v(l)-v(l+1))/(play(l)-play(l+1))             
1046
1047       
1048       elseif(l.eq.llm) then
1049        omgup=min(omega(l),0.0)
1050        alpha = rkappa*temp(l)*(1.+q(l,1)*rv/rd)/(play(l)*(1.+q(l,1)))
1051        d_t_va(l)= alpha*(omgup)-                                          &
1052
1053!bug?     &              omgup*(temp(l-1)-temp(l))/(play(l-1)-plev(l))
1054     &              omgup*(temp(l-1)-temp(l))/(play(l-1)-play(l))
1055        d_q_va(l,:)= -omgup*(q(l-1,:)-q(l,:))/(play(l-1)-play(l))
1056        d_u_va(l)= -omgup*(u(l-1)-u(l))/(play(l-1)-play(l))
1057        d_v_va(l)= -omgup*(v(l-1)-v(l))/(play(l-1)-play(l))
1058       
1059       else
1060        omgup=min(omega(l),0.0)
1061        omgdown=max(omega(l+1),0.0)
1062        alpha = rkappa*temp(l)*(1.+q(l,1)*rv/rd)/(play(l)*(1.+q(l,1)))
1063        d_t_va(l)= alpha*(omgup+omgdown)-omgdown*(temp(l)-temp(l+1))       &
1064     &              /(play(l)-play(l+1))-                                  &
1065!bug?     &              omgup*(temp(l-1)-temp(l))/(play(l-1)-plev(l))
1066     &              omgup*(temp(l-1)-temp(l))/(play(l-1)-play(l))
1067!      print*, '  ??? '
1068
1069        d_q_va(l,:)= -omgdown*(q(l,:)-q(l+1,:))                            &
1070     &              /(play(l)-play(l+1))-                                  &
1071     &              omgup*(q(l-1,:)-q(l,:))/(play(l-1)-play(l)) 
1072        d_u_va(l)= -omgdown*(u(l)-u(l+1))                                  &
1073     &              /(play(l)-play(l+1))-                                  &
1074     &              omgup*(u(l-1)-u(l))/(play(l-1)-play(l)) 
1075        d_v_va(l)= -omgdown*(v(l)-v(l+1))                                  &
1076     &              /(play(l)-play(l+1))-                                  &
1077     &              omgup*(v(l-1)-v(l))/(play(l-1)-play(l))
1078       
1079      endif
1080         
1081      enddo
1082!fin itlmd
1083        return
1084        end
1085!       SUBROUTINE lstendH(llm,omega,d_t_va,d_q_va,d_u_va,d_v_va,
1086       SUBROUTINE lstendH(llm,nqtot,omega,d_t_va,d_q_va,                        &
1087     &                q,temp,u,v,play)
1088!itlmd
1089!----------------------------------------------------------------------
1090!   Calcul de l'advection verticale (ascendance et subsidence) de
1091!   température et d'humidité. Hypothèse : ce qui rentre de l'extérieur
1092!   a les mêmes caractéristiques que l'air de la colonne 1D (WTG) ou
1093!   sans WTG rajouter une advection horizontale
1094!---------------------------------------------------------------------- 
1095        implicit none
1096#include "YOMCST.h"
1097!        argument
1098        integer llm,nqtot
1099        real  omega(llm+1),d_t_va(llm), d_q_va(llm,nqtot)
1100!        real  d_u_va(llm), d_v_va(llm)
1101        real  q(llm,nqtot),temp(llm)
1102        real  u(llm),v(llm)
1103        real  play(llm)
1104        real cor(llm)
1105!        real dph(llm),dudp(llm),dvdp(llm),dqdp(llm),dtdp(llm)
1106        real dph(llm),dqdp(llm),dtdp(llm)
1107! interne
1108        integer k
1109        real omdn,omup
1110
1111!        dudp=0.
1112!        dvdp=0.
1113        dqdp=0.
1114        dtdp=0.
1115!        d_u_va=0.
1116!        d_v_va=0.
1117
1118      cor(:) = rkappa*temp*(1.+q(:,1)*rv/rd)/(play*(1.+q(:,1)))
1119
1120
1121      do k=2,llm-1
1122
1123       dph  (k-1) = (play()- play(k-1  ))
1124!       dudp (k-1) = (u   ()- u   (k-1  ))/dph(k-1)
1125!       dvdp (k-1) = (v   ()- v   (k-1  ))/dph(k-1)
1126       dqdp (k-1) = (q   (k,1)- q   (k-1,1))/dph(k-1)
1127       dtdp (k-1) = (temp()- temp(k-1  ))/dph(k-1)
1128
1129      enddo
1130
1131!      dudp (  llm  ) = dudp ( llm-1 )
1132!      dvdp (  llm  ) = dvdp ( llm-1 )
1133      dqdp (  llm  ) = dqdp ( llm-1 )
1134      dtdp (  llm  ) = dtdp ( llm-1 )
1135
1136      do k=2,llm-1
1137      omdn=max(0.0,omega(k+1))
1138      omup=min(0.0,omega( k ))
1139
1140!      d_u_va(k)  = -omdn*dudp(k)-omup*dudp(k-1)
1141!      d_v_va(k)  = -omdn*dvdp(k)-omup*dvdp(k-1)
1142      d_q_va(k,1)= -omdn*dqdp(k)-omup*dqdp(k-1)
1143      d_t_va(k)  = -omdn*dtdp(k)-omup*dtdp(k-1)+(omup+omdn)*cor(k)
1144      enddo
1145
1146      omdn=max(0.0,omega( 2 ))
1147      omup=min(0.0,omega(llm))
1148!      d_u_va( 1 )   = -omdn*dudp( 1 )
1149!      d_u_va(llm)   = -omup*dudp(llm)
1150!      d_v_va( 1 )   = -omdn*dvdp( 1 )
1151!      d_v_va(llm)   = -omup*dvdp(llm)
1152      d_q_va( 1 ,1) = -omdn*dqdp( 1 )
1153      d_q_va(llm,1) = -omup*dqdp(llm)
1154      d_t_va( 1 )   = -omdn*dtdp( 1 )+omdn*cor( 1 )
1155      d_t_va(llm)   = -omup*dtdp(llm)!+omup*cor(llm)
1156
1157!      if(abs(rlat(1))>10.) then
1158!     Calculate the tendency due agestrophic motions
1159!      du_age = fcoriolis*(v-vg)
1160!      dv_age = fcoriolis*(ug-u)
1161!      endif
1162
1163!       call writefield_phy('d_t_va',d_t_va,llm)
1164
1165          return
1166         end
1167
1168!======================================================================
1169      SUBROUTINE read_togacoare(fich_toga,nlev_toga,nt_toga                     &
1170     &             ,ts_toga,plev_toga,t_toga,q_toga,u_toga,v_toga,w_toga        &
1171     &             ,ht_toga,vt_toga,hq_toga,vq_toga)
1172      implicit none
1173
1174!-------------------------------------------------------------------------
1175! Read TOGA-COARE forcing data
1176!-------------------------------------------------------------------------
1177
1178      integer nlev_toga,nt_toga
1179      real ts_toga(nt_toga),plev_toga(nlev_toga,nt_toga)
1180      real t_toga(nlev_toga,nt_toga),q_toga(nlev_toga,nt_toga)
1181      real u_toga(nlev_toga,nt_toga),v_toga(nlev_toga,nt_toga)
1182      real w_toga(nlev_toga,nt_toga)
1183      real ht_toga(nlev_toga,nt_toga),vt_toga(nlev_toga,nt_toga)
1184      real hq_toga(nlev_toga,nt_toga),vq_toga(nlev_toga,nt_toga)
1185      character*80 fich_toga
1186
1187      integer k,ip
1188      real bid
1189
1190      integer iy,im,id,ih
1191     
1192       real plev_min
1193
1194       plev_min = 55.  ! pas de tendance de vap. d eau au-dessus de 55 hPa
1195
1196      open(21,file=trim(fich_toga),form='formatted')
1197      read(21,'(a)') 
1198      do ip = 1, nt_toga
1199      read(21,'(a)') 
1200      read(21,'(a)') 
1201      read(21,223) iy, im, id, ih, bid, ts_toga(ip), bid,bid,bid,bid
1202      read(21,'(a)') 
1203      read(21,'(a)') 
1204
1205       do k = 1, nlev_toga
1206         read(21,230) plev_toga(k,ip), t_toga(k,ip), q_toga(k,ip)          &
1207     &       ,u_toga(k,ip), v_toga(k,ip), w_toga(k,ip)                     &
1208     &       ,ht_toga(k,ip), vt_toga(k,ip), hq_toga(k,ip), vq_toga(k,ip)
1209
1210! conversion in SI units:
1211         t_toga(k,ip)=t_toga(k,ip)+273.15     ! K
1212         q_toga(k,ip)=q_toga(k,ip)*0.001      ! kg/kg
1213         w_toga(k,ip)=w_toga(k,ip)*100./3600. ! Pa/s
1214! no water vapour tendency above 55 hPa
1215         if (plev_toga(k,ip) .lt. plev_min) then
1216          q_toga(k,ip) = 0.
1217          hq_toga(k,ip) = 0.
1218          vq_toga(k,ip) =0.
1219         endif
1220       enddo
1221
1222         ts_toga(ip)=ts_toga(ip)+273.15       ! K
1223       enddo
1224       close(21)
1225
1226  223 format(4i3,6f8.2)
1227  230 format(6f9.3,4e11.3)
1228
1229          return
1230          end
1231
1232!-------------------------------------------------------------------------
1233      SUBROUTINE read_sandu(fich_sandu,nlev_sandu,nt_sandu,ts_sandu)
1234      implicit none
1235
1236!-------------------------------------------------------------------------
1237! Read I.SANDU case forcing data
1238!-------------------------------------------------------------------------
1239
1240      integer nlev_sandu,nt_sandu
1241      real ts_sandu(nt_sandu)
1242      character*80 fich_sandu
1243
1244      integer ip
1245      integer iy,im,id,ih
1246
1247      real plev_min
1248
1249      print*,'nlev_sandu',nlev_sandu
1250      plev_min = 55000.  ! pas de tendance de vap. d eau au-dessus de 55 hPa
1251
1252      open(21,file=trim(fich_sandu),form='formatted')
1253      read(21,'(a)')
1254      do ip = 1, nt_sandu
1255      read(21,'(a)')
1256      read(21,'(a)')
1257      read(21,223) iy, im, id, ih, ts_sandu(ip)
1258      print *,'ts=',iy,im,id,ih,ip,ts_sandu(ip)
1259      enddo
1260      close(21)
1261
1262  223 format(4i3,f8.2)
1263
1264          return
1265          end
1266
1267!=====================================================================
1268!-------------------------------------------------------------------------
1269      SUBROUTINE read_astex(fich_astex,nlev_astex,nt_astex,div_astex,      &
1270     & ts_astex,ug_astex,vg_astex,ufa_astex,vfa_astex)
1271      implicit none
1272
1273!-------------------------------------------------------------------------
1274! Read Astex case forcing data
1275!-------------------------------------------------------------------------
1276
1277      integer nlev_astex,nt_astex
1278      real div_astex(nt_astex),ts_astex(nt_astex),ug_astex(nt_astex)
1279      real vg_astex(nt_astex),ufa_astex(nt_astex),vfa_astex(nt_astex)
1280      character*80 fich_astex
1281
1282      integer ip
1283      integer iy,im,id,ih
1284
1285       real plev_min
1286
1287      print*,'nlev_astex',nlev_astex
1288       plev_min = 55000.  ! pas de tendance de vap. d eau au-dessus de 55 hPa
1289
1290      open(21,file=trim(fich_astex),form='formatted')
1291      read(21,'(a)')
1292      read(21,'(a)')
1293      do ip = 1, nt_astex
1294      read(21,'(a)')
1295      read(21,'(a)')
1296      read(21,223) iy, im, id, ih, div_astex(ip),ts_astex(ip),             &
1297     &ug_astex(ip),vg_astex(ip),ufa_astex(ip),vfa_astex(ip)
1298      ts_astex(ip)=ts_astex(ip)+273.15
1299      print *,'ts=',iy,im,id,ih,ip,div_astex(ip),ts_astex(ip),             &
1300     &ug_astex(ip),vg_astex(ip),ufa_astex(ip),vg_astex(ip)
1301      enddo
1302      close(21)
1303
1304  223 format(4i3,e13.2,f7.2,f7.3,f7.2,f7.3,f7.2)
1305
1306          return
1307          end
1308!=====================================================================
1309      subroutine read_twpice(fich_twpice,nlevel,ntime                       &
1310     &     ,T_srf,plev,T,q,u,v,omega                                       &
1311     &     ,T_adv_h,T_adv_v,q_adv_h,q_adv_v)
1312
1313!program reading forcings of the TWP-ICE experiment
1314
1315!      use netcdf
1316
1317      implicit none
1318
1319#include "netcdf.inc"
1320
1321      integer ntime,nlevel
1322      integer l,k
1323      character*80 :: fich_twpice
1324      real*8 time(ntime)
1325      real*8 lat, lon, alt, phis
1326      real*8 lev(nlevel)
1327      real*8 plev(nlevel,ntime)
1328
1329      real*8 T(nlevel,ntime)
1330      real*8 q(nlevel,ntime),u(nlevel,ntime)
1331      real*8 v(nlevel,ntime)
1332      real*8 omega(nlevel,ntime), div(nlevel,ntime)
1333      real*8 T_adv_h(nlevel,ntime)
1334      real*8 T_adv_v(nlevel,ntime), q_adv_h(nlevel,ntime)
1335      real*8 q_adv_v(nlevel,ntime)
1336      real*8 s(nlevel,ntime), s_adv_h(nlevel,ntime)
1337      real*8 s_adv_v(nlevel,ntime)
1338      real*8 p_srf_aver(ntime), p_srf_center(ntime)
1339      real*8 T_srf(ntime)
1340
1341      integer nid, ierr
1342      integer nbvar3d
1343      parameter(nbvar3d=20)
1344      integer var3didin(nbvar3d)
1345
1346      ierr = NF_OPEN(fich_twpice,NF_NOWRITE,nid)
1347      if (ierr.NE.NF_NOERR) then
1348         write(*,*) 'ERROR: Pb opening forcings cdf file '
1349         write(*,*) NF_STRERROR(ierr)
1350         stop ""
1351      endif
1352
1353      ierr=NF_INQ_VARID(nid,"lat",var3didin(1))
1354         if(ierr/=NF_NOERR) then
1355           write(*,*) NF_STRERROR(ierr)
1356           stop 'lat'
1357         endif
1358     
1359       ierr=NF_INQ_VARID(nid,"lon",var3didin(2))
1360         if(ierr/=NF_NOERR) then
1361           write(*,*) NF_STRERROR(ierr)
1362           stop 'lon'
1363         endif
1364
1365       ierr=NF_INQ_VARID(nid,"alt",var3didin(3))
1366         if(ierr/=NF_NOERR) then
1367           write(*,*) NF_STRERROR(ierr)
1368           stop 'alt'
1369         endif
1370
1371      ierr=NF_INQ_VARID(nid,"phis",var3didin(4))
1372         if(ierr/=NF_NOERR) then
1373           write(*,*) NF_STRERROR(ierr)
1374           stop 'phis'
1375         endif
1376
1377      ierr=NF_INQ_VARID(nid,"T",var3didin(5))
1378         if(ierr/=NF_NOERR) then
1379           write(*,*) NF_STRERROR(ierr)
1380           stop 'T'
1381         endif
1382
1383      ierr=NF_INQ_VARID(nid,"q",var3didin(6))
1384         if(ierr/=NF_NOERR) then
1385           write(*,*) NF_STRERROR(ierr)
1386           stop 'q'
1387         endif
1388
1389      ierr=NF_INQ_VARID(nid,"u",var3didin(7))
1390         if(ierr/=NF_NOERR) then
1391           write(*,*) NF_STRERROR(ierr)
1392           stop 'u'
1393         endif
1394
1395      ierr=NF_INQ_VARID(nid,"v",var3didin(8))
1396         if(ierr/=NF_NOERR) then
1397           write(*,*) NF_STRERROR(ierr)
1398           stop 'v'
1399         endif
1400
1401      ierr=NF_INQ_VARID(nid,"omega",var3didin(9))
1402         if(ierr/=NF_NOERR) then
1403           write(*,*) NF_STRERROR(ierr)
1404           stop 'omega'
1405         endif
1406
1407      ierr=NF_INQ_VARID(nid,"div",var3didin(10))
1408         if(ierr/=NF_NOERR) then
1409           write(*,*) NF_STRERROR(ierr)
1410           stop 'div'
1411         endif
1412
1413      ierr=NF_INQ_VARID(nid,"T_adv_h",var3didin(11))
1414         if(ierr/=NF_NOERR) then
1415           write(*,*) NF_STRERROR(ierr)
1416           stop 'T_adv_h'
1417         endif
1418
1419      ierr=NF_INQ_VARID(nid,"T_adv_v",var3didin(12))
1420         if(ierr/=NF_NOERR) then
1421           write(*,*) NF_STRERROR(ierr)
1422           stop 'T_adv_v'
1423         endif
1424
1425      ierr=NF_INQ_VARID(nid,"q_adv_h",var3didin(13))
1426         if(ierr/=NF_NOERR) then
1427           write(*,*) NF_STRERROR(ierr)
1428           stop 'q_adv_h'
1429         endif
1430
1431      ierr=NF_INQ_VARID(nid,"q_adv_v",var3didin(14))
1432         if(ierr/=NF_NOERR) then
1433           write(*,*) NF_STRERROR(ierr)
1434           stop 'q_adv_v'
1435         endif
1436
1437      ierr=NF_INQ_VARID(nid,"s",var3didin(15))
1438         if(ierr/=NF_NOERR) then
1439           write(*,*) NF_STRERROR(ierr)
1440           stop 's'
1441         endif
1442
1443      ierr=NF_INQ_VARID(nid,"s_adv_h",var3didin(16))
1444         if(ierr/=NF_NOERR) then
1445           write(*,*) NF_STRERROR(ierr)
1446           stop 's_adv_h'
1447         endif
1448   
1449      ierr=NF_INQ_VARID(nid,"s_adv_v",var3didin(17))
1450         if(ierr/=NF_NOERR) then
1451           write(*,*) NF_STRERROR(ierr)
1452           stop 's_adv_v'
1453         endif
1454
1455      ierr=NF_INQ_VARID(nid,"p_srf_aver",var3didin(18))
1456         if(ierr/=NF_NOERR) then
1457           write(*,*) NF_STRERROR(ierr)
1458           stop 'p_srf_aver'
1459         endif
1460
1461      ierr=NF_INQ_VARID(nid,"p_srf_center",var3didin(19))
1462         if(ierr/=NF_NOERR) then
1463           write(*,*) NF_STRERROR(ierr)
1464           stop 'p_srf_center'
1465         endif
1466
1467      ierr=NF_INQ_VARID(nid,"T_srf",var3didin(20))
1468         if(ierr/=NF_NOERR) then
1469           write(*,*) NF_STRERROR(ierr)
1470           stop 'T_srf'
1471         endif
1472
1473!dimensions lecture
1474      call catchaxis(nid,ntime,nlevel,time,lev,ierr)
1475
1476!pressure
1477       do l=1,ntime
1478       do k=1,nlevel
1479          plev(k,l)=lev(k)
1480       enddo
1481       enddo
1482         
1483#ifdef NC_DOUBLE
1484         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(1),lat)
1485#else
1486         ierr = NF_GET_VAR_REAL(nid,var3didin(1),lat)
1487#endif
1488         if(ierr/=NF_NOERR) then
1489            write(*,*) NF_STRERROR(ierr)
1490            stop "getvarup"
1491         endif
1492!         write(*,*)'lecture lat ok',lat
1493
1494#ifdef NC_DOUBLE
1495         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(2),lon)
1496#else
1497         ierr = NF_GET_VAR_REAL(nid,var3didin(2),lon)
1498#endif
1499         if(ierr/=NF_NOERR) then
1500            write(*,*) NF_STRERROR(ierr)
1501            stop "getvarup"
1502         endif
1503!         write(*,*)'lecture lon ok',lon
1504 
1505#ifdef NC_DOUBLE
1506         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(3),alt)
1507#else
1508         ierr = NF_GET_VAR_REAL(nid,var3didin(3),alt)
1509#endif
1510         if(ierr/=NF_NOERR) then
1511            write(*,*) NF_STRERROR(ierr)
1512            stop "getvarup"
1513         endif
1514!          write(*,*)'lecture alt ok',alt
1515 
1516#ifdef NC_DOUBLE
1517         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(4),phis)
1518#else
1519         ierr = NF_GET_VAR_REAL(nid,var3didin(4),phis)
1520#endif
1521         if(ierr/=NF_NOERR) then
1522            write(*,*) NF_STRERROR(ierr)
1523            stop "getvarup"
1524         endif
1525!          write(*,*)'lecture phis ok',phis
1526         
1527#ifdef NC_DOUBLE
1528         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(5),T)
1529#else
1530         ierr = NF_GET_VAR_REAL(nid,var3didin(5),T)
1531#endif
1532         if(ierr/=NF_NOERR) then
1533            write(*,*) NF_STRERROR(ierr)
1534            stop "getvarup"
1535         endif
1536!         write(*,*)'lecture T ok'
1537
1538#ifdef NC_DOUBLE
1539         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(6),q)
1540#else
1541         ierr = NF_GET_VAR_REAL(nid,var3didin(6),q)
1542#endif
1543         if(ierr/=NF_NOERR) then
1544            write(*,*) NF_STRERROR(ierr)
1545            stop "getvarup"
1546         endif
1547!         write(*,*)'lecture q ok'
1548!q in kg/kg
1549       do l=1,ntime
1550       do k=1,nlevel
1551          q(k,l)=q(k,l)/1000.
1552       enddo
1553       enddo
1554#ifdef NC_DOUBLE
1555         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(7),u)
1556#else
1557         ierr = NF_GET_VAR_REAL(nid,var3didin(7),u)
1558#endif
1559         if(ierr/=NF_NOERR) then
1560            write(*,*) NF_STRERROR(ierr)
1561            stop "getvarup"
1562         endif
1563!         write(*,*)'lecture u ok'
1564
1565#ifdef NC_DOUBLE
1566         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(8),v)
1567#else
1568         ierr = NF_GET_VAR_REAL(nid,var3didin(8),v)
1569#endif
1570         if(ierr/=NF_NOERR) then
1571            write(*,*) NF_STRERROR(ierr)
1572            stop "getvarup"
1573         endif
1574!         write(*,*)'lecture v ok'
1575
1576#ifdef NC_DOUBLE
1577         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(9),omega)
1578#else
1579         ierr = NF_GET_VAR_REAL(nid,var3didin(9),omega)
1580#endif
1581         if(ierr/=NF_NOERR) then
1582            write(*,*) NF_STRERROR(ierr)
1583            stop "getvarup"
1584         endif
1585!         write(*,*)'lecture omega ok'
1586!omega in mb/hour
1587       do l=1,ntime
1588       do k=1,nlevel
1589          omega(k,l)=omega(k,l)*100./3600.
1590       enddo
1591       enddo
1592
1593#ifdef NC_DOUBLE
1594         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(10),div)
1595#else
1596         ierr = NF_GET_VAR_REAL(nid,var3didin(10),div)
1597#endif
1598         if(ierr/=NF_NOERR) then
1599            write(*,*) NF_STRERROR(ierr)
1600            stop "getvarup"
1601         endif
1602!         write(*,*)'lecture div ok'
1603
1604#ifdef NC_DOUBLE
1605         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(11),T_adv_h)
1606#else
1607         ierr = NF_GET_VAR_REAL(nid,var3didin(11),T_adv_h)
1608#endif
1609         if(ierr/=NF_NOERR) then
1610            write(*,*) NF_STRERROR(ierr)
1611            stop "getvarup"
1612         endif
1613!         write(*,*)'lecture T_adv_h ok'
1614!T adv in K/s
1615       do l=1,ntime
1616       do k=1,nlevel
1617          T_adv_h(k,l)=T_adv_h(k,l)/3600.
1618       enddo
1619       enddo
1620
1621
1622#ifdef NC_DOUBLE
1623         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(12),T_adv_v)
1624#else
1625         ierr = NF_GET_VAR_REAL(nid,var3didin(12),T_adv_v)
1626#endif
1627         if(ierr/=NF_NOERR) then
1628            write(*,*) NF_STRERROR(ierr)
1629            stop "getvarup"
1630         endif
1631!         write(*,*)'lecture T_adv_v ok'
1632!T adv in K/s
1633       do l=1,ntime
1634       do k=1,nlevel
1635          T_adv_v(k,l)=T_adv_v(k,l)/3600.
1636       enddo
1637       enddo
1638
1639#ifdef NC_DOUBLE
1640         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(13),q_adv_h)
1641#else
1642         ierr = NF_GET_VAR_REAL(nid,var3didin(13),q_adv_h)
1643#endif
1644         if(ierr/=NF_NOERR) then
1645            write(*,*) NF_STRERROR(ierr)
1646            stop "getvarup"
1647         endif
1648!         write(*,*)'lecture q_adv_h ok'
1649!q adv in kg/kg/s
1650       do l=1,ntime
1651       do k=1,nlevel
1652          q_adv_h(k,l)=q_adv_h(k,l)/1000./3600.
1653       enddo
1654       enddo
1655
1656
1657#ifdef NC_DOUBLE
1658         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(14),q_adv_v)
1659#else
1660         ierr = NF_GET_VAR_REAL(nid,var3didin(14),q_adv_v)
1661#endif
1662         if(ierr/=NF_NOERR) then
1663            write(*,*) NF_STRERROR(ierr)
1664            stop "getvarup"
1665         endif
1666!         write(*,*)'lecture q_adv_v ok'
1667!q adv in kg/kg/s
1668       do l=1,ntime
1669       do k=1,nlevel
1670          q_adv_v(k,l)=q_adv_v(k,l)/1000./3600.
1671       enddo
1672       enddo
1673
1674
1675#ifdef NC_DOUBLE
1676         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(15),s)
1677#else
1678         ierr = NF_GET_VAR_REAL(nid,var3didin(15),s)
1679#endif
1680         if(ierr/=NF_NOERR) then
1681            write(*,*) NF_STRERROR(ierr)
1682            stop "getvarup"
1683         endif
1684
1685#ifdef NC_DOUBLE
1686         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(16),s_adv_h)
1687#else
1688         ierr = NF_GET_VAR_REAL(nid,var3didin(16),s_adv_h)
1689#endif
1690         if(ierr/=NF_NOERR) then
1691            write(*,*) NF_STRERROR(ierr)
1692            stop "getvarup"
1693         endif
1694
1695#ifdef NC_DOUBLE
1696         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(17),s_adv_v)
1697#else
1698         ierr = NF_GET_VAR_REAL(nid,var3didin(17),s_adv_v)
1699#endif
1700         if(ierr/=NF_NOERR) then
1701            write(*,*) NF_STRERROR(ierr)
1702            stop "getvarup"
1703         endif
1704
1705#ifdef NC_DOUBLE
1706         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(18),p_srf_aver)
1707#else
1708         ierr = NF_GET_VAR_REAL(nid,var3didin(18),p_srf_aver)
1709#endif
1710         if(ierr/=NF_NOERR) then
1711            write(*,*) NF_STRERROR(ierr)
1712            stop "getvarup"
1713         endif
1714
1715#ifdef NC_DOUBLE
1716         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(19),p_srf_center)
1717#else
1718         ierr = NF_GET_VAR_REAL(nid,var3didin(19),p_srf_center)
1719#endif
1720         if(ierr/=NF_NOERR) then
1721            write(*,*) NF_STRERROR(ierr)
1722            stop "getvarup"
1723         endif
1724
1725#ifdef NC_DOUBLE
1726         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(20),T_srf)
1727#else
1728         ierr = NF_GET_VAR_REAL(nid,var3didin(20),T_srf)
1729#endif
1730         if(ierr/=NF_NOERR) then
1731            write(*,*) NF_STRERROR(ierr)
1732            stop "getvarup"
1733         endif
1734!         write(*,*)'lecture T_srf ok', T_srf
1735
1736         return 
1737         end subroutine read_twpice
1738!=====================================================================
1739         subroutine catchaxis(nid,ttm,llm,time,lev,ierr)
1740
1741!         use netcdf
1742
1743         implicit none
1744#include "netcdf.inc"
1745         integer nid,ttm,llm
1746         real*8 time(ttm)
1747         real*8 lev(llm)
1748         integer ierr
1749
1750         integer timevar,levvar
1751         integer timelen,levlen
1752         integer timedimin,levdimin
1753
1754! Control & lecture on dimensions
1755! ===============================
1756         ierr=NF_INQ_DIMID(nid,"time",timedimin)
1757         ierr=NF_INQ_VARID(nid,"time",timevar)
1758         if (ierr.NE.NF_NOERR) then
1759            write(*,*) 'ERROR: Field <time> is missing'
1760            stop "" 
1761         endif
1762         ierr=NF_INQ_DIMLEN(nid,timedimin,timelen)
1763
1764         ierr=NF_INQ_DIMID(nid,"lev",levdimin)
1765         ierr=NF_INQ_VARID(nid,"lev",levvar)
1766         if (ierr.NE.NF_NOERR) then
1767             write(*,*) 'ERROR: Field <lev> is lacking'
1768             stop "" 
1769         endif
1770         ierr=NF_INQ_DIMLEN(nid,levdimin,levlen)
1771
1772         if((timelen/=ttm).or.(levlen/=llm)) then
1773            write(*,*) 'ERROR: Not the good lenght for axis'
1774            write(*,*) 'longitude: ',timelen,ttm+1
1775            write(*,*) 'latitude: ',levlen,llm
1776            stop "" 
1777         endif
1778
1779!#ifdef NC_DOUBLE
1780         ierr = NF_GET_VAR_DOUBLE(nid,timevar,time)
1781         ierr = NF_GET_VAR_DOUBLE(nid,levvar,lev)
1782!#else
1783!        ierr = NF_GET_VAR_REAL(nid,timevar,time)
1784!        ierr = NF_GET_VAR_REAL(nid,levvar,lev)
1785!#endif
1786
1787       return
1788       end
1789!=====================================================================
1790
1791       SUBROUTINE interp_sandu_vertical(play,nlev_sandu,plev_prof          &
1792     &         ,t_prof,thl_prof,q_prof,u_prof,v_prof,w_prof                &
1793     &         ,omega_prof,o3mmr_prof                                      &
1794     &         ,t_mod,thl_mod,q_mod,u_mod,v_mod,w_mod                      &
1795     &         ,omega_mod,o3mmr_mod,mxcalc)
1796
1797       implicit none
1798
1799#include "dimensions.h"
1800
1801!-------------------------------------------------------------------------
1802! Vertical interpolation of SANDUREF forcing data onto model levels
1803!-------------------------------------------------------------------------
1804
1805       integer nlevmax
1806       parameter (nlevmax=41)
1807       integer nlev_sandu,mxcalc
1808!       real play(llm), plev_prof(nlevmax)
1809!       real t_prof(nlevmax),q_prof(nlevmax)
1810!       real u_prof(nlevmax),v_prof(nlevmax), w_prof(nlevmax)
1811!       real ht_prof(nlevmax),vt_prof(nlevmax)
1812!       real hq_prof(nlevmax),vq_prof(nlevmax)
1813
1814       real play(llm), plev_prof(nlev_sandu)
1815       real t_prof(nlev_sandu),thl_prof(nlev_sandu),q_prof(nlev_sandu)
1816       real u_prof(nlev_sandu),v_prof(nlev_sandu), w_prof(nlev_sandu)
1817       real omega_prof(nlev_sandu),o3mmr_prof(nlev_sandu)
1818
1819       real t_mod(llm),thl_mod(llm),q_mod(llm)
1820       real u_mod(llm),v_mod(llm), w_mod(llm)
1821       real omega_mod(llm),o3mmr_mod(llm)
1822
1823       integer l,k,k1,k2
1824       real frac,frac1,frac2,fact
1825
1826       do l = 1, llm
1827
1828        if (play(l).ge.plev_prof(nlev_sandu)) then
1829
1830        mxcalc=l
1831         k1=0
1832         k2=0
1833
1834         if (play(l).le.plev_prof(1)) then
1835
1836         do k = 1, nlev_sandu-1
1837          if (play(l).le.plev_prof(k).and. play(l).gt.plev_prof(k+1)) then
1838            k1=k
1839            k2=k+1
1840          endif
1841         enddo
1842
1843         if (k1.eq.0 .or. k2.eq.0) then
1844          write(*,*) 'PB! k1, k2 = ',k1,k2
1845          write(*,*) 'l,play(l) = ',l,play(l)/100
1846         do k = 1, nlev_sandu-1
1847          write(*,*) 'k,plev_prof(k) = ',k,plev_prof(k)/100
1848         enddo
1849         endif
1850
1851         frac = (plev_prof(k2)-play(l))/(plev_prof(k2)-plev_prof(k1))
1852         t_mod(l)= t_prof(k2) - frac*(t_prof(k2)-t_prof(k1))
1853         thl_mod(l)= thl_prof(k2) - frac*(thl_prof(k2)-thl_prof(k1))
1854         q_mod(l)= q_prof(k2) - frac*(q_prof(k2)-q_prof(k1))
1855         u_mod(l)= u_prof(k2) - frac*(u_prof(k2)-u_prof(k1))
1856         v_mod(l)= v_prof(k2) - frac*(v_prof(k2)-v_prof(k1))
1857         w_mod(l)= w_prof(k2) - frac*(w_prof(k2)-w_prof(k1))
1858         omega_mod(l)=omega_prof(k2)-frac*(omega_prof(k2)-omega_prof(k1))
1859         o3mmr_mod(l)=o3mmr_prof(k2)-frac*(o3mmr_prof(k2)-o3mmr_prof(k1))
1860
1861         else !play>plev_prof(1)
1862
1863         k1=1
1864         k2=2
1865         frac1 = (play(l)-plev_prof(k2))/(plev_prof(k1)-plev_prof(k2))
1866         frac2 = (play(l)-plev_prof(k1))/(plev_prof(k1)-plev_prof(k2))
1867         t_mod(l)= frac1*t_prof(k1) - frac2*t_prof(k2)
1868         thl_mod(l)= frac1*thl_prof(k1) - frac2*thl_prof(k2)
1869         q_mod(l)= frac1*q_prof(k1) - frac2*q_prof(k2)
1870         u_mod(l)= frac1*u_prof(k1) - frac2*u_prof(k2)
1871         v_mod(l)= frac1*v_prof(k1) - frac2*v_prof(k2)
1872         w_mod(l)= frac1*w_prof(k1) - frac2*w_prof(k2)
1873         omega_mod(l)= frac1*omega_prof(k1) - frac2*omega_prof(k2)
1874         o3mmr_mod(l)= frac1*o3mmr_prof(k1) - frac2*o3mmr_prof(k2)
1875
1876         endif ! play.le.plev_prof(1)
1877
1878        else ! above max altitude of forcing file
1879
1880!jyg
1881         fact=20.*(plev_prof(nlev_sandu)-play(l))/plev_prof(nlev_sandu) !jyg
1882         fact = max(fact,0.)                                           !jyg
1883         fact = exp(-fact)                                             !jyg
1884         t_mod(l)= t_prof(nlev_sandu)                                   !jyg
1885         thl_mod(l)= thl_prof(nlev_sandu)                                   !jyg
1886         q_mod(l)= q_prof(nlev_sandu)*fact                              !jyg
1887         u_mod(l)= u_prof(nlev_sandu)*fact                              !jyg
1888         v_mod(l)= v_prof(nlev_sandu)*fact                              !jyg
1889         w_mod(l)= w_prof(nlev_sandu)*fact                              !jyg
1890         omega_mod(l)= omega_prof(nlev_sandu)*fact                      !jyg
1891         o3mmr_mod(l)= o3mmr_prof(nlev_sandu)*fact                      !jyg
1892
1893        endif ! play
1894
1895       enddo ! l
1896
1897       do l = 1,llm
1898!      print *,'t_mod(l),thl_mod(l),q_mod(l),u_mod(l),v_mod(l) ',
1899!    $        l,t_mod(l),thl_mod(l),q_mod(l),u_mod(l),v_mod(l)
1900       enddo
1901
1902          return
1903          end
1904!=====================================================================
1905       SUBROUTINE interp_astex_vertical(play,nlev_astex,plev_prof          &
1906     &         ,t_prof,thl_prof,qv_prof,ql_prof,qt_prof,u_prof,v_prof      &
1907     &         ,w_prof,tke_prof,o3mmr_prof                                 &
1908     &         ,t_mod,thl_mod,qv_mod,ql_mod,qt_mod,u_mod,v_mod,w_mod       &
1909     &         ,tke_mod,o3mmr_mod,mxcalc)
1910
1911       implicit none
1912
1913#include "dimensions.h"
1914
1915!-------------------------------------------------------------------------
1916! Vertical interpolation of Astex forcing data onto model levels
1917!-------------------------------------------------------------------------
1918
1919       integer nlevmax
1920       parameter (nlevmax=41)
1921       integer nlev_astex,mxcalc
1922!       real play(llm), plev_prof(nlevmax)
1923!       real t_prof(nlevmax),qv_prof(nlevmax)
1924!       real u_prof(nlevmax),v_prof(nlevmax), w_prof(nlevmax)
1925!       real ht_prof(nlevmax),vt_prof(nlevmax)
1926!       real hq_prof(nlevmax),vq_prof(nlevmax)
1927
1928       real play(llm), plev_prof(nlev_astex)
1929       real t_prof(nlev_astex),thl_prof(nlev_astex),qv_prof(nlev_astex)
1930       real u_prof(nlev_astex),v_prof(nlev_astex), w_prof(nlev_astex)
1931       real o3mmr_prof(nlev_astex),ql_prof(nlev_astex)
1932       real qt_prof(nlev_astex),tke_prof(nlev_astex)
1933
1934       real t_mod(llm),thl_mod(llm),qv_mod(llm)
1935       real u_mod(llm),v_mod(llm), w_mod(llm),tke_mod(llm)
1936       real o3mmr_mod(llm),ql_mod(llm),qt_mod(llm)
1937
1938       integer l,k,k1,k2
1939       real frac,frac1,frac2,fact
1940
1941       do l = 1, llm
1942
1943        if (play(l).ge.plev_prof(nlev_astex)) then
1944
1945        mxcalc=l
1946         k1=0
1947         k2=0
1948
1949         if (play(l).le.plev_prof(1)) then
1950
1951         do k = 1, nlev_astex-1
1952          if (play(l).le.plev_prof(k).and. play(l).gt.plev_prof(k+1)) then
1953            k1=k
1954            k2=k+1
1955          endif
1956         enddo
1957
1958         if (k1.eq.0 .or. k2.eq.0) then
1959          write(*,*) 'PB! k1, k2 = ',k1,k2
1960          write(*,*) 'l,play(l) = ',l,play(l)/100
1961         do k = 1, nlev_astex-1
1962          write(*,*) 'k,plev_prof(k) = ',k,plev_prof(k)/100
1963         enddo
1964         endif
1965
1966         frac = (plev_prof(k2)-play(l))/(plev_prof(k2)-plev_prof(k1))
1967         t_mod(l)= t_prof(k2) - frac*(t_prof(k2)-t_prof(k1))
1968         thl_mod(l)= thl_prof(k2) - frac*(thl_prof(k2)-thl_prof(k1))
1969         qv_mod(l)= qv_prof(k2) - frac*(qv_prof(k2)-qv_prof(k1))
1970         ql_mod(l)= ql_prof(k2) - frac*(ql_prof(k2)-ql_prof(k1))
1971         qt_mod(l)= qt_prof(k2) - frac*(qt_prof(k2)-qt_prof(k1))
1972         u_mod(l)= u_prof(k2) - frac*(u_prof(k2)-u_prof(k1))
1973         v_mod(l)= v_prof(k2) - frac*(v_prof(k2)-v_prof(k1))
1974         w_mod(l)= w_prof(k2) - frac*(w_prof(k2)-w_prof(k1))
1975         tke_mod(l)= tke_prof(k2) - frac*(tke_prof(k2)-tke_prof(k1))
1976         o3mmr_mod(l)=o3mmr_prof(k2)-frac*(o3mmr_prof(k2)-o3mmr_prof(k1))
1977
1978         else !play>plev_prof(1)
1979
1980         k1=1
1981         k2=2
1982         frac1 = (play(l)-plev_prof(k2))/(plev_prof(k1)-plev_prof(k2))
1983         frac2 = (play(l)-plev_prof(k1))/(plev_prof(k1)-plev_prof(k2))
1984         t_mod(l)= frac1*t_prof(k1) - frac2*t_prof(k2)
1985         thl_mod(l)= frac1*thl_prof(k1) - frac2*thl_prof(k2)
1986         qv_mod(l)= frac1*qv_prof(k1) - frac2*qv_prof(k2)
1987         ql_mod(l)= frac1*ql_prof(k1) - frac2*ql_prof(k2)
1988         qt_mod(l)= frac1*qt_prof(k1) - frac2*qt_prof(k2)
1989         u_mod(l)= frac1*u_prof(k1) - frac2*u_prof(k2)
1990         v_mod(l)= frac1*v_prof(k1) - frac2*v_prof(k2)
1991         w_mod(l)= frac1*w_prof(k1) - frac2*w_prof(k2)
1992         tke_mod(l)= frac1*tke_prof(k1) - frac2*tke_prof(k2)
1993         o3mmr_mod(l)= frac1*o3mmr_prof(k1) - frac2*o3mmr_prof(k2)
1994
1995         endif ! play.le.plev_prof(1)
1996
1997        else ! above max altitude of forcing file
1998
1999!jyg
2000         fact=20.*(plev_prof(nlev_astex)-play(l))/plev_prof(nlev_astex) !jyg
2001         fact = max(fact,0.)                                           !jyg
2002         fact = exp(-fact)                                             !jyg
2003         t_mod(l)= t_prof(nlev_astex)                                   !jyg
2004         thl_mod(l)= thl_prof(nlev_astex)                                   !jyg
2005         qv_mod(l)= qv_prof(nlev_astex)*fact                              !jyg
2006         ql_mod(l)= ql_prof(nlev_astex)*fact                              !jyg
2007         qt_mod(l)= qt_prof(nlev_astex)*fact                              !jyg
2008         u_mod(l)= u_prof(nlev_astex)*fact                              !jyg
2009         v_mod(l)= v_prof(nlev_astex)*fact                              !jyg
2010         w_mod(l)= w_prof(nlev_astex)*fact                              !jyg
2011         tke_mod(l)= tke_prof(nlev_astex)*fact                              !jyg
2012         o3mmr_mod(l)= o3mmr_prof(nlev_astex)*fact                      !jyg
2013
2014        endif ! play
2015
2016       enddo ! l
2017
2018       do l = 1,llm
2019!      print *,'t_mod(l),thl_mod(l),qv_mod(l),u_mod(l),v_mod(l) ',
2020!    $        l,t_mod(l),thl_mod(l),qv_mod(l),u_mod(l),v_mod(l)
2021       enddo
2022
2023          return
2024          end
2025
2026!======================================================================
2027      SUBROUTINE read_rico(fich_rico,nlev_rico,ps_rico,play                &
2028     &             ,ts_rico,t_rico,q_rico,u_rico,v_rico,w_rico             &
2029     &             ,dth_dyn,dqh_dyn)
2030      implicit none
2031
2032!-------------------------------------------------------------------------
2033! Read RICO forcing data
2034!-------------------------------------------------------------------------
2035#include "dimensions.h"
2036
2037
2038      integer nlev_rico
2039      real ts_rico,ps_rico
2040      real t_rico(llm),q_rico(llm)
2041      real u_rico(llm),v_rico(llm)
2042      real w_rico(llm)
2043      real dth_dyn(llm)
2044      real dqh_dyn(llm)
2045     
2046
2047      real play(llm),zlay(llm)
2048     
2049
2050      real prico(nlev_rico),zrico(nlev_rico)
2051
2052      character*80 fich_rico
2053
2054      integer k,l
2055
2056     
2057      print*,fich_rico
2058      open(21,file=trim(fich_rico),form='formatted')
2059        do k=1,llm
2060      zlay(k)=0.
2061         enddo
2062     
2063        read(21,*) ps_rico,ts_rico
2064        prico(1)=ps_rico
2065        zrico(1)=0.0
2066      do l=2,nlev_rico
2067        read(21,*) k,prico(l),zrico(l)
2068      enddo
2069       close(21)
2070
2071      do k=1,llm
2072        do l=1,80
2073          if(prico(l)>play(k)) then
2074              if(play(k)>prico(l+1)) then
2075                zlay(k)=zrico(l)+(play(k)-prico(l)) *                      &
2076     &              (zrico(l+1)-zrico(l))/(prico(l+1)-prico(l))
2077              else 
2078                zlay(k)=zrico(l)+(play(k)-prico(80))*                      &
2079     &              (zrico(81)-zrico(80))/(prico(81)-prico(80))
2080              endif
2081          endif
2082        enddo
2083        print*,k,zlay(k)
2084        ! U
2085        if(0 < zlay(k) .and. zlay(k) < 4000) then
2086          u_rico(k)=-9.9 + (-1.9 + 9.9)*zlay(k)/4000
2087        elseif(4000 < zlay(k) .and. zlay(k) < 12000) then
2088       u_rico(k)=  -1.9 + (30.0 + 1.9) /                                   &
2089     &          (12000 - 4000) * (zlay(k) - 4000)
2090        elseif(12000 < zlay(k) .and. zlay(k) < 13000) then
2091          u_rico(k)=30.0
2092        elseif(13000 < zlay(k) .and. zlay(k) < 20000) then
2093          u_rico(k)=30.0 - (30.0) /                                        &
2094     & (20000 - 13000) * (zlay(k) - 13000)
2095        else
2096          u_rico(k)=0.0
2097        endif
2098
2099!Q_v
2100        if(0 < zlay(k) .and. zlay(k) < 740) then
2101          q_rico(k)=16.0 + (13.8 - 16.0) / (740) * zlay(k)
2102        elseif(740 < zlay(k) .and. zlay(k) < 3260) then
2103          q_rico(k)=13.8 + (2.4 - 13.8) /                                   &
2104     &          (3260 - 740) * (zlay(k) - 740)
2105        elseif(3260 < zlay(k) .and. zlay(k) < 4000) then
2106          q_rico(k)=2.4 + (1.8 - 2.4) /                                    &
2107     &               (4000 - 3260) * (zlay(k) - 3260)
2108        elseif(4000 < zlay(k) .and. zlay(k) < 9000) then
2109          q_rico(k)=1.8 + (0 - 1.8) /                                      &
2110     &             (10000 - 4000) * (zlay(k) - 4000)
2111        else
2112          q_rico(k)=0.0
2113        endif
2114
2115!T
2116        if(0 < zlay(k) .and. zlay(k) < 740) then
2117          t_rico(k)=299.2 + (292.0 - 299.2) / (740) * zlay(k)
2118        elseif(740 < zlay(k) .and. zlay(k) < 4000) then
2119          t_rico(k)=292.0 + (278.0 - 292.0) /                              &                       
2120     &       (4000 - 740) * (zlay(k) - 740)
2121        elseif(4000 < zlay(k) .and. zlay(k) < 15000) then
2122          t_rico(k)=278.0 + (203.0 - 278.0) /                              &
2123     &       (15000 - 4000) * (zlay(k) - 4000)
2124        elseif(15000 < zlay(k) .and. zlay(k) < 17500) then
2125          t_rico(k)=203.0 + (194.0 - 203.0) /                              & 
2126     &       (17500 - 15000)* (zlay(k) - 15000)
2127        elseif(17500 < zlay(k) .and. zlay(k) < 20000) then
2128          t_rico(k)=194.0 + (206.0 - 194.0) /                              &
2129     &       (20000 - 17500)* (zlay(k) - 17500)
2130        elseif(20000 < zlay(k) .and. zlay(k) < 60000) then
2131          t_rico(k)=206.0 + (270.0 - 206.0) /                              & 
2132     &        (60000 - 20000)* (zlay(k) - 20000)
2133        endif
2134
2135! W
2136        if(0 < zlay(k) .and. zlay(k) < 2260 ) then
2137          w_rico(k)=- (0.005/2260) * zlay(k)
2138        elseif(2260 < zlay(k) .and. zlay(k) < 4000 ) then
2139          w_rico(k)=- 0.005
2140        elseif(4000 < zlay(k) .and. zlay(k) < 5000 ) then
2141       w_rico(k)=- 0.005 + (0.005/ (5000 - 4000)) * (zlay(k) - 4000)
2142        else
2143          w_rico(k)=0.0
2144        endif
2145
2146! dThrz+dTsw0+dTlw0
2147        if(0 < zlay(k) .and. zlay(k) < 4000) then
2148          dth_dyn(k)=- 2.51 / 86400 + (-2.18 + 2.51 )/                     &
2149     &               (86400*4000) * zlay(k)
2150        elseif(4000 < zlay(k) .and. zlay(k) < 5000) then
2151          dth_dyn(k)=- 2.18 / 86400 + ( 2.18 ) /                           &
2152     &           (86400*(5000 - 4000)) * (zlay(k) - 4000)
2153        else
2154          dth_dyn(k)=0.0
2155        endif
2156! dQhrz
2157        if(0 < zlay(k) .and. zlay(k) < 3000) then
2158          dqh_dyn(k)=-1.0 / 86400 + (0.345 + 1.0)/                         &
2159     &                    (86400*3000) * (zlay(k))
2160        elseif(3000 < zlay(k) .and. zlay(k) < 4000) then
2161          dqh_dyn(k)=0.345 / 86400
2162        elseif(4000 < zlay(k) .and. zlay(k) < 5000) then
2163          dqh_dyn(k)=0.345 / 86400 +                                       &
2164     &   (-0.345)/(86400 * (5000 - 4000)) * (zlay(k)-4000)
2165        else
2166          dqh_dyn(k)=0.0
2167        endif
2168
2169!?        if(play(k)>6e4) then
2170!?          ratqs0(1,k)=ratqsbas*(plev(1)-play(k))/(plev(1)-6e4)
2171!?        elseif((play(k)>3e4).and.(play(k)<6e4)) then
2172!?          ratqs0(1,k)=ratqsbas+(ratqshaut-ratqsbas)&
2173!?                          *(6e4-play(k))/(6e4-3e4)
2174!?        else
2175!?          ratqs0(1,k)=ratqshaut
2176!?        endif
2177
2178      enddo
2179
2180      do k=1,llm
2181      q_rico(k)=q_rico(k)/1e3 
2182      dqh_dyn(k)=dqh_dyn(k)/1e3
2183      v_rico(k)=-3.8
2184      enddo
2185
2186          return
2187          end
2188
2189!======================================================================
2190        SUBROUTINE interp_sandu_time(day,day1,annee_ref                    &
2191     &             ,year_ini_sandu,day_ini_sandu,nt_sandu,dt_sandu         &
2192     &             ,nlev_sandu,ts_sandu,ts_prof)
2193        implicit none
2194
2195!---------------------------------------------------------------------------------------
2196! Time interpolation of a 2D field to the timestep corresponding to day
2197!
2198! day: current julian day (e.g. 717538.2)
2199! day1: first day of the simulation
2200! nt_sandu: total nb of data in the forcing (e.g. 13 for Sanduref)
2201! dt_sandu: total time interval (in sec) between 2 forcing data (e.g. 6h for Sanduref)
2202!---------------------------------------------------------------------------------------
2203! inputs:
2204        integer annee_ref
2205        integer nt_sandu,nlev_sandu
2206        integer year_ini_sandu
2207        real day, day1,day_ini_sandu,dt_sandu
2208        real ts_sandu(nt_sandu)
2209! outputs:
2210        real ts_prof
2211! local:
2212        integer it_sandu1, it_sandu2
2213        real timeit,time_sandu1,time_sandu2,frac
2214! Check that initial day of the simulation consistent with SANDU period:
2215       if (annee_ref.ne.2006 ) then
2216        print*,'Pour SANDUREF, annee_ref doit etre 2006 '
2217        print*,'Changer annee_ref dans run.def'
2218        stop
2219       endif
2220!      if (annee_ref.eq.2006 .and. day1.lt.day_ini_sandu) then
2221!       print*,'SANDUREF debute le 15 Juillet 2006 (jour julien=196)'
2222!       print*,'Changer dayref dans run.def'
2223!       stop
2224!      endif
2225
2226! Determine timestep relative to the 1st day of TOGA-COARE:
2227!       timeit=(day-day1)*86400.
2228!       if (annee_ref.eq.1992) then
2229!        timeit=(day-day_ini_sandu)*86400.
2230!       else
2231!        timeit=(day+61.-1.)*86400. ! 61 days between Nov01 and Dec31 1992
2232!       endif
2233      timeit=(day-day_ini_sandu)*86400
2234
2235! Determine the closest observation times:
2236       it_sandu1=INT(timeit/dt_sandu)+1
2237       it_sandu2=it_sandu1 + 1
2238       time_sandu1=(it_sandu1-1)*dt_sandu
2239       time_sandu2=(it_sandu2-1)*dt_sandu
2240       print *,'timeit day day_ini_sandu',timeit,day,day_ini_sandu
2241       print *,'it_sandu1,it_sandu2,time_sandu1,time_sandu2',              &
2242     &          it_sandu1,it_sandu2,time_sandu1,time_sandu2
2243
2244       if (it_sandu1 .ge. nt_sandu) then
2245        write(*,*) 'PB-stop: day, it_sandu1, it_sandu2, timeit: '          &
2246     &        ,day,it_sandu1,it_sandu2,timeit/86400.
2247        stop
2248       endif
2249
2250! time interpolation:
2251       frac=(time_sandu2-timeit)/(time_sandu2-time_sandu1)
2252       frac=max(frac,0.0)
2253
2254       ts_prof = ts_sandu(it_sandu2)                                       &
2255     &          -frac*(ts_sandu(it_sandu2)-ts_sandu(it_sandu1))
2256
2257         print*,                                                           &
2258     &'day,annee_ref,day_ini_sandu,timeit,it_sandu1,it_sandu2,SST:',       &
2259     &day,annee_ref,day_ini_sandu,timeit/86400.,it_sandu1,                  &
2260     &it_sandu2,ts_prof
2261
2262        return
2263        END
2264!=====================================================================
2265!-------------------------------------------------------------------------
2266      SUBROUTINE read_armcu(fich_armcu,nlev_armcu,nt_armcu,                &
2267     & sens,flat,adv_theta,rad_theta,adv_qt)
2268      implicit none
2269
2270!-------------------------------------------------------------------------
2271! Read ARM_CU case forcing data
2272!-------------------------------------------------------------------------
2273
2274      integer nlev_armcu,nt_armcu
2275      real sens(nt_armcu),flat(nt_armcu)
2276      real adv_theta(nt_armcu),rad_theta(nt_armcu),adv_qt(nt_armcu)
2277      character*80 fich_armcu
2278
2279      integer ip
2280
2281      integer iy,im,id,ih,in
2282
2283      print*,'nlev_armcu',nlev_armcu
2284
2285      open(21,file=trim(fich_armcu),form='formatted')
2286      read(21,'(a)')
2287      do ip = 1, nt_armcu
2288      read(21,'(a)')
2289      read(21,'(a)')
2290      read(21,223) iy, im, id, ih, in, sens(ip),flat(ip),                  &
2291     &             adv_theta(ip),rad_theta(ip),adv_qt(ip)
2292      print *,'forcages=',iy,im,id,ih,in, sens(ip),flat(ip),               &
2293     &             adv_theta(ip),rad_theta(ip),adv_qt(ip)
2294      enddo
2295      close(21)
2296
2297  223 format(5i3,5f8.3)
2298
2299          return
2300          end
2301
2302!=====================================================================
2303       SUBROUTINE interp_toga_vertical(play,nlev_toga,plev_prof            &
2304     &         ,t_prof,q_prof,u_prof,v_prof,w_prof                         &
2305     &         ,ht_prof,vt_prof,hq_prof,vq_prof                            &
2306     &         ,t_mod,q_mod,u_mod,v_mod,w_mod                              &
2307     &         ,ht_mod,vt_mod,hq_mod,vq_mod,mxcalc)
2308 
2309       implicit none
2310 
2311#include "dimensions.h"
2312
2313!-------------------------------------------------------------------------
2314! Vertical interpolation of TOGA-COARE forcing data onto model levels
2315!-------------------------------------------------------------------------
2316 
2317       integer nlevmax
2318       parameter (nlevmax=41)
2319       integer nlev_toga,mxcalc
2320!       real play(llm), plev_prof(nlevmax) 
2321!       real t_prof(nlevmax),q_prof(nlevmax)
2322!       real u_prof(nlevmax),v_prof(nlevmax), w_prof(nlevmax)
2323!       real ht_prof(nlevmax),vt_prof(nlevmax)
2324!       real hq_prof(nlevmax),vq_prof(nlevmax)
2325 
2326       real play(llm), plev_prof(nlev_toga) 
2327       real t_prof(nlev_toga),q_prof(nlev_toga)
2328       real u_prof(nlev_toga),v_prof(nlev_toga), w_prof(nlev_toga)
2329       real ht_prof(nlev_toga),vt_prof(nlev_toga)
2330       real hq_prof(nlev_toga),vq_prof(nlev_toga)
2331 
2332       real t_mod(llm),q_mod(llm)
2333       real u_mod(llm),v_mod(llm), w_mod(llm)
2334       real ht_mod(llm),vt_mod(llm)
2335       real hq_mod(llm),vq_mod(llm)
2336 
2337       integer l,k,k1,k2
2338       real frac,frac1,frac2,fact
2339 
2340       do l = 1, llm
2341
2342        if (play(l).ge.plev_prof(nlev_toga)) then
2343 
2344        mxcalc=l
2345         k1=0
2346         k2=0
2347
2348         if (play(l).le.plev_prof(1)) then
2349
2350         do k = 1, nlev_toga-1
2351          if (play(l).le.plev_prof(k).and. play(l).gt.plev_prof(k+1)) then
2352            k1=k
2353            k2=k+1
2354          endif
2355         enddo
2356
2357         if (k1.eq.0 .or. k2.eq.0) then
2358          write(*,*) 'PB! k1, k2 = ',k1,k2
2359          write(*,*) 'l,play(l) = ',l,play(l)/100
2360         do k = 1, nlev_toga-1
2361          write(*,*) 'k,plev_prof(k) = ',k,plev_prof(k)/100
2362         enddo
2363         endif
2364
2365         frac = (plev_prof(k2)-play(l))/(plev_prof(k2)-plev_prof(k1))
2366         t_mod(l)= t_prof(k2) - frac*(t_prof(k2)-t_prof(k1))
2367         q_mod(l)= q_prof(k2) - frac*(q_prof(k2)-q_prof(k1))
2368         u_mod(l)= u_prof(k2) - frac*(u_prof(k2)-u_prof(k1))
2369         v_mod(l)= v_prof(k2) - frac*(v_prof(k2)-v_prof(k1))
2370         w_mod(l)= w_prof(k2) - frac*(w_prof(k2)-w_prof(k1))
2371         ht_mod(l)= ht_prof(k2) - frac*(ht_prof(k2)-ht_prof(k1))
2372         vt_mod(l)= vt_prof(k2) - frac*(vt_prof(k2)-vt_prof(k1))
2373         hq_mod(l)= hq_prof(k2) - frac*(hq_prof(k2)-hq_prof(k1))
2374         vq_mod(l)= vq_prof(k2) - frac*(vq_prof(k2)-vq_prof(k1))
2375     
2376         else !play>plev_prof(1)
2377
2378         k1=1
2379         k2=2
2380         frac1 = (play(l)-plev_prof(k2))/(plev_prof(k1)-plev_prof(k2))
2381         frac2 = (play(l)-plev_prof(k1))/(plev_prof(k1)-plev_prof(k2))
2382         t_mod(l)= frac1*t_prof(k1) - frac2*t_prof(k2)
2383         q_mod(l)= frac1*q_prof(k1) - frac2*q_prof(k2)
2384         u_mod(l)= frac1*u_prof(k1) - frac2*u_prof(k2)
2385         v_mod(l)= frac1*v_prof(k1) - frac2*v_prof(k2)
2386         w_mod(l)= frac1*w_prof(k1) - frac2*w_prof(k2)
2387         ht_mod(l)= frac1*ht_prof(k1) - frac2*ht_prof(k2)
2388         vt_mod(l)= frac1*vt_prof(k1) - frac2*vt_prof(k2)
2389         hq_mod(l)= frac1*hq_prof(k1) - frac2*hq_prof(k2)
2390         vq_mod(l)= frac1*vq_prof(k1) - frac2*vq_prof(k2)
2391
2392         endif ! play.le.plev_prof(1)
2393
2394        else ! above max altitude of forcing file
2395 
2396!jyg
2397         fact=20.*(plev_prof(nlev_toga)-play(l))/plev_prof(nlev_toga) !jyg
2398         fact = max(fact,0.)                                           !jyg
2399         fact = exp(-fact)                                             !jyg
2400         t_mod(l)= t_prof(nlev_toga)                                   !jyg
2401         q_mod(l)= q_prof(nlev_toga)*fact                              !jyg
2402         u_mod(l)= u_prof(nlev_toga)*fact                              !jyg
2403         v_mod(l)= v_prof(nlev_toga)*fact                              !jyg
2404         w_mod(l)= 0.0                                                 !jyg
2405         ht_mod(l)= ht_prof(nlev_toga)                                 !jyg
2406         vt_mod(l)= vt_prof(nlev_toga)                                 !jyg
2407         hq_mod(l)= hq_prof(nlev_toga)*fact                            !jyg
2408         vq_mod(l)= vq_prof(nlev_toga)*fact                            !jyg
2409 
2410        endif ! play
2411 
2412       enddo ! l
2413
2414!       do l = 1,llm
2415!       print *,'t_mod(l),q_mod(l),ht_mod(l),hq_mod(l) ',
2416!     $        l,t_mod(l),q_mod(l),ht_mod(l),hq_mod(l)
2417!       enddo
2418 
2419          return
2420          end
2421 
2422!======================================================================
2423        SUBROUTINE interp_astex_time(day,day1,annee_ref                    &
2424     &             ,year_ini_astex,day_ini_astex,nt_astex,dt_astex         &
2425     &             ,nlev_astex,div_astex,ts_astex,ug_astex,vg_astex        &
2426     &             ,ufa_astex,vfa_astex,div_prof,ts_prof,ug_prof,vg_prof   &
2427     &             ,ufa_prof,vfa_prof)
2428        implicit none
2429
2430!---------------------------------------------------------------------------------------
2431! Time interpolation of a 2D field to the timestep corresponding to day
2432!
2433! day: current julian day (e.g. 717538.2)
2434! day1: first day of the simulation
2435! nt_astex: total nb of data in the forcing (e.g. 41 for Astex)
2436! dt_astex: total time interval (in sec) between 2 forcing data (e.g. 1h for Astex)
2437!---------------------------------------------------------------------------------------
2438
2439! inputs:
2440        integer annee_ref
2441        integer nt_astex,nlev_astex
2442        integer year_ini_astex
2443        real day, day1,day_ini_astex,dt_astex
2444        real div_astex(nt_astex),ts_astex(nt_astex),ug_astex(nt_astex)
2445        real vg_astex(nt_astex),ufa_astex(nt_astex),vfa_astex(nt_astex)
2446! outputs:
2447        real div_prof,ts_prof,ug_prof,vg_prof,ufa_prof,vfa_prof
2448! local:
2449        integer it_astex1, it_astex2
2450        real timeit,time_astex1,time_astex2,frac
2451
2452! Check that initial day of the simulation consistent with ASTEX period:
2453       if (annee_ref.ne.1992 ) then
2454        print*,'Pour Astex, annee_ref doit etre 1992 '
2455        print*,'Changer annee_ref dans run.def'
2456        stop
2457       endif
2458       if (annee_ref.eq.1992 .and. day1.lt.day_ini_astex) then
2459        print*,'Astex debute le 13 Juin 1992 (jour julien=165)'
2460        print*,'Changer dayref dans run.def'
2461        stop
2462       endif
2463
2464! Determine timestep relative to the 1st day of TOGA-COARE:
2465!       timeit=(day-day1)*86400.
2466!       if (annee_ref.eq.1992) then
2467!        timeit=(day-day_ini_astex)*86400.
2468!       else
2469!        timeit=(day+2.-1.)*86400. ! 2 days between Jun13 and Jun15 1992
2470!       endif
2471      timeit=(day-day_ini_astex)*86400
2472
2473! Determine the closest observation times:
2474       it_astex1=INT(timeit/dt_astex)+1
2475       it_astex2=it_astex1 + 1
2476       time_astex1=(it_astex1-1)*dt_astex
2477       time_astex2=(it_astex2-1)*dt_astex
2478       print *,'timeit day day_ini_astex',timeit,day,day_ini_astex
2479       print *,'it_astex1,it_astex2,time_astex1,time_astex2',              &
2480     &          it_astex1,it_astex2,time_astex1,time_astex2
2481
2482       if (it_astex1 .ge. nt_astex) then
2483        write(*,*) 'PB-stop: day, it_astex1, it_astex2, timeit: '          &
2484     &        ,day,it_astex1,it_astex2,timeit/86400.
2485        stop
2486       endif
2487
2488! time interpolation:
2489       frac=(time_astex2-timeit)/(time_astex2-time_astex1)
2490       frac=max(frac,0.0)
2491
2492       div_prof = div_astex(it_astex2)                                     &
2493     &          -frac*(div_astex(it_astex2)-div_astex(it_astex1))
2494       ts_prof = ts_astex(it_astex2)                                        &
2495     &          -frac*(ts_astex(it_astex2)-ts_astex(it_astex1))
2496       ug_prof = ug_astex(it_astex2)                                       &
2497     &          -frac*(ug_astex(it_astex2)-ug_astex(it_astex1))
2498       vg_prof = vg_astex(it_astex2)                                       &
2499     &          -frac*(vg_astex(it_astex2)-vg_astex(it_astex1))
2500       ufa_prof = ufa_astex(it_astex2)                                     &
2501     &          -frac*(ufa_astex(it_astex2)-ufa_astex(it_astex1))
2502       vfa_prof = vfa_astex(it_astex2)                                     &
2503     &          -frac*(vfa_astex(it_astex2)-vfa_astex(it_astex1))
2504
2505         print*,                                                           &
2506     &'day,annee_ref,day_ini_astex,timeit,it_astex1,it_astex2,SST:',       &
2507     &day,annee_ref,day_ini_astex,timeit/86400.,it_astex1,                 &
2508     &it_astex2,div_prof,ts_prof,ug_prof,vg_prof,ufa_prof,vfa_prof
2509
2510        return
2511        END
2512
2513!======================================================================
2514        SUBROUTINE interp_toga_time(day,day1,annee_ref                     &
2515     &             ,year_ini_toga,day_ini_toga,nt_toga,dt_toga,nlev_toga   &
2516     &             ,ts_toga,plev_toga,t_toga,q_toga,u_toga,v_toga,w_toga   &
2517     &             ,ht_toga,vt_toga,hq_toga,vq_toga                        &
2518     &             ,ts_prof,plev_prof,t_prof,q_prof,u_prof,v_prof,w_prof   &
2519     &             ,ht_prof,vt_prof,hq_prof,vq_prof)
2520        implicit none
2521
2522!---------------------------------------------------------------------------------------
2523! Time interpolation of a 2D field to the timestep corresponding to day
2524!
2525! day: current julian day (e.g. 717538.2)
2526! day1: first day of the simulation
2527! nt_toga: total nb of data in the forcing (e.g. 480 for TOGA-COARE)
2528! dt_toga: total time interval (in sec) between 2 forcing data (e.g. 6h for TOGA-COARE)
2529!---------------------------------------------------------------------------------------
2530
2531#include "compar1d.h"
2532
2533! inputs:
2534        integer annee_ref
2535        integer nt_toga,nlev_toga
2536        integer year_ini_toga
2537        real day, day1,day_ini_toga,dt_toga
2538        real ts_toga(nt_toga)
2539        real plev_toga(nlev_toga,nt_toga),t_toga(nlev_toga,nt_toga)
2540        real q_toga(nlev_toga,nt_toga),u_toga(nlev_toga,nt_toga)
2541        real v_toga(nlev_toga,nt_toga),w_toga(nlev_toga,nt_toga)
2542        real ht_toga(nlev_toga,nt_toga),vt_toga(nlev_toga,nt_toga)
2543        real hq_toga(nlev_toga,nt_toga),vq_toga(nlev_toga,nt_toga)
2544! outputs:
2545        real ts_prof
2546        real plev_prof(nlev_toga),t_prof(nlev_toga)
2547        real q_prof(nlev_toga),u_prof(nlev_toga)
2548        real v_prof(nlev_toga),w_prof(nlev_toga)
2549        real ht_prof(nlev_toga),vt_prof(nlev_toga)
2550        real hq_prof(nlev_toga),vq_prof(nlev_toga)
2551! local:
2552        integer it_toga1, it_toga2,k
2553        real timeit,time_toga1,time_toga2,frac
2554
2555
2556        if (forcing_type.eq.2) then
2557! Check that initial day of the simulation consistent with TOGA-COARE period:
2558       if (annee_ref.ne.1992 .and. annee_ref.ne.1993) then
2559        print*,'Pour TOGA-COARE, annee_ref doit etre 1992 ou 1993'
2560        print*,'Changer annee_ref dans run.def'
2561        stop
2562       endif
2563       if (annee_ref.eq.1992 .and. day1.lt.day_ini_toga) then
2564        print*,'TOGA-COARE a débuté le 1er Nov 1992 (jour julien=306)'
2565        print*,'Changer dayref dans run.def'
2566        stop
2567       endif
2568       if (annee_ref.eq.1993 .and. day1.gt.day_ini_toga+119) then
2569        print*,'TOGA-COARE a fini le 28 Feb 1993 (jour julien=59)'
2570        print*,'Changer dayref ou nday dans run.def'
2571        stop
2572       endif
2573
2574       else if (forcing_type.eq.4) then
2575
2576! Check that initial day of the simulation consistent with TWP-ICE period:
2577       if (annee_ref.ne.2006) then
2578        print*,'Pour TWP-ICE, annee_ref doit etre 2006'
2579        print*,'Changer annee_ref dans run.def'
2580        stop
2581       endif
2582       if (annee_ref.eq.2006 .and. day1.lt.day_ini_toga) then
2583        print*,'TWP-ICE a debute le 17 Jan 2006 (jour julien=17)'
2584        print*,'Changer dayref dans run.def'
2585        stop
2586       endif
2587       if (annee_ref.eq.2006 .and. day1.gt.day_ini_toga+26) then
2588        print*,'TWP-ICE a fini le 12 Feb 2006 (jour julien=43)'
2589        print*,'Changer dayref ou nday dans run.def'
2590        stop
2591       endif
2592
2593       endif
2594
2595! Determine timestep relative to the 1st day of TOGA-COARE:
2596!       timeit=(day-day1)*86400.
2597!       if (annee_ref.eq.1992) then
2598!        timeit=(day-day_ini_toga)*86400.
2599!       else
2600!        timeit=(day+61.-1.)*86400. ! 61 days between Nov01 and Dec31 1992
2601!       endif
2602      timeit=(day-day_ini_toga)*86400
2603
2604! Determine the closest observation times:
2605       it_toga1=INT(timeit/dt_toga)+1
2606       it_toga2=it_toga1 + 1
2607       time_toga1=(it_toga1-1)*dt_toga
2608       time_toga2=(it_toga2-1)*dt_toga
2609
2610       if (it_toga1 .ge. nt_toga) then
2611        write(*,*) 'PB-stop: day, it_toga1, it_toga2, timeit: '            &
2612     &        ,day,it_toga1,it_toga2,timeit/86400.
2613        stop
2614       endif
2615
2616! time interpolation:
2617       frac=(time_toga2-timeit)/(time_toga2-time_toga1)
2618       frac=max(frac,0.0)
2619
2620       ts_prof = ts_toga(it_toga2)                                         &
2621     &          -frac*(ts_toga(it_toga2)-ts_toga(it_toga1))
2622
2623!        print*,
2624!     :'day,annee_ref,day_ini_toga,timeit,it_toga1,it_toga2,SST:',
2625!     :day,annee_ref,day_ini_toga,timeit/86400.,it_toga1,it_toga2,ts_prof
2626
2627       do k=1,nlev_toga
2628        plev_prof(k) = 100.*(plev_toga(k,it_toga2)                         &
2629     &          -frac*(plev_toga(k,it_toga2)-plev_toga(k,it_toga1)))
2630        t_prof(k) = t_toga(k,it_toga2)                                     &
2631     &          -frac*(t_toga(k,it_toga2)-t_toga(k,it_toga1))
2632        q_prof(k) = q_toga(k,it_toga2)                                     &
2633     &          -frac*(q_toga(k,it_toga2)-q_toga(k,it_toga1))
2634        u_prof(k) = u_toga(k,it_toga2)                                     &
2635     &          -frac*(u_toga(k,it_toga2)-u_toga(k,it_toga1))
2636        v_prof(k) = v_toga(k,it_toga2)                                     &
2637     &          -frac*(v_toga(k,it_toga2)-v_toga(k,it_toga1))
2638        w_prof(k) = w_toga(k,it_toga2)                                     &
2639     &          -frac*(w_toga(k,it_toga2)-w_toga(k,it_toga1))
2640        ht_prof(k) = ht_toga(k,it_toga2)                                   &
2641     &          -frac*(ht_toga(k,it_toga2)-ht_toga(k,it_toga1))
2642        vt_prof(k) = vt_toga(k,it_toga2)                                   &
2643     &          -frac*(vt_toga(k,it_toga2)-vt_toga(k,it_toga1))
2644        hq_prof(k) = hq_toga(k,it_toga2)                                   &
2645     &          -frac*(hq_toga(k,it_toga2)-hq_toga(k,it_toga1))
2646        vq_prof(k) = vq_toga(k,it_toga2)                                   &
2647     &          -frac*(vq_toga(k,it_toga2)-vq_toga(k,it_toga1))
2648        enddo
2649
2650        return
2651        END
2652
2653!======================================================================
2654        SUBROUTINE interp_armcu_time(day,day1,annee_ref                    &
2655     &             ,year_ini_armcu,day_ini_armcu,nt_armcu,dt_armcu         &
2656     &             ,nlev_armcu,fs_armcu,fl_armcu,at_armcu,rt_armcu         &
2657     &             ,aqt_armcu,fs_prof,fl_prof,at_prof,rt_prof,aqt_prof)
2658        implicit none
2659
2660!---------------------------------------------------------------------------------------
2661! Time interpolation of a 2D field to the timestep corresponding to day
2662!
2663! day: current julian day (e.g. 717538.2)
2664! day1: first day of the simulation
2665! nt_armcu: total nb of data in the forcing (e.g. 31 for armcu)
2666! dt_armcu: total time interval (in sec) between 2 forcing data (e.g. 1/2h for armcu)
2667! fs= sensible flux
2668! fl= latent flux
2669! at,rt,aqt= advective and radiative tendencies
2670!---------------------------------------------------------------------------------------
2671
2672! inputs:
2673        integer annee_ref
2674        integer nt_armcu,nlev_armcu
2675        integer year_ini_armcu
2676        real day, day1,day_ini_armcu,dt_armcu
2677        real fs_armcu(nt_armcu),fl_armcu(nt_armcu),at_armcu(nt_armcu)
2678        real rt_armcu(nt_armcu),aqt_armcu(nt_armcu)
2679! outputs:
2680        real fs_prof,fl_prof,at_prof,rt_prof,aqt_prof
2681! local:
2682        integer it_armcu1, it_armcu2,k
2683        real timeit,time_armcu1,time_armcu2,frac
2684
2685! Check that initial day of the simulation consistent with ARMCU period:
2686       if (annee_ref.ne.1997 ) then
2687        print*,'Pour ARMCU, annee_ref doit etre 1997 '
2688        print*,'Changer annee_ref dans run.def'
2689        stop
2690       endif
2691
2692      timeit=(day-day_ini_armcu)*86400
2693
2694! Determine the closest observation times:
2695       it_armcu1=INT(timeit/dt_armcu)+1
2696       it_armcu2=it_armcu1 + 1
2697       time_armcu1=(it_armcu1-1)*dt_armcu
2698       time_armcu2=(it_armcu2-1)*dt_armcu
2699       print *,'timeit day day_ini_armcu',timeit,day,day_ini_armcu
2700       print *,'it_armcu1,it_armcu2,time_armcu1,time_armcu2',              &
2701     &          it_armcu1,it_armcu2,time_armcu1,time_armcu2
2702
2703       if (it_armcu1 .ge. nt_armcu) then
2704        write(*,*) 'PB-stop: day, it_armcu1, it_armcu2, timeit: '          &
2705     &        ,day,it_armcu1,it_armcu2,timeit/86400.
2706        stop
2707       endif
2708
2709! time interpolation:
2710       frac=(time_armcu2-timeit)/(time_armcu2-time_armcu1)
2711       frac=max(frac,0.0)
2712
2713       fs_prof = fs_armcu(it_armcu2)                                       &
2714     &          -frac*(fs_armcu(it_armcu2)-fs_armcu(it_armcu1))
2715       fl_prof = fl_armcu(it_armcu2)                                       &
2716     &          -frac*(fl_armcu(it_armcu2)-fl_armcu(it_armcu1))
2717       at_prof = at_armcu(it_armcu2)                                       &
2718     &          -frac*(at_armcu(it_armcu2)-at_armcu(it_armcu1))
2719       rt_prof = rt_armcu(it_armcu2)                                       &
2720     &          -frac*(rt_armcu(it_armcu2)-rt_armcu(it_armcu1))
2721       aqt_prof = aqt_armcu(it_armcu2)                                       &
2722     &          -frac*(aqt_armcu(it_armcu2)-aqt_armcu(it_armcu1))
2723
2724         print*,                                                           &
2725     &'day,annee_ref,day_ini_armcu,timeit,it_armcu1,it_armcu2,SST:',       &
2726     &day,annee_ref,day_ini_armcu,timeit/86400.,it_armcu1,                 &
2727     &it_armcu2,fs_prof,fl_prof,at_prof,rt_prof,aqt_prof
2728
2729        return
2730        END
2731
2732!=====================================================================
2733      subroutine readprofiles(nlev_max,kmax,ntrac,height,                  &
2734     &           thlprof,qtprof,uprof,                                     &
2735     &           vprof,e12prof,ugprof,vgprof,                              &
2736     &           wfls,dqtdxls,dqtdyls,dqtdtls,                             &
2737     &           thlpcar,tracer,nt1,nt2)
2738      implicit none
2739
2740        integer nlev_max,kmax,kmax2,ntrac
2741        logical :: llesread = .true.
2742
2743        real height(nlev_max),thlprof(nlev_max),qtprof(nlev_max),          &
2744     &       uprof(nlev_max),vprof(nlev_max),e12prof(nlev_max),            &
2745     &       ugprof(nlev_max),vgprof(nlev_max),wfls(nlev_max),             &
2746     &       dqtdxls(nlev_max),dqtdyls(nlev_max),dqtdtls(nlev_max),        &
2747     &           thlpcar(nlev_max),tracer(nlev_max,ntrac)
2748
2749        integer, parameter :: ilesfile=1
2750        integer :: ierr,k,itrac,nt1,nt2
2751
2752        if(.not.(llesread)) return
2753
2754       open (ilesfile,file='prof.inp.001',status='old',iostat=ierr)
2755        if (ierr /= 0) stop 'ERROR:Prof.inp does not exist'
2756        read (ilesfile,*) kmax
2757        do k=1,kmax
2758          read (ilesfile,*) height(k),thlprof(k),qtprof (k),               &
2759     &                      uprof (k),vprof  (k),e12prof(k)
2760        enddo
2761        close(ilesfile)
2762
2763       open(ilesfile,file='lscale.inp.001',status='old',iostat=ierr)
2764        if (ierr /= 0) stop 'ERROR:Lscale.inp does not exist'
2765        read (ilesfile,*) kmax2
2766        if (kmax .ne. kmax2) then
2767          print *, 'fichiers prof.inp et lscale.inp incompatibles :'
2768          print *, 'nbre de niveaux : ',kmax,' et ',kmax2
2769          stop 'lecture profiles'
2770        endif
2771        do k=1,kmax
2772          read (ilesfile,*) height(k),ugprof(k),vgprof(k),wfls(k),         &
2773     &                      dqtdxls(k),dqtdyls(k),dqtdtls(k),thlpcar(k)
2774        end do
2775        close(ilesfile)
2776
2777       open(ilesfile,file='trac.inp.001',status='old',iostat=ierr)
2778        if (ierr /= 0) then
2779            print*,'WARNING : trac.inp does not exist'
2780        else
2781        read (ilesfile,*) kmax2,nt1,nt2
2782        if (nt2>ntrac) then
2783          stop'Augmenter le nombre de traceurs dans traceur.def'
2784        endif
2785        if (kmax .ne. kmax2) then
2786          print *, 'fichiers prof.inp et lscale.inp incompatibles :'
2787          print *, 'nbre de niveaux : ',kmax,' et ',kmax2
2788          stop 'lecture profiles'
2789        endif
2790        do k=1,kmax
2791          read (ilesfile,*) height(k),(tracer(k,itrac),itrac=nt1,nt2)
2792        end do
2793        close(ilesfile)
2794        endif
2795
2796        return
2797        end
2798!======================================================================
2799      subroutine readprofile_sandu(nlev_max,kmax,height,pprof,tprof,       &
2800     &       thlprof,qprof,uprof,vprof,wprof,omega,o3mmr)
2801!======================================================================
2802      implicit none
2803
2804        integer nlev_max,kmax
2805        logical :: llesread = .true.
2806
2807        real height(nlev_max),pprof(nlev_max),tprof(nlev_max)
2808        real thlprof(nlev_max)
2809        real qprof(nlev_max),uprof(nlev_max),vprof(nlev_max)
2810        real wprof(nlev_max),omega(nlev_max),o3mmr(nlev_max)
2811
2812        integer, parameter :: ilesfile=1
2813        integer :: k,ierr
2814
2815        if(.not.(llesread)) return
2816
2817       open (ilesfile,file='prof.inp.001',status='old',iostat=ierr)
2818        if (ierr /= 0) stop 'ERROR:Prof.inp does not exist'
2819        read (ilesfile,*) kmax
2820        do k=1,kmax
2821          read (ilesfile,*) height(k),pprof(k),  tprof(k),thlprof(k),      &
2822     &                      qprof (k),uprof(k),  vprof(k),  wprof(k),      &
2823     &                      omega (k),o3mmr(k)
2824        enddo
2825        close(ilesfile)
2826
2827        return
2828        end
2829
2830!======================================================================
2831      subroutine readprofile_astex(nlev_max,kmax,height,pprof,tprof,       &
2832     &    thlprof,qvprof,qlprof,qtprof,uprof,vprof,wprof,tkeprof,o3mmr)
2833!======================================================================
2834      implicit none
2835
2836        integer nlev_max,kmax
2837        logical :: llesread = .true.
2838
2839        real height(nlev_max),pprof(nlev_max),tprof(nlev_max),             &
2840     &  thlprof(nlev_max),qlprof(nlev_max),qtprof(nlev_max),               &
2841     &  qvprof(nlev_max),uprof(nlev_max),vprof(nlev_max),                  &
2842     &  wprof(nlev_max),tkeprof(nlev_max),o3mmr(nlev_max)
2843
2844        integer, parameter :: ilesfile=1
2845        integer :: ierr,k
2846
2847        if(.not.(llesread)) return
2848
2849       open (ilesfile,file='prof.inp.001',status='old',iostat=ierr)
2850        if (ierr /= 0) stop 'ERROR:Prof.inp does not exist'
2851        read (ilesfile,*) kmax
2852        do k=1,kmax
2853          read (ilesfile,*) height(k),pprof(k),  tprof(k),thlprof(k),      &
2854     &                qvprof (k),qlprof (k),qtprof (k),                    &
2855     &                uprof(k),  vprof(k),  wprof(k),tkeprof(k),o3mmr(k)
2856        enddo
2857        close(ilesfile)
2858
2859        return
2860        end
2861
2862
2863
2864!======================================================================
2865      subroutine readprofile_armcu(nlev_max,kmax,height,pprof,uprof,       &
2866     &       vprof,thetaprof,tprof,qvprof,rvprof,aprof,bprof)
2867!======================================================================
2868      implicit none
2869
2870        integer nlev_max,kmax
2871        logical :: llesread = .true.
2872
2873        real height(nlev_max),pprof(nlev_max),tprof(nlev_max)
2874        real thetaprof(nlev_max),rvprof(nlev_max)
2875        real qvprof(nlev_max),uprof(nlev_max),vprof(nlev_max)
2876        real aprof(nlev_max+1),bprof(nlev_max+1)
2877
2878        integer, parameter :: ilesfile=1
2879        integer, parameter :: ifile=2
2880        integer :: ierr,jtot,k
2881
2882        if(.not.(llesread)) return
2883
2884! Read profiles at full levels
2885       IF(nlev_max.EQ.19) THEN
2886       open (ilesfile,file='prof.inp.19',status='old',iostat=ierr)
2887       print *,'On ouvre prof.inp.19'
2888       ELSE
2889       open (ilesfile,file='prof.inp.40',status='old',iostat=ierr)
2890       print *,'On ouvre prof.inp.40'
2891       ENDIF
2892        if (ierr /= 0) stop 'ERROR:Prof.inp does not exist'
2893        read (ilesfile,*) kmax
2894        do k=1,kmax
2895          read (ilesfile,*) height(k)    ,pprof(k),  uprof(k), vprof(k),   &
2896     &                      thetaprof(k) ,tprof(k), qvprof(k),rvprof(k)
2897        enddo
2898        close(ilesfile)
2899
2900! Vertical coordinates half levels for eta-coordinates (plev = alpha + beta * psurf) 
2901       IF(nlev_max.EQ.19) THEN
2902       open (ifile,file='proh.inp.19',status='old',iostat=ierr)
2903       print *,'On ouvre proh.inp.19'
2904       if (ierr /= 0) stop 'ERROR:Proh.inp.19 does not exist'
2905       ELSE
2906       open (ifile,file='proh.inp.40',status='old',iostat=ierr)
2907       print *,'On ouvre proh.inp.40'
2908       if (ierr /= 0) stop 'ERROR:Proh.inp.40 does not exist'
2909       ENDIF
2910        read (ifile,*) kmax
2911        do k=1,kmax
2912          read (ifile,*) jtot,aprof(k),bprof(k)
2913        enddo
2914        close(ifile)
2915
2916        return
2917        end
2918!=====================================================================
2919      subroutine read_amma(fich_amma,nlevel,ntime                          &
2920     &     ,zz,pp,temp,qv,u,v,dw                                           &
2921     &     ,dt,dq,sens,flat)
2922
2923!program reading forcings of the AMMA case study
2924
2925
2926      implicit none
2927
2928#include "netcdf.inc"
2929
2930      integer ntime,nlevel
2931      character*80 :: fich_amma
2932      real*8 zz(nlevel)
2933
2934      real*8 temp(nlevel),pp(nlevel)
2935      real*8 qv(nlevel),u(nlevel)
2936      real*8 v(nlevel)
2937      real*8 dw(nlevel,ntime)
2938      real*8 dt(nlevel,ntime)
2939      real*8 dq(nlevel,ntime)
2940      real*8 flat(ntime),sens(ntime)
2941
2942      integer nid, ierr
2943      integer nbvar3d
2944      parameter(nbvar3d=30)
2945      integer var3didin(nbvar3d)
2946
2947      ierr = NF_OPEN(fich_amma,NF_NOWRITE,nid)
2948      if (ierr.NE.NF_NOERR) then
2949         write(*,*) 'ERROR: Pb opening forcings nc file '
2950         write(*,*) NF_STRERROR(ierr)
2951         stop ""
2952      endif
2953
2954
2955       ierr=NF_INQ_VARID(nid,"zz",var3didin(1)) 
2956         if(ierr/=NF_NOERR) then
2957           write(*,*) NF_STRERROR(ierr)
2958           stop 'lev'
2959         endif
2960
2961
2962      ierr=NF_INQ_VARID(nid,"temp",var3didin(2))
2963         if(ierr/=NF_NOERR) then
2964           write(*,*) NF_STRERROR(ierr)
2965           stop 'temp'
2966         endif
2967
2968      ierr=NF_INQ_VARID(nid,"qv",var3didin(3))
2969         if(ierr/=NF_NOERR) then
2970           write(*,*) NF_STRERROR(ierr)
2971           stop 'qv'
2972         endif
2973
2974      ierr=NF_INQ_VARID(nid,"u",var3didin(4))
2975         if(ierr/=NF_NOERR) then
2976           write(*,*) NF_STRERROR(ierr)
2977           stop 'u'
2978         endif
2979
2980      ierr=NF_INQ_VARID(nid,"v",var3didin(5))
2981         if(ierr/=NF_NOERR) then
2982           write(*,*) NF_STRERROR(ierr)
2983           stop 'v'
2984         endif
2985
2986      ierr=NF_INQ_VARID(nid,"dw",var3didin(6))
2987         if(ierr/=NF_NOERR) then
2988           write(*,*) NF_STRERROR(ierr)
2989           stop 'dw'
2990         endif
2991
2992      ierr=NF_INQ_VARID(nid,"dt",var3didin(7))
2993         if(ierr/=NF_NOERR) then
2994           write(*,*) NF_STRERROR(ierr)
2995           stop 'dt'
2996         endif
2997
2998      ierr=NF_INQ_VARID(nid,"dq",var3didin(8))
2999         if(ierr/=NF_NOERR) then
3000           write(*,*) NF_STRERROR(ierr)
3001           stop 'dq'
3002         endif
3003     
3004      ierr=NF_INQ_VARID(nid,"sens",var3didin(9))
3005         if(ierr/=NF_NOERR) then
3006           write(*,*) NF_STRERROR(ierr)
3007           stop 'sens'
3008         endif
3009
3010      ierr=NF_INQ_VARID(nid,"flat",var3didin(10))
3011         if(ierr/=NF_NOERR) then
3012           write(*,*) NF_STRERROR(ierr)
3013           stop 'flat'
3014         endif
3015
3016      ierr=NF_INQ_VARID(nid,"pp",var3didin(11))
3017         if(ierr/=NF_NOERR) then
3018           write(*,*) NF_STRERROR(ierr)
3019           stop 'pp'
3020      endif
3021
3022!dimensions lecture
3023!      call catchaxis(nid,ntime,nlevel,time,z,ierr)
3024 
3025#ifdef NC_DOUBLE
3026         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(1),zz)
3027#else
3028         ierr = NF_GET_VAR_REAL(nid,var3didin(1),zz)
3029#endif
3030         if(ierr/=NF_NOERR) then
3031            write(*,*) NF_STRERROR(ierr)
3032            stop "getvarup"
3033         endif
3034!          write(*,*)'lecture z ok',zz
3035
3036#ifdef NC_DOUBLE
3037         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(2),temp)
3038#else
3039         ierr = NF_GET_VAR_REAL(nid,var3didin(2),temp)
3040#endif
3041         if(ierr/=NF_NOERR) then
3042            write(*,*) NF_STRERROR(ierr)
3043            stop "getvarup"
3044         endif
3045!          write(*,*)'lecture th ok',temp
3046
3047#ifdef NC_DOUBLE
3048         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(3),qv)
3049#else
3050         ierr = NF_GET_VAR_REAL(nid,var3didin(3),qv)
3051#endif
3052         if(ierr/=NF_NOERR) then
3053            write(*,*) NF_STRERROR(ierr)
3054            stop "getvarup"
3055         endif
3056!          write(*,*)'lecture qv ok',qv
3057 
3058#ifdef NC_DOUBLE
3059         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(4),u)
3060#else
3061         ierr = NF_GET_VAR_REAL(nid,var3didin(4),u)
3062#endif
3063         if(ierr/=NF_NOERR) then
3064            write(*,*) NF_STRERROR(ierr)
3065            stop "getvarup"
3066         endif
3067!          write(*,*)'lecture u ok',u
3068
3069#ifdef NC_DOUBLE
3070         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(5),v)
3071#else
3072         ierr = NF_GET_VAR_REAL(nid,var3didin(5),v)
3073#endif
3074         if(ierr/=NF_NOERR) then
3075            write(*,*) NF_STRERROR(ierr)
3076            stop "getvarup"
3077         endif
3078!          write(*,*)'lecture v ok',v
3079
3080#ifdef NC_DOUBLE
3081         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(6),dw)
3082#else
3083         ierr = NF_GET_VAR_REAL(nid,var3didin(6),dw)
3084#endif
3085         if(ierr/=NF_NOERR) then
3086            write(*,*) NF_STRERROR(ierr)
3087            stop "getvarup"
3088         endif
3089!          write(*,*)'lecture w ok',dw
3090
3091#ifdef NC_DOUBLE
3092         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(7),dt)
3093#else
3094         ierr = NF_GET_VAR_REAL(nid,var3didin(7),dt)
3095#endif
3096         if(ierr/=NF_NOERR) then
3097            write(*,*) NF_STRERROR(ierr)
3098            stop "getvarup"
3099         endif
3100!          write(*,*)'lecture dt ok',dt
3101
3102#ifdef NC_DOUBLE
3103         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(8),dq)
3104#else
3105         ierr = NF_GET_VAR_REAL(nid,var3didin(8),dq)
3106#endif
3107         if(ierr/=NF_NOERR) then
3108            write(*,*) NF_STRERROR(ierr)
3109            stop "getvarup"
3110         endif
3111!          write(*,*)'lecture dq ok',dq
3112
3113#ifdef NC_DOUBLE
3114         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(9),sens)
3115#else
3116         ierr = NF_GET_VAR_REAL(nid,var3didin(9),sens)
3117#endif
3118         if(ierr/=NF_NOERR) then
3119            write(*,*) NF_STRERROR(ierr)
3120            stop "getvarup"
3121         endif
3122!          write(*,*)'lecture sens ok',sens
3123
3124#ifdef NC_DOUBLE
3125         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(10),flat)
3126#else
3127         ierr = NF_GET_VAR_REAL(nid,var3didin(10),flat)
3128#endif
3129         if(ierr/=NF_NOERR) then
3130            write(*,*) NF_STRERROR(ierr)
3131            stop "getvarup"
3132         endif
3133!          write(*,*)'lecture flat ok',flat
3134
3135#ifdef NC_DOUBLE
3136         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(11),pp)
3137#else
3138         ierr = NF_GET_VAR_REAL(nid,var3didin(11),pp)
3139#endif
3140         if(ierr/=NF_NOERR) then
3141            write(*,*) NF_STRERROR(ierr)
3142            stop "getvarup"
3143         endif
3144!          write(*,*)'lecture pp ok',pp
3145
3146         return 
3147         end subroutine read_amma
3148!======================================================================
3149        SUBROUTINE interp_amma_time(day,day1,annee_ref                     &
3150     &         ,year_ini_amma,day_ini_amma,nt_amma,dt_amma,nlev_amma       &
3151     &         ,vitw_amma,ht_amma,hq_amma,lat_amma,sens_amma               &
3152     &         ,vitw_prof,ht_prof,hq_prof,lat_prof,sens_prof)
3153        implicit none
3154
3155!---------------------------------------------------------------------------------------
3156! Time interpolation of a 2D field to the timestep corresponding to day
3157!
3158! day: current julian day (e.g. 717538.2)
3159! day1: first day of the simulation
3160! nt_amma: total nb of data in the forcing (e.g. 48 for AMMA)
3161! dt_amma: total time interval (in sec) between 2 forcing data (e.g. 30min for AMMA)
3162!---------------------------------------------------------------------------------------
3163
3164#include "compar1d.h"
3165
3166! inputs:
3167        integer annee_ref
3168        integer nt_amma,nlev_amma
3169        integer year_ini_amma
3170        real day, day1,day_ini_amma,dt_amma
3171        real vitw_amma(nlev_amma,nt_amma)
3172        real ht_amma(nlev_amma,nt_amma)
3173        real hq_amma(nlev_amma,nt_amma)
3174        real lat_amma(nt_amma)
3175        real sens_amma(nt_amma)
3176! outputs:
3177        real vitw_prof(nlev_amma)
3178        real ht_prof(nlev_amma)
3179        real hq_prof(nlev_amma)
3180        real lat_prof,sens_prof
3181! local:
3182        integer it_amma1, it_amma2,k
3183        real timeit,time_amma1,time_amma2,frac
3184
3185
3186        if (forcing_type.eq.6) then
3187! Check that initial day of the simulation consistent with AMMA case:
3188       if (annee_ref.ne.2006) then
3189        print*,'Pour AMMA, annee_ref doit etre 2006'
3190        print*,'Changer annee_ref dans run.def'
3191        stop
3192       endif
3193       if (annee_ref.eq.2006 .and. day1.lt.day_ini_amma) then
3194        print*,'AMMA a débuté le 10 juillet 2006',day1,day_ini_amma
3195        print*,'Changer dayref dans run.def'
3196        stop
3197       endif
3198       if (annee_ref.eq.2006 .and. day1.gt.day_ini_amma+1) then
3199        print*,'AMMA a fini le 11 juillet'
3200        print*,'Changer dayref ou nday dans run.def'
3201        stop
3202       endif
3203       endif
3204
3205! Determine timestep relative to the 1st day of AMMA:
3206!       timeit=(day-day1)*86400.
3207!       if (annee_ref.eq.1992) then
3208!        timeit=(day-day_ini_toga)*86400.
3209!       else
3210!        timeit=(day+61.-1.)*86400. ! 61 days between Nov01 and Dec31 1992
3211!       endif
3212      timeit=(day-day_ini_amma)*86400
3213
3214! Determine the closest observation times:
3215!       it_amma1=INT(timeit/dt_amma)+1
3216!       it_amma2=it_amma1 + 1
3217!       time_amma1=(it_amma1-1)*dt_amma
3218!       time_amma2=(it_amma2-1)*dt_amma
3219
3220       it_amma1=INT(timeit/dt_amma)+1
3221       IF (it_amma1 .EQ. nt_amma) THEN
3222       it_amma2=it_amma1
3223       ELSE
3224       it_amma2=it_amma1 + 1
3225       ENDIF
3226       time_amma1=(it_amma1-1)*dt_amma
3227       time_amma2=(it_amma2-1)*dt_amma
3228
3229       if (it_amma1 .gt. nt_amma) then
3230        write(*,*) 'PB-stop: day, it_amma1, it_amma2, timeit: '            &
3231     &        ,day,day_ini_amma,it_amma1,it_amma2,timeit/86400.
3232        stop
3233       endif
3234
3235! time interpolation:
3236       frac=(time_amma2-timeit)/(time_amma2-time_amma1)
3237       frac=max(frac,0.0)
3238
3239       lat_prof = lat_amma(it_amma2)                                       &
3240     &          -frac*(lat_amma(it_amma2)-lat_amma(it_amma1)) 
3241       sens_prof = sens_amma(it_amma2)                                     &
3242     &          -frac*(sens_amma(it_amma2)-sens_amma(it_amma1))
3243
3244       do k=1,nlev_amma
3245        vitw_prof(k) = vitw_amma(k,it_amma2)                               &
3246     &          -frac*(vitw_amma(k,it_amma2)-vitw_amma(k,it_amma1))
3247        ht_prof(k) = ht_amma(k,it_amma2)                                   &
3248     &          -frac*(ht_amma(k,it_amma2)-ht_amma(k,it_amma1))
3249        hq_prof(k) = hq_amma(k,it_amma2)                                   &
3250     &          -frac*(hq_amma(k,it_amma2)-hq_amma(k,it_amma1))
3251        enddo
3252
3253        return
3254        END
3255
3256!=====================================================================
3257      subroutine read_fire(fich_fire,nlevel,ntime                          &
3258     &     ,zz,thl,qt,u,v,tke                                              &
3259     &     ,ug,vg,wls,dqtdx,dqtdy,dqtdt,thl_rad)
3260
3261!program reading forcings of the FIRE case study
3262
3263
3264      implicit none
3265
3266#include "netcdf.inc"
3267
3268      integer ntime,nlevel
3269      character*80 :: fich_fire
3270      real*8 zz(nlevel)
3271
3272      real*8 thl(nlevel)
3273      real*8 qt(nlevel),u(nlevel)
3274      real*8 v(nlevel),tke(nlevel)
3275      real*8 ug(nlevel,ntime),vg(nlevel,ntime),wls(nlevel,ntime)
3276      real*8 dqtdx(nlevel,ntime),dqtdy(nlevel,ntime)
3277      real*8 dqtdt(nlevel,ntime),thl_rad(nlevel,ntime)
3278
3279      integer nid, ierr
3280      integer nbvar3d
3281      parameter(nbvar3d=30)
3282      integer var3didin(nbvar3d)
3283
3284      ierr = NF_OPEN(fich_fire,NF_NOWRITE,nid)
3285      if (ierr.NE.NF_NOERR) then
3286         write(*,*) 'ERROR: Pb opening forcings nc file '
3287         write(*,*) NF_STRERROR(ierr)
3288         stop ""
3289      endif
3290
3291
3292       ierr=NF_INQ_VARID(nid,"zz",var3didin(1)) 
3293         if(ierr/=NF_NOERR) then
3294           write(*,*) NF_STRERROR(ierr)
3295           stop 'lev'
3296         endif
3297
3298
3299      ierr=NF_INQ_VARID(nid,"thetal",var3didin(2))
3300         if(ierr/=NF_NOERR) then
3301           write(*,*) NF_STRERROR(ierr)
3302           stop 'temp'
3303         endif
3304
3305      ierr=NF_INQ_VARID(nid,"qt",var3didin(3))
3306         if(ierr/=NF_NOERR) then
3307           write(*,*) NF_STRERROR(ierr)
3308           stop 'qv'
3309         endif
3310
3311      ierr=NF_INQ_VARID(nid,"u",var3didin(4))
3312         if(ierr/=NF_NOERR) then
3313           write(*,*) NF_STRERROR(ierr)
3314           stop 'u'
3315         endif
3316
3317      ierr=NF_INQ_VARID(nid,"v",var3didin(5))
3318         if(ierr/=NF_NOERR) then
3319           write(*,*) NF_STRERROR(ierr)
3320           stop 'v'
3321         endif
3322
3323      ierr=NF_INQ_VARID(nid,"tke",var3didin(6))
3324         if(ierr/=NF_NOERR) then
3325           write(*,*) NF_STRERROR(ierr)
3326           stop 'tke'
3327         endif
3328
3329      ierr=NF_INQ_VARID(nid,"ugeo",var3didin(7))
3330         if(ierr/=NF_NOERR) then
3331           write(*,*) NF_STRERROR(ierr)
3332           stop 'ug'
3333         endif
3334
3335      ierr=NF_INQ_VARID(nid,"vgeo",var3didin(8))
3336         if(ierr/=NF_NOERR) then
3337           write(*,*) NF_STRERROR(ierr)
3338           stop 'vg'
3339         endif
3340     
3341      ierr=NF_INQ_VARID(nid,"wls",var3didin(9))
3342         if(ierr/=NF_NOERR) then
3343           write(*,*) NF_STRERROR(ierr)
3344           stop 'wls'
3345         endif
3346
3347      ierr=NF_INQ_VARID(nid,"dqtdx",var3didin(10))
3348         if(ierr/=NF_NOERR) then
3349           write(*,*) NF_STRERROR(ierr)
3350           stop 'dqtdx'
3351         endif
3352
3353      ierr=NF_INQ_VARID(nid,"dqtdy",var3didin(11))
3354         if(ierr/=NF_NOERR) then
3355           write(*,*) NF_STRERROR(ierr)
3356           stop 'dqtdy'
3357      endif
3358
3359      ierr=NF_INQ_VARID(nid,"dqtdt",var3didin(12))
3360         if(ierr/=NF_NOERR) then
3361           write(*,*) NF_STRERROR(ierr)
3362           stop 'dqtdt'
3363      endif
3364
3365      ierr=NF_INQ_VARID(nid,"thl_rad",var3didin(13))
3366         if(ierr/=NF_NOERR) then
3367           write(*,*) NF_STRERROR(ierr)
3368           stop 'thl_rad'
3369      endif
3370!dimensions lecture
3371!      call catchaxis(nid,ntime,nlevel,time,z,ierr)
3372 
3373#ifdef NC_DOUBLE
3374         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(1),zz)
3375#else
3376         ierr = NF_GET_VAR_REAL(nid,var3didin(1),zz)
3377#endif
3378         if(ierr/=NF_NOERR) then
3379            write(*,*) NF_STRERROR(ierr)
3380            stop "getvarup"
3381         endif
3382!          write(*,*)'lecture z ok',zz
3383
3384#ifdef NC_DOUBLE
3385         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(2),thl)
3386#else
3387         ierr = NF_GET_VAR_REAL(nid,var3didin(2),thl)
3388#endif
3389         if(ierr/=NF_NOERR) then
3390            write(*,*) NF_STRERROR(ierr)
3391            stop "getvarup"
3392         endif
3393!          write(*,*)'lecture thl ok',thl
3394
3395#ifdef NC_DOUBLE
3396         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(3),qt)
3397#else
3398         ierr = NF_GET_VAR_REAL(nid,var3didin(3),qt)
3399#endif
3400         if(ierr/=NF_NOERR) then
3401            write(*,*) NF_STRERROR(ierr)
3402            stop "getvarup"
3403         endif
3404!          write(*,*)'lecture qt ok',qt
3405 
3406#ifdef NC_DOUBLE
3407         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(4),u)
3408#else
3409         ierr = NF_GET_VAR_REAL(nid,var3didin(4),u)
3410#endif
3411         if(ierr/=NF_NOERR) then
3412            write(*,*) NF_STRERROR(ierr)
3413            stop "getvarup"
3414         endif
3415!          write(*,*)'lecture u ok',u
3416
3417#ifdef NC_DOUBLE
3418         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(5),v)
3419#else
3420         ierr = NF_GET_VAR_REAL(nid,var3didin(5),v)
3421#endif
3422         if(ierr/=NF_NOERR) then
3423            write(*,*) NF_STRERROR(ierr)
3424            stop "getvarup"
3425         endif
3426!          write(*,*)'lecture v ok',v
3427
3428#ifdef NC_DOUBLE
3429         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(6),tke)
3430#else
3431         ierr = NF_GET_VAR_REAL(nid,var3didin(6),tke)
3432#endif
3433         if(ierr/=NF_NOERR) then
3434            write(*,*) NF_STRERROR(ierr)
3435            stop "getvarup"
3436         endif
3437!          write(*,*)'lecture tke ok',tke
3438
3439#ifdef NC_DOUBLE
3440         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(7),ug)
3441#else
3442         ierr = NF_GET_VAR_REAL(nid,var3didin(7),ug)
3443#endif
3444         if(ierr/=NF_NOERR) then
3445            write(*,*) NF_STRERROR(ierr)
3446            stop "getvarup"
3447         endif
3448!          write(*,*)'lecture ug ok',ug
3449
3450#ifdef NC_DOUBLE
3451         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(8),vg)
3452#else
3453         ierr = NF_GET_VAR_REAL(nid,var3didin(8),vg)
3454#endif
3455         if(ierr/=NF_NOERR) then
3456            write(*,*) NF_STRERROR(ierr)
3457            stop "getvarup"
3458         endif
3459!          write(*,*)'lecture vg ok',vg
3460
3461#ifdef NC_DOUBLE
3462         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(9),wls)
3463#else
3464         ierr = NF_GET_VAR_REAL(nid,var3didin(9),wls)
3465#endif
3466         if(ierr/=NF_NOERR) then
3467            write(*,*) NF_STRERROR(ierr)
3468            stop "getvarup"
3469         endif
3470!          write(*,*)'lecture wls ok',wls
3471
3472#ifdef NC_DOUBLE
3473         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(10),dqtdx)
3474#else
3475         ierr = NF_GET_VAR_REAL(nid,var3didin(10),dqtdx)
3476#endif
3477         if(ierr/=NF_NOERR) then
3478            write(*,*) NF_STRERROR(ierr)
3479            stop "getvarup"
3480         endif
3481!          write(*,*)'lecture dqtdx ok',dqtdx
3482
3483#ifdef NC_DOUBLE
3484         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(11),dqtdy)
3485#else
3486         ierr = NF_GET_VAR_REAL(nid,var3didin(11),dqtdy)
3487#endif
3488         if(ierr/=NF_NOERR) then
3489            write(*,*) NF_STRERROR(ierr)
3490            stop "getvarup"
3491         endif
3492!          write(*,*)'lecture dqtdy ok',dqtdy
3493
3494#ifdef NC_DOUBLE
3495         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(12),dqtdt)
3496#else
3497         ierr = NF_GET_VAR_REAL(nid,var3didin(12),dqtdt)
3498#endif
3499         if(ierr/=NF_NOERR) then
3500            write(*,*) NF_STRERROR(ierr)
3501            stop "getvarup"
3502         endif
3503!          write(*,*)'lecture dqtdt ok',dqtdt
3504
3505#ifdef NC_DOUBLE
3506         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(13),thl_rad)
3507#else
3508         ierr = NF_GET_VAR_REAL(nid,var3didin(13),thl_rad)
3509#endif
3510         if(ierr/=NF_NOERR) then
3511            write(*,*) NF_STRERROR(ierr)
3512            stop "getvarup"
3513         endif
3514!          write(*,*)'lecture thl_rad ok',thl_rad
3515
3516         return 
3517         end subroutine read_fire
Note: See TracBrowser for help on using the repository browser.