source: LMDZ5/branches/testing/libf/phy1d/1DUTILS.h_no_writelim @ 1856

Last change on this file since 1856 was 1795, checked in by Ehouarn Millour, 11 years ago

Version testing basee sur la r1794


Testing release based on r1794

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