source: LMDZ5/trunk/libf/phylmd/dyn1d/1DUTILS.h @ 2602

Last change on this file since 2602 was 2601, checked in by Ehouarn Millour, 8 years ago

Cleanup in the dynamics: turn temps.h into module temps_mod.F90
EM

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