source: LMDZ5/trunk/libf/phy1d/1DUTILS.h @ 1914

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

Modifs pour le cas FIRE
MPL


FIRE case mods
MPL

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