LMDZ
1DUTILS.h
Go to the documentation of this file.
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 
438  IMPLICIT NONE
439 !=======================================================
440 ! Ecriture du fichier de redemarrage sous format NetCDF
441 !=======================================================
442 ! Declarations:
443 ! -------------
444  include "dimensions.h"
445  include "comconst.h"
446  include "temps.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 
582  IMPLICIT NONE
583 !=======================================================
584 ! Ecriture du fichier de redemarrage sous format NetCDF
585 !=======================================================
586 ! Declarations:
587 ! -------------
588  include "dimensions.h"
589  include "comconst.h"
590  include "temps.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))
956 101 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(k )- play(k-1 ))
1215 ! dudp (k-1) = (u (k )- u (k-1 ))/dph(k-1)
1216 ! dvdp (k-1) = (v (k )- 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(k )- 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)' print*,'Changer dayref dans run.def' stop endif if (annee_ref.eq.1993 .and. day1.gt.day_ini_toga+119) then print*,'TOGA-COARE a fini le 28 Feb 1993 (jour julien=59)' print*,'Changer dayref ou nday dans run.def' stop endif else if (forcing_type.eq.4) then ! Check that initial day of the simulation consistent with TWP-ICE period: if (annee_ref.ne.2006) then print*,'Pour TWP-ICE, annee_ref doit etre 2006' print*,'Changer annee_ref dans run.def' stop endif if (annee_ref.eq.2006 .and. day1.lt.day_ini_toga) then print*,'TWP-ICE a debute le 17 Jan 2006 (jour julien=17)' print*,'Changer dayref dans run.def' stop endif if (annee_ref.eq.2006 .and. day1.gt.day_ini_toga+26) then print*,'TWP-ICE a fini le 12 Feb 2006 (jour julien=43)' print*,'Changer dayref ou nday dans run.def' stop endif endif ! Determine timestep relative to the 1st day of TOGA-COARE: ! timeit=(day-day1)*86400. ! if (annee_ref.eq.1992) then ! timeit=(day-day_ini_toga)*86400. ! else ! timeit=(day+61.-1.)*86400. ! 61 days between Nov01 and Dec31 1992 ! endif timeit=(day-day_ini_toga)*86400 ! Determine the closest observation times: it_toga1=INT(timeit/dt_toga)+1 it_toga2=it_toga1 + 1 time_toga1=(it_toga1-1)*dt_toga time_toga2=(it_toga2-1)*dt_toga if (it_toga1 .ge. nt_toga) then write(*,*) 'PB-stop: day, it_toga1, it_toga2, timeit: ' & & ,day,it_toga1,it_toga2,timeit/86400. stop endif ! time interpolation: frac=(time_toga2-timeit)/(time_toga2-time_toga1) frac=max(frac,0.0) ts_prof = ts_toga(it_toga2) & & -frac*(ts_toga(it_toga2)-ts_toga(it_toga1)) ! print*, ! :'day,annee_ref,day_ini_toga,timeit,it_toga1,it_toga2,SST:', ! :day,annee_ref,day_ini_toga,timeit/86400.,it_toga1,it_toga2,ts_prof do k=1,nlev_toga plev_prof(k) = 100.*(plev_toga(k,it_toga2) & & -frac*(plev_toga(k,it_toga2)-plev_toga(k,it_toga1))) t_prof(k) = t_toga(k,it_toga2) & & -frac*(t_toga(k,it_toga2)-t_toga(k,it_toga1)) q_prof(k) = q_toga(k,it_toga2) & & -frac*(q_toga(k,it_toga2)-q_toga(k,it_toga1)) u_prof(k) = u_toga(k,it_toga2) & & -frac*(u_toga(k,it_toga2)-u_toga(k,it_toga1)) v_prof(k) = v_toga(k,it_toga2) & & -frac*(v_toga(k,it_toga2)-v_toga(k,it_toga1)) w_prof(k) = w_toga(k,it_toga2) & & -frac*(w_toga(k,it_toga2)-w_toga(k,it_toga1)) ht_prof(k) = ht_toga(k,it_toga2) & & -frac*(ht_toga(k,it_toga2)-ht_toga(k,it_toga1)) vt_prof(k) = vt_toga(k,it_toga2) & & -frac*(vt_toga(k,it_toga2)-vt_toga(k,it_toga1)) hq_prof(k) = hq_toga(k,it_toga2) & & -frac*(hq_toga(k,it_toga2)-hq_toga(k,it_toga1)) vq_prof(k) = vq_toga(k,it_toga2) & & -frac*(vq_toga(k,it_toga2)-vq_toga(k,it_toga1)) enddo return END !====================================================================== SUBROUTINE interp_dice_time(day,day1,annee_ref & & ,year_ini_dice,day_ini_dice,nt_dice,dt_dice & & ,nlev_dice,shf_dice,lhf_dice,lwup_dice,swup_dice & & ,tg_dice,ustar_dice,psurf_dice,ug_dice,vg_dice & & ,ht_dice,hq_dice,hu_dice,hv_dice,w_dice,omega_dice & & ,shf_prof,lhf_prof,lwup_prof,swup_prof,tg_prof & & ,ustar_prof,psurf_prof,ug_prof,vg_prof & & ,ht_prof,hq_prof,hu_prof,hv_prof,w_prof,omega_prof) implicit none !--------------------------------------------------------------------------------------- ! Time interpolation of a 2D field to the timestep corresponding to day ! ! day: current julian day (e.g. 717538.2) ! day1: first day of the simulation ! nt_dice: total nb of data in the forcing (e.g. 145 for Dice) ! dt_dice: total time interval (in sec) between 2 forcing data (e.g. 30min. for Dice) !--------------------------------------------------------------------------------------- #include "compar1d.h" ! inputs: integer annee_ref integer nt_dice,nlev_dice integer year_ini_dice real day, day1,day_ini_dice,dt_dice real shf_dice(nt_dice),lhf_dice(nt_dice),lwup_dice(nt_dice) real swup_dice(nt_dice),tg_dice(nt_dice),ustar_dice(nt_dice) real psurf_dice(nt_dice),ug_dice(nt_dice),vg_dice(nt_dice) real ht_dice(nlev_dice,nt_dice),hq_dice(nlev_dice,nt_dice) real hu_dice(nlev_dice,nt_dice),hv_dice(nlev_dice,nt_dice) real w_dice(nlev_dice,nt_dice),omega_dice(nlev_dice,nt_dice) ! outputs: real tg_prof,shf_prof,lhf_prof,lwup_prof,swup_prof real ustar_prof,psurf_prof,ug_prof,vg_prof real ht_prof(nlev_dice),hq_prof(nlev_dice) real hu_prof(nlev_dice),hv_prof(nlev_dice) real w_prof(nlev_dice),omega_prof(nlev_dice) ! local: integer it_dice1, it_dice2,k real timeit,time_dice1,time_dice2,frac if (forcing_type.eq.7) then ! Check that initial day of the simulation consistent with Dice period: print *,'annee_ref=',annee_ref print *,'day1=',day1 print *,'day_ini_dice=',day_ini_dice if (annee_ref.ne.1999) then print*,'Pour Dice, annee_ref doit etre 1999' print*,'Changer annee_ref dans run.def' stop endif if (annee_ref.eq.1999 .and. day1.gt.day_ini_dice) then print*,'Dice a debute le 23 Oct 1999 (jour julien=296)' print*,'Changer dayref dans run.def',day1,day_ini_dice stop endif if (annee_ref.eq.1999 .and. day1.gt.day_ini_dice+2) then print*,'Dice a fini le 25 Oct 1999 (jour julien=298)' print*,'Changer dayref ou nday dans run.def',day1,day_ini_dice stop endif endif ! Determine timestep relative to the 1st day of TOGA-COARE: ! timeit=(day-day1)*86400. ! if (annee_ref.eq.1992) then ! timeit=(day-day_ini_dice)*86400. ! else ! timeit=(day+61.-1.)*86400. ! 61 days between Nov01 and Dec31 1992 ! endif timeit=(day-day_ini_dice)*86400 ! Determine the closest observation times: it_dice1=INT(timeit/dt_dice)+1 it_dice2=it_dice1 + 1 time_dice1=(it_dice1-1)*dt_dice time_dice2=(it_dice2-1)*dt_dice if (it_dice1 .ge. nt_dice) then write(*,*) 'PB-stop: day, it_dice1, it_dice2, timeit: ',day,it_dice1,it_dice2,timeit/86400. stop endif ! time interpolation: frac=(time_dice2-timeit)/(time_dice2-time_dice1) frac=max(frac,0.0) shf_prof = shf_dice(it_dice2)-frac*(shf_dice(it_dice2)-shf_dice(it_dice1)) lhf_prof = lhf_dice(it_dice2)-frac*(lhf_dice(it_dice2)-lhf_dice(it_dice1)) lwup_prof = lwup_dice(it_dice2)-frac*(lwup_dice(it_dice2)-lwup_dice(it_dice1)) swup_prof = swup_dice(it_dice2)-frac*(swup_dice(it_dice2)-swup_dice(it_dice1)) tg_prof = tg_dice(it_dice2)-frac*(tg_dice(it_dice2)-tg_dice(it_dice1)) ustar_prof = ustar_dice(it_dice2)-frac*(ustar_dice(it_dice2)-ustar_dice(it_dice1)) psurf_prof = psurf_dice(it_dice2)-frac*(psurf_dice(it_dice2)-psurf_dice(it_dice1)) ug_prof = ug_dice(it_dice2)-frac*(ug_dice(it_dice2)-ug_dice(it_dice1)) vg_prof = vg_dice(it_dice2)-frac*(vg_dice(it_dice2)-vg_dice(it_dice1)) ! print*, ! :'day,annee_ref,day_ini_dice,timeit,it_dice1,it_dice2,SST:', ! :day,annee_ref,day_ini_dice,timeit/86400.,it_dice1,it_dice2,ts_prof do k=1,nlev_dice ht_prof(k) = ht_dice(k,it_dice2)-frac*(ht_dice(k,it_dice2)-ht_dice(k,it_dice1)) hq_prof(k) = hq_dice(k,it_dice2)-frac*(hq_dice(k,it_dice2)-hq_dice(k,it_dice1)) hu_prof(k) = hu_dice(k,it_dice2)-frac*(hu_dice(k,it_dice2)-hu_dice(k,it_dice1)) hv_prof(k) = hv_dice(k,it_dice2)-frac*(hv_dice(k,it_dice2)-hv_dice(k,it_dice1)) w_prof(k) = w_dice(k,it_dice2)-frac*(w_dice(k,it_dice2)-w_dice(k,it_dice1)) omega_prof(k) = omega_dice(k,it_dice2)-frac*(omega_dice(k,it_dice2)-omega_dice(k,it_dice1)) enddo return END !====================================================================== SUBROUTINE interp_armcu_time(day,day1,annee_ref & & ,year_ini_armcu,day_ini_armcu,nt_armcu,dt_armcu & & ,nlev_armcu,fs_armcu,fl_armcu,at_armcu,rt_armcu & & ,aqt_armcu,fs_prof,fl_prof,at_prof,rt_prof,aqt_prof) implicit none !--------------------------------------------------------------------------------------- ! Time interpolation of a 2D field to the timestep corresponding to day ! ! day: current julian day (e.g. 717538.2) ! day1: first day of the simulation ! nt_armcu: total nb of data in the forcing (e.g. 31 for armcu) ! dt_armcu: total time interval (in sec) between 2 forcing data (e.g. 1/2h for armcu) ! fs= sensible flux ! fl= latent flux ! at,rt,aqt= advective and radiative tendencies !--------------------------------------------------------------------------------------- ! inputs: integer annee_ref integer nt_armcu,nlev_armcu integer year_ini_armcu real day, day1,day_ini_armcu,dt_armcu real fs_armcu(nt_armcu),fl_armcu(nt_armcu),at_armcu(nt_armcu) real rt_armcu(nt_armcu),aqt_armcu(nt_armcu) ! outputs: real fs_prof,fl_prof,at_prof,rt_prof,aqt_prof ! local: integer it_armcu1, it_armcu2,k real timeit,time_armcu1,time_armcu2,frac ! Check that initial day of the simulation consistent with ARMCU period: if (annee_ref.ne.1997 ) then print*,'Pour ARMCU, annee_ref doit etre 1997 ' print*,'Changer annee_ref dans run.def' stop endif timeit=(day-day_ini_armcu)*86400 ! Determine the closest observation times: it_armcu1=INT(timeit/dt_armcu)+1 it_armcu2=it_armcu1 + 1 time_armcu1=(it_armcu1-1)*dt_armcu time_armcu2=(it_armcu2-1)*dt_armcu print *,'timeit day day_ini_armcu',timeit,day,day_ini_armcu print *,'it_armcu1,it_armcu2,time_armcu1,time_armcu2', & & it_armcu1,it_armcu2,time_armcu1,time_armcu2 if (it_armcu1 .ge. nt_armcu) then write(*,*) 'PB-stop: day, it_armcu1, it_armcu2, timeit: ' & & ,day,it_armcu1,it_armcu2,timeit/86400. stop endif ! time interpolation: frac=(time_armcu2-timeit)/(time_armcu2-time_armcu1) frac=max(frac,0.0) fs_prof = fs_armcu(it_armcu2) & & -frac*(fs_armcu(it_armcu2)-fs_armcu(it_armcu1)) fl_prof = fl_armcu(it_armcu2) & & -frac*(fl_armcu(it_armcu2)-fl_armcu(it_armcu1)) at_prof = at_armcu(it_armcu2) & & -frac*(at_armcu(it_armcu2)-at_armcu(it_armcu1)) rt_prof = rt_armcu(it_armcu2) & & -frac*(rt_armcu(it_armcu2)-rt_armcu(it_armcu1)) aqt_prof = aqt_armcu(it_armcu2) & & -frac*(aqt_armcu(it_armcu2)-aqt_armcu(it_armcu1)) print*, & &'day,annee_ref,day_ini_armcu,timeit,it_armcu1,it_armcu2,SST:', & &day,annee_ref,day_ini_armcu,timeit/86400.,it_armcu1, & &it_armcu2,fs_prof,fl_prof,at_prof,rt_prof,aqt_prof return END !===================================================================== subroutine readprofiles(nlev_max,kmax,ntrac,height, & & thlprof,qtprof,uprof, & & vprof,e12prof,ugprof,vgprof, & & wfls,dqtdxls,dqtdyls,dqtdtls, & & thlpcar,tracer,nt1,nt2) implicit none integer nlev_max,kmax,kmax2,ntrac logical :: llesread = .true. real height(nlev_max),thlprof(nlev_max),qtprof(nlev_max), & & uprof(nlev_max),vprof(nlev_max),e12prof(nlev_max), & & ugprof(nlev_max),vgprof(nlev_max),wfls(nlev_max), & & dqtdxls(nlev_max),dqtdyls(nlev_max),dqtdtls(nlev_max), & & thlpcar(nlev_max),tracer(nlev_max,ntrac) integer, parameter :: ilesfile=1 integer :: ierr,k,itrac,nt1,nt2 if(.not.(llesread)) return open (ilesfile,file='prof.inp.001',status='old',iostat=ierr) if (ierr /= 0) stop 'ERROR:Prof.inp does not exist' read (ilesfile,*) kmax do k=1,kmax read (ilesfile,*) height(k),thlprof(k),qtprof (k), & & uprof (k),vprof (k),e12prof(k) enddo close(ilesfile) open(ilesfile,file='lscale.inp.001',status='old',iostat=ierr) if (ierr /= 0) stop 'ERROR:Lscale.inp does not exist' read (ilesfile,*) kmax2 if (kmax .ne. kmax2) then print *, 'fichiers prof.inp et lscale.inp incompatibles :' print *, 'nbre de niveaux : ',kmax,' et ',kmax2 stop 'lecture profiles' endif do k=1,kmax read (ilesfile,*) height(k),ugprof(k),vgprof(k),wfls(k), & & dqtdxls(k),dqtdyls(k),dqtdtls(k),thlpcar(k) end do close(ilesfile) open(ilesfile,file='trac.inp.001',status='old',iostat=ierr) if (ierr /= 0) then print*,'WARNING : trac.inp does not exist' else read (ilesfile,*) kmax2,nt1,nt2 if (nt2>ntrac) then stop'Augmenter le nombre de traceurs dans traceur.def' endif if (kmax .ne. kmax2) then print *, 'fichiers prof.inp et lscale.inp incompatibles :' print *, 'nbre de niveaux : ',kmax,' et ',kmax2 stop 'lecture profiles' endif do k=1,kmax read (ilesfile,*) height(k),(tracer(k,itrac),itrac=nt1,nt2) end do close(ilesfile) endif return end !====================================================================== subroutine readprofile_sandu(nlev_max,kmax,height,pprof,tprof, & & thlprof,qprof,uprof,vprof,wprof,omega,o3mmr) !====================================================================== implicit none integer nlev_max,kmax logical :: llesread = .true. real height(nlev_max),pprof(nlev_max),tprof(nlev_max) real thlprof(nlev_max) real qprof(nlev_max),uprof(nlev_max),vprof(nlev_max) real wprof(nlev_max),omega(nlev_max),o3mmr(nlev_max) integer, parameter :: ilesfile=1 integer :: k,ierr if(.not.(llesread)) return open (ilesfile,file='prof.inp.001',status='old',iostat=ierr) if (ierr /= 0) stop 'ERROR:Prof.inp does not exist' read (ilesfile,*) kmax do k=1,kmax read (ilesfile,*) height(k),pprof(k), tprof(k),thlprof(k), & & qprof (k),uprof(k), vprof(k), wprof(k), & & omega (k),o3mmr(k) enddo close(ilesfile) return end !====================================================================== subroutine readprofile_astex(nlev_max,kmax,height,pprof,tprof, & & thlprof,qvprof,qlprof,qtprof,uprof,vprof,wprof,tkeprof,o3mmr) !====================================================================== implicit none integer nlev_max,kmax logical :: llesread = .true. real height(nlev_max),pprof(nlev_max),tprof(nlev_max), & & thlprof(nlev_max),qlprof(nlev_max),qtprof(nlev_max), & & qvprof(nlev_max),uprof(nlev_max),vprof(nlev_max), & & wprof(nlev_max),tkeprof(nlev_max),o3mmr(nlev_max) integer, parameter :: ilesfile=1 integer :: ierr,k if(.not.(llesread)) return open (ilesfile,file='prof.inp.001',status='old',iostat=ierr) if (ierr /= 0) stop 'ERROR:Prof.inp does not exist' read (ilesfile,*) kmax do k=1,kmax read (ilesfile,*) height(k),pprof(k), tprof(k),thlprof(k), & & qvprof (k),qlprof (k),qtprof (k), & & uprof(k), vprof(k), wprof(k),tkeprof(k),o3mmr(k) enddo close(ilesfile) return end !====================================================================== subroutine readprofile_armcu(nlev_max,kmax,height,pprof,uprof, & & vprof,thetaprof,tprof,qvprof,rvprof,aprof,bprof) !====================================================================== implicit none integer nlev_max,kmax logical :: llesread = .true. real height(nlev_max),pprof(nlev_max),tprof(nlev_max) real thetaprof(nlev_max),rvprof(nlev_max) real qvprof(nlev_max),uprof(nlev_max),vprof(nlev_max) real aprof(nlev_max+1),bprof(nlev_max+1) integer, parameter :: ilesfile=1 integer, parameter :: ifile=2 integer :: ierr,jtot,k if(.not.(llesread)) return ! Read profiles at full levels IF(nlev_max.EQ.19) THEN open (ilesfile,file='prof.inp.19',status='old',iostat=ierr) print *,'On ouvre prof.inp.19' ELSE open (ilesfile,file='prof.inp.40',status='old',iostat=ierr) print *,'On ouvre prof.inp.40' ENDIF if (ierr /= 0) stop 'ERROR:Prof.inp does not exist' read (ilesfile,*) kmax do k=1,kmax read (ilesfile,*) height(k) ,pprof(k), uprof(k), vprof(k), & & thetaprof(k) ,tprof(k), qvprof(k),rvprof(k) enddo close(ilesfile) ! Vertical coordinates half levels for eta-coordinates (plev = alpha + beta * psurf) IF(nlev_max.EQ.19) THEN open (ifile,file='proh.inp.19',status='old',iostat=ierr) print *,'On ouvre proh.inp.19' if (ierr /= 0) stop 'ERROR:Proh.inp.19 does not exist' ELSE open (ifile,file='proh.inp.40',status='old',iostat=ierr) print *,'On ouvre proh.inp.40' if (ierr /= 0) stop 'ERROR:Proh.inp.40 does not exist' ENDIF read (ifile,*) kmax do k=1,kmax read (ifile,*) jtot,aprof(k),bprof(k) enddo close(ifile) return end !===================================================================== subroutine read_fire(fich_fire,nlevel,ntime & & ,zz,thl,qt,u,v,tke & & ,ug,vg,wls,dqtdx,dqtdy,dqtdt,thl_rad) !program reading forcings of the FIRE case study implicit none #include "netcdf.inc" integer ntime,nlevel character*80 :: fich_fire real*8 zz(nlevel) real*8 thl(nlevel) real*8 qt(nlevel),u(nlevel) real*8 v(nlevel),tke(nlevel) real*8 ug(nlevel,ntime),vg(nlevel,ntime),wls(nlevel,ntime) real*8 dqtdx(nlevel,ntime),dqtdy(nlevel,ntime) real*8 dqtdt(nlevel,ntime),thl_rad(nlevel,ntime) integer nid, ierr integer nbvar3d parameter(nbvar3d=30) integer var3didin(nbvar3d) ierr = NF_OPEN(fich_fire,NF_NOWRITE,nid) if (ierr.NE.NF_NOERR) then write(*,*) 'ERROR: Pb opening forcings nc file ' write(*,*) NF_STRERROR(ierr) stop "" endif ierr=NF_INQ_VARID(nid,"zz",var3didin(1)) if(ierr/=NF_NOERR) then write(*,*) NF_STRERROR(ierr) stop 'lev' endif ierr=NF_INQ_VARID(nid,"thetal",var3didin(2)) if(ierr/=NF_NOERR) then write(*,*) NF_STRERROR(ierr) stop 'temp' endif ierr=NF_INQ_VARID(nid,"qt",var3didin(3)) if(ierr/=NF_NOERR) then write(*,*) NF_STRERROR(ierr) stop 'qv' endif ierr=NF_INQ_VARID(nid,"u",var3didin(4)) if(ierr/=NF_NOERR) then write(*,*) NF_STRERROR(ierr) stop 'u' endif ierr=NF_INQ_VARID(nid,"v",var3didin(5)) if(ierr/=NF_NOERR) then write(*,*) NF_STRERROR(ierr) stop 'v' endif ierr=NF_INQ_VARID(nid,"tke",var3didin(6)) if(ierr/=NF_NOERR) then write(*,*) NF_STRERROR(ierr) stop 'tke' endif ierr=NF_INQ_VARID(nid,"ugeo",var3didin(7)) if(ierr/=NF_NOERR) then write(*,*) NF_STRERROR(ierr) stop 'ug' endif ierr=NF_INQ_VARID(nid,"vgeo",var3didin(8)) if(ierr/=NF_NOERR) then write(*,*) NF_STRERROR(ierr) stop 'vg' endif ierr=NF_INQ_VARID(nid,"wls",var3didin(9)) if(ierr/=NF_NOERR) then write(*,*) NF_STRERROR(ierr) stop 'wls' endif ierr=NF_INQ_VARID(nid,"dqtdx",var3didin(10)) if(ierr/=NF_NOERR) then write(*,*) NF_STRERROR(ierr) stop 'dqtdx' endif ierr=NF_INQ_VARID(nid,"dqtdy",var3didin(11)) if(ierr/=NF_NOERR) then write(*,*) NF_STRERROR(ierr) stop 'dqtdy' endif ierr=NF_INQ_VARID(nid,"dqtdt",var3didin(12)) if(ierr/=NF_NOERR) then write(*,*) NF_STRERROR(ierr) stop 'dqtdt' endif ierr=NF_INQ_VARID(nid,"thl_rad",var3didin(13)) if(ierr/=NF_NOERR) then write(*,*) NF_STRERROR(ierr) stop 'thl_rad' endif !dimensions lecture ! call catchaxis(nid,ntime,nlevel,time,z,ierr) #ifdef NC_DOUBLE ierr = NF_GET_VAR_DOUBLE(nid,var3didin(1),zz) #else ierr = NF_GET_VAR_REAL(nid,var3didin(1),zz) #endif if(ierr/=NF_NOERR) then write(*,*) NF_STRERROR(ierr) stop "getvarup" endif ! write(*,*)'lecture z ok',zz #ifdef NC_DOUBLE ierr = NF_GET_VAR_DOUBLE(nid,var3didin(2),thl) #else ierr = NF_GET_VAR_REAL(nid,var3didin(2),thl) #endif if(ierr/=NF_NOERR) then write(*,*) NF_STRERROR(ierr) stop "getvarup" endif ! write(*,*)'lecture thl ok',thl #ifdef NC_DOUBLE ierr = NF_GET_VAR_DOUBLE(nid,var3didin(3),qt) #else ierr = NF_GET_VAR_REAL(nid,var3didin(3),qt) #endif if(ierr/=NF_NOERR) then write(*,*) NF_STRERROR(ierr) stop "getvarup" endif ! write(*,*)'lecture qt ok',qt #ifdef NC_DOUBLE ierr = NF_GET_VAR_DOUBLE(nid,var3didin(4),u) #else ierr = NF_GET_VAR_REAL(nid,var3didin(4),u) #endif if(ierr/=NF_NOERR) then write(*,*) NF_STRERROR(ierr) stop "getvarup" endif ! write(*,*)'lecture u ok',u #ifdef NC_DOUBLE ierr = NF_GET_VAR_DOUBLE(nid,var3didin(5),v) #else ierr = NF_GET_VAR_REAL(nid,var3didin(5),v) #endif if(ierr/=NF_NOERR) then write(*,*) NF_STRERROR(ierr) stop "getvarup" endif ! write(*,*)'lecture v ok',v #ifdef NC_DOUBLE ierr = NF_GET_VAR_DOUBLE(nid,var3didin(6),tke) #else ierr = NF_GET_VAR_REAL(nid,var3didin(6),tke) #endif if(ierr/=NF_NOERR) then write(*,*) NF_STRERROR(ierr) stop "getvarup" endif ! write(*,*)'lecture tke ok',tke #ifdef NC_DOUBLE ierr = NF_GET_VAR_DOUBLE(nid,var3didin(7),ug) #else ierr = NF_GET_VAR_REAL(nid,var3didin(7),ug) #endif if(ierr/=NF_NOERR) then write(*,*) NF_STRERROR(ierr) stop "getvarup" endif ! write(*,*)'lecture ug ok',ug #ifdef NC_DOUBLE ierr = NF_GET_VAR_DOUBLE(nid,var3didin(8),vg) #else ierr = NF_GET_VAR_REAL(nid,var3didin(8),vg) #endif if(ierr/=NF_NOERR) then write(*,*) NF_STRERROR(ierr) stop "getvarup" endif ! write(*,*)'lecture vg ok',vg #ifdef NC_DOUBLE ierr = NF_GET_VAR_DOUBLE(nid,var3didin(9),wls) #else ierr = NF_GET_VAR_REAL(nid,var3didin(9),wls) #endif if(ierr/=NF_NOERR) then write(*,*) NF_STRERROR(ierr) stop "getvarup" endif ! write(*,*)'lecture wls ok',wls #ifdef NC_DOUBLE ierr = NF_GET_VAR_DOUBLE(nid,var3didin(10),dqtdx) #else ierr = NF_GET_VAR_REAL(nid,var3didin(10),dqtdx) #endif if(ierr/=NF_NOERR) then write(*,*) NF_STRERROR(ierr) stop "getvarup" endif ! write(*,*)'lecture dqtdx ok',dqtdx #ifdef NC_DOUBLE ierr = NF_GET_VAR_DOUBLE(nid,var3didin(11),dqtdy) #else ierr = NF_GET_VAR_REAL(nid,var3didin(11),dqtdy) #endif if(ierr/=NF_NOERR) then write(*,*) NF_STRERROR(ierr) stop "getvarup" endif ! write(*,*)'lecture dqtdy ok',dqtdy #ifdef NC_DOUBLE ierr = NF_GET_VAR_DOUBLE(nid,var3didin(12),dqtdt) #else ierr = NF_GET_VAR_REAL(nid,var3didin(12),dqtdt) #endif if(ierr/=NF_NOERR) then write(*,*) NF_STRERROR(ierr) stop "getvarup" endif ! write(*,*)'lecture dqtdt ok',dqtdt #ifdef NC_DOUBLE ierr = NF_GET_VAR_DOUBLE(nid,var3didin(13),thl_rad) #else ierr = NF_GET_VAR_REAL(nid,var3didin(13),thl_rad) #endif if(ierr/=NF_NOERR) then write(*,*) NF_STRERROR(ierr) stop "getvarup" endif ! write(*,*)'lecture thl_rad ok',thl_rad return end subroutine read_fire !===================================================================== subroutine read_dice(fich_dice,nlevel,ntime & & ,zz,pres,th,qv,u,v,o3 & & ,shf,lhf,lwup,swup,tg,ustar,psurf,ug,vg & & ,hadvt,hadvq,hadvu,hadvv,w,omega) !program reading initial profils and forcings of the Dice case study implicit none #include "netcdf.inc" integer ntime,nlevel integer l,k character*80 :: fich_dice real*8 time(ntime) real*8 zz(nlevel) real*8 th(nlevel),pres(nlevel) real*8 qv(nlevel),u(nlevel),v(nlevel),o3(nlevel) real*8 shf(ntime),lhf(ntime),lwup(ntime),swup(ntime),tg(ntime) real*8 ustar(ntime),psurf(ntime),ug(ntime),vg(ntime) real*8 hadvt(nlevel,ntime),hadvq(nlevel,ntime),hadvu(nlevel,ntime) real*8 hadvv(nlevel,ntime),w(nlevel,ntime),omega(nlevel,ntime) integer nid, ierr integer nbvar3d parameter(nbvar3d=30) integer var3didin(nbvar3d) ierr = NF_OPEN(fich_dice,NF_NOWRITE,nid) if (ierr.NE.NF_NOERR) then write(*,*) 'ERROR: Pb opening forcings nc file ' write(*,*) NF_STRERROR(ierr) stop "" endif ierr=NF_INQ_VARID(nid,"height",var3didin(1)) if(ierr/=NF_NOERR) then write(*,*) NF_STRERROR(ierr) stop 'height' endif ierr=NF_INQ_VARID(nid,"pf",var3didin(11)) if(ierr/=NF_NOERR) then write(*,*) NF_STRERROR(ierr) stop 'pf' endif ierr=NF_INQ_VARID(nid,"theta",var3didin(12)) if(ierr/=NF_NOERR) then write(*,*) NF_STRERROR(ierr) stop 'theta' endif ierr=NF_INQ_VARID(nid,"qv",var3didin(13)) if(ierr/=NF_NOERR) then write(*,*) NF_STRERROR(ierr) stop 'qv' endif ierr=NF_INQ_VARID(nid,"u",var3didin(14)) if(ierr/=NF_NOERR) then write(*,*) NF_STRERROR(ierr) stop 'u' endif ierr=NF_INQ_VARID(nid,"v",var3didin(15)) if(ierr/=NF_NOERR) then write(*,*) NF_STRERROR(ierr) stop 'v' endif ierr=NF_INQ_VARID(nid,"o3mmr",var3didin(16)) if(ierr/=NF_NOERR) then write(*,*) NF_STRERROR(ierr) stop 'o3' endif ierr=NF_INQ_VARID(nid,"shf",var3didin(2)) if(ierr/=NF_NOERR) then write(*,*) NF_STRERROR(ierr) stop 'shf' endif ierr=NF_INQ_VARID(nid,"lhf",var3didin(3)) if(ierr/=NF_NOERR) then write(*,*) NF_STRERROR(ierr) stop 'lhf' endif ierr=NF_INQ_VARID(nid,"lwup",var3didin(4)) if(ierr/=NF_NOERR) then write(*,*) NF_STRERROR(ierr) stop 'lwup' endif ierr=NF_INQ_VARID(nid,"swup",var3didin(5)) if(ierr/=NF_NOERR) then write(*,*) NF_STRERROR(ierr) stop 'dqtdx' endif ierr=NF_INQ_VARID(nid,"Tg",var3didin(6)) if(ierr/=NF_NOERR) then write(*,*) NF_STRERROR(ierr) stop 'Tg' endif ierr=NF_INQ_VARID(nid,"ustar",var3didin(7)) if(ierr/=NF_NOERR) then write(*,*) NF_STRERROR(ierr) stop 'ustar' endif ierr=NF_INQ_VARID(nid,"psurf",var3didin(8)) if(ierr/=NF_NOERR) then write(*,*) NF_STRERROR(ierr) stop 'psurf' endif ierr=NF_INQ_VARID(nid,"Ug",var3didin(9)) if(ierr/=NF_NOERR) then write(*,*) NF_STRERROR(ierr) stop 'Ug' endif ierr=NF_INQ_VARID(nid,"Vg",var3didin(10)) if(ierr/=NF_NOERR) then write(*,*) NF_STRERROR(ierr) stop 'Vg' endif ierr=NF_INQ_VARID(nid,"hadvT",var3didin(17)) if(ierr/=NF_NOERR) then write(*,*) NF_STRERROR(ierr) stop 'hadvT' endif ierr=NF_INQ_VARID(nid,"hadvq",var3didin(18)) if(ierr/=NF_NOERR) then write(*,*) NF_STRERROR(ierr) stop 'hadvq' endif ierr=NF_INQ_VARID(nid,"hadvu",var3didin(19)) if(ierr/=NF_NOERR) then write(*,*) NF_STRERROR(ierr) stop 'hadvu' endif ierr=NF_INQ_VARID(nid,"hadvv",var3didin(20)) if(ierr/=NF_NOERR) then write(*,*) NF_STRERROR(ierr) stop 'hadvv' endif ierr=NF_INQ_VARID(nid,"w",var3didin(21)) if(ierr/=NF_NOERR) then write(*,*) NF_STRERROR(ierr) stop 'w' endif ierr=NF_INQ_VARID(nid,"omega",var3didin(22)) if(ierr/=NF_NOERR) then write(*,*) NF_STRERROR(ierr) stop 'omega' endif !dimensions lecture ! call catchaxis(nid,ntime,nlevel,time,z,ierr) #ifdef NC_DOUBLE ierr = NF_GET_VAR_DOUBLE(nid,var3didin(1),zz) #else ierr = NF_GET_VAR_REAL(nid,var3didin(1),zz) #endif if(ierr/=NF_NOERR) then write(*,*) NF_STRERROR(ierr) stop "getvarup" endif ! write(*,*)'lecture zz ok',zz #ifdef NC_DOUBLE ierr = NF_GET_VAR_DOUBLE(nid,var3didin(11),pres) #else ierr = NF_GET_VAR_REAL(nid,var3didin(11),pres) #endif if(ierr/=NF_NOERR) then write(*,*) NF_STRERROR(ierr) stop "getvarup" endif ! write(*,*)'lecture pres ok',pres #ifdef NC_DOUBLE ierr = NF_GET_VAR_DOUBLE(nid,var3didin(12),th) #else ierr = NF_GET_VAR_REAL(nid,var3didin(12),th) #endif if(ierr/=NF_NOERR) then write(*,*) NF_STRERROR(ierr) stop "getvarup" endif ! write(*,*)'lecture th ok',th #ifdef NC_DOUBLE ierr = NF_GET_VAR_DOUBLE(nid,var3didin(13),qv) #else ierr = NF_GET_VAR_REAL(nid,var3didin(13),qv) #endif if(ierr/=NF_NOERR) then write(*,*) NF_STRERROR(ierr) stop "getvarup" endif ! write(*,*)'lecture qv ok',qv #ifdef NC_DOUBLE ierr = NF_GET_VAR_DOUBLE(nid,var3didin(14),u) #else ierr = NF_GET_VAR_REAL(nid,var3didin(14),u) #endif if(ierr/=NF_NOERR) then write(*,*) NF_STRERROR(ierr) stop "getvarup" endif ! write(*,*)'lecture u ok',u #ifdef NC_DOUBLE ierr = NF_GET_VAR_DOUBLE(nid,var3didin(15),v) #else ierr = NF_GET_VAR_REAL(nid,var3didin(15),v) #endif if(ierr/=NF_NOERR) then write(*,*) NF_STRERROR(ierr) stop "getvarup" endif ! write(*,*)'lecture v ok',v #ifdef NC_DOUBLE ierr = NF_GET_VAR_DOUBLE(nid,var3didin(16),o3) #else ierr = NF_GET_VAR_REAL(nid,var3didin(16),o3) #endif if(ierr/=NF_NOERR) then write(*,*) NF_STRERROR(ierr) stop "getvarup" endif ! write(*,*)'lecture o3 ok',o3 #ifdef NC_DOUBLE ierr = NF_GET_VAR_DOUBLE(nid,var3didin(2),shf) #else ierr = NF_GET_VAR_REAL(nid,var3didin(2),shf) #endif if(ierr/=NF_NOERR) then write(*,*) NF_STRERROR(ierr) stop "getvarup" endif ! write(*,*)'lecture shf ok',shf #ifdef NC_DOUBLE ierr = NF_GET_VAR_DOUBLE(nid,var3didin(3),lhf) #else ierr = NF_GET_VAR_REAL(nid,var3didin(3),lhf) #endif if(ierr/=NF_NOERR) then write(*,*) NF_STRERROR(ierr) stop "getvarup" endif ! write(*,*)'lecture lhf ok',lhf #ifdef NC_DOUBLE ierr = NF_GET_VAR_DOUBLE(nid,var3didin(4),lwup) #else ierr = NF_GET_VAR_REAL(nid,var3didin(4),lwup) #endif if(ierr/=NF_NOERR) then write(*,*) NF_STRERROR(ierr) stop "getvarup" endif ! write(*,*)'lecture lwup ok',lwup #ifdef NC_DOUBLE ierr = NF_GET_VAR_DOUBLE(nid,var3didin(5),swup) #else ierr = NF_GET_VAR_REAL(nid,var3didin(5),swup) #endif if(ierr/=NF_NOERR) then write(*,*) NF_STRERROR(ierr) stop "getvarup" endif ! write(*,*)'lecture swup ok',swup #ifdef NC_DOUBLE ierr = NF_GET_VAR_DOUBLE(nid,var3didin(6),tg) #else ierr = NF_GET_VAR_REAL(nid,var3didin(6),tg) #endif if(ierr/=NF_NOERR) then write(*,*) NF_STRERROR(ierr) stop "getvarup" endif ! write(*,*)'lecture tg ok',tg #ifdef NC_DOUBLE ierr = NF_GET_VAR_DOUBLE(nid,var3didin(7),ustar) #else ierr = NF_GET_VAR_REAL(nid,var3didin(7),ustar) #endif if(ierr/=NF_NOERR) then write(*,*) NF_STRERROR(ierr) stop "getvarup" endif ! write(*,*)'lecture ustar ok',ustar #ifdef NC_DOUBLE ierr = NF_GET_VAR_DOUBLE(nid,var3didin(8),psurf) #else ierr = NF_GET_VAR_REAL(nid,var3didin(8),psurf) #endif if(ierr/=NF_NOERR) then write(*,*) NF_STRERROR(ierr) stop "getvarup" endif ! write(*,*)'lecture psurf ok',psurf #ifdef NC_DOUBLE ierr = NF_GET_VAR_DOUBLE(nid,var3didin(9),ug) #else ierr = NF_GET_VAR_REAL(nid,var3didin(9),ug) #endif if(ierr/=NF_NOERR) then write(*,*) NF_STRERROR(ierr) stop "getvarup" endif ! write(*,*)'lecture ug ok',ug #ifdef NC_DOUBLE ierr = NF_GET_VAR_DOUBLE(nid,var3didin(10),vg) #else ierr = NF_GET_VAR_REAL(nid,var3didin(10),vg) #endif if(ierr/=NF_NOERR) then write(*,*) NF_STRERROR(ierr) stop "getvarup" endif ! write(*,*)'lecture vg ok',vg #ifdef NC_DOUBLE ierr = NF_GET_VAR_DOUBLE(nid,var3didin(17),hadvt) #else ierr = NF_GET_VAR_REAL(nid,var3didin(17),hadvt) #endif if(ierr/=NF_NOERR) then write(*,*) NF_STRERROR(ierr) stop "getvarup" endif ! write(*,*)'lecture hadvt ok',hadvt #ifdef NC_DOUBLE ierr = NF_GET_VAR_DOUBLE(nid,var3didin(18),hadvq) #else ierr = NF_GET_VAR_REAL(nid,var3didin(18),hadvq) #endif if(ierr/=NF_NOERR) then write(*,*) NF_STRERROR(ierr) stop "getvarup" endif ! write(*,*)'lecture hadvq ok',hadvq #ifdef NC_DOUBLE ierr = NF_GET_VAR_DOUBLE(nid,var3didin(19),hadvu) #else ierr = NF_GET_VAR_REAL(nid,var3didin(19),hadvu) #endif if(ierr/=NF_NOERR) then write(*,*) NF_STRERROR(ierr) stop "getvarup" endif ! write(*,*)'lecture hadvu ok',hadvu #ifdef NC_DOUBLE ierr = NF_GET_VAR_DOUBLE(nid,var3didin(20),hadvv) #else ierr = NF_GET_VAR_REAL(nid,var3didin(20),hadvv) #endif if(ierr/=NF_NOERR) then write(*,*) NF_STRERROR(ierr) stop "getvarup" endif ! write(*,*)'lecture hadvv ok',hadvv #ifdef NC_DOUBLE ierr = NF_GET_VAR_DOUBLE(nid,var3didin(21),w) #else ierr = NF_GET_VAR_REAL(nid,var3didin(21),w) #endif if(ierr/=NF_NOERR) then write(*,*) NF_STRERROR(ierr) stop "getvarup" endif ! write(*,*)'lecture w ok',w #ifdef NC_DOUBLE ierr = NF_GET_VAR_DOUBLE(nid,var3didin(22),omega) #else ierr = NF_GET_VAR_REAL(nid,var3didin(22),omega) #endif if(ierr/=NF_NOERR) then write(*,*) NF_STRERROR(ierr) stop "getvarup" endif ! write(*,*)'lecture omega ok',omega return end subroutine read_dice !===================================================================== ! Reads CIRC input files SUBROUTINE read_circ(nlev_circ,cf,lwp,iwp,reliq,reice,t,z,p,pm,h2o,o3,sza) parameter (ncm_1=49180) #include "YOMCST.h" real albsfc(ncm_1), albsfc_w(ncm_1) real cf(nlev_circ), icefra(nlev_circ), deice(nlev_circ), & reliq(nlev_circ), reice(nlev_circ), lwp(nlev_circ), iwp(nlev_circ) real t(nlev_circ+1), z(nlev_circ+1), dz(nlev_circ), p(nlev_circ+1) real aer_beta(nlev_circ), waer(nlev_circ), gaer(nlev_circ) real pm(nlev_circ), tm(nlev_circ), h2o(nlev_circ), o3(nlev_circ) real co2(nlev_circ), n2o(nlev_circ), co(nlev_circ), ch4(nlev_circ), & o2(nlev_circ), ccl4(nlev_circ), f11(nlev_circ), f12(nlev_circ) ! za= zenital angle ! sza= cosinus angle zenital real wavn(ncm_1), ssf(ncm_1),za,sza integer nlev ! Open the files open (11, file='Tsfc_sza_nlev_case.txt', status='old') open (12, file='level_input_case.txt', status='old') open (13, file='layer_input_case.txt', status='old') open (14, file='aerosol_input_case.txt', status='old') open (15, file='cloud_input_case.txt', status='old') open (16, file='sfcalbedo_input_case.txt', status='old') ! Read scalar information do iskip=1,5 read (11, *) enddo read (11, '(i8)') nlev read (11, '(f10.2)') tsfc read (11, '(f10.2)') za read (11, '(f10.4)') sw_dn_toa sza=cos(za/180.*RPI) print *,'nlev,tsfc,sza,sw_dn_toa,RPI',nlev,tsfc,sza,sw_dn_toa,RPI close(11) ! Read level information read (12, *) do il=1,nlev read (12, 302) ilev, z(il), p(il), t(il) z(il)=z(il)*1000. ! z donne en km p(il)=p(il)*100. ! p donne en mb enddo 302 format (i8, f8.3, 2f9.2) close(12) ! Read layer information (midpoint values) do iskip=1,3 read (13, *) enddo do il=1,nlev-1 read (13, 303) ilev,pm(il),tm(il),h2o(il),co2(il),o3(il), & n2o(il),co(il),ch4(il),o2(il),ccl4(il), & f11(il),f12(il) pm(il)=pm(il)*100. enddo 303 format (i8, 2f9.2, 10(2x,e13.7)) close(13) ! Read aerosol layer information do iskip=1,3 read (14, *) enddo read (14, '(f10.2)') aer_alpha read (14, *) read (14, *) do il=1,nlev-1 read (14, 304) ilev, aer_beta(il), waer(il), gaer(il) enddo 304 format (i8, f9.5, 2f8.3) close(14) ! Read cloud information do iskip=1,3 read (15, *) enddo do il=1,nlev-1 read (15, 305) ilev, cf(il), lwp(il), iwp(il), reliq(il), reice(il) lwp(il)=lwp(il)/1000. ! lwp donne en g/kg iwp(il)=iwp(il)/1000. ! iwp donne en g/kg reliq(il)=reliq(il)/1000000. ! reliq donne en microns reice(il)=reice(il)/1000000. ! reice donne en microns enddo 305 format (i8, f8.3, 4f9.2) close(15) ! Read surface albedo (weighted & unweighted) and spectral solar irradiance do iskip=1,6 read (16, *) enddo do icm_1=1,ncm_1 read (16, 306) wavn(icm_1), albsfc(icm_1), albsfc_w(icm_1), ssf(icm_1) enddo 306 format(f10.1, 2f12.5, f14.8) close(16) return end subroutine read_circ !===================================================================== ! Reads RTMIP input files SUBROUTINE read_rtmip(nlev_rtmip,play,plev,t,h2o,o3) #include "YOMCST.h" real t(nlev_rtmip), pt(nlev_rtmip),pb(nlev_rtmip),h2o(nlev_rtmip), o3(nlev_rtmip) real temp(nlev_rtmip), play(nlev_rtmip),ovap(nlev_rtmip), oz(nlev_rtmip),plev(nlev_rtmip+1) integer nlev ! Open the files open (11, file='low_resolution_profile.txt', status='old') ! Read level information read (11, *) do il=1,nlev_rtmip read (11, 302) pt(il), pb(il), t(il),h2o(il),o3(il) enddo do il=1,nlev_rtmip play(il)=pt(nlev_rtmip-il+1)*100. ! p donne en mb temp(il)=t(nlev_rtmip-il+1) ovap(il)=h2o(nlev_rtmip-il+1) oz(il)=o3(nlev_rtmip-il+1) enddo do il=1,39 plev(il)=play(il)+(play(il+1)-play(il))/2. print *,'il p t ovap oz=',il,plev(il),temp(il),ovap(il),oz(il) enddo plev(41)=101300. 302 format (e16.10,3x,e16.10,3x,e16.10,3x,e12.6,3x,e12.6) close(12) return end subroutine read_rtmip !===================================================================== ! Subroutines for nudging Subroutine Nudge_RHT_init (paprs,pplay,t,q,t_targ,rh_targ) ! ======================================================== USE dimphy implicit none ! ======================================================== REAL paprs(klon,klevp1) REAL pplay(klon,klev) ! ! Variables d'etat REAL t(klon,klev) REAL q(klon,klev) ! ! Profiles cible REAL t_targ(klon,klev) REAL rh_targ(klon,klev) ! INTEGER k,i REAL zx_qs ! Declaration des constantes et des fonctions thermodynamiques ! include "YOMCST.h" include "YOETHF.h" ! ! ---------------------------------------- ! Statement functions include "FCTTRE.h" ! ---------------------------------------- ! DO k = 1,klev DO i = 1,klon t_targ(i,k) = t(i,k) IF (t(i,k).LT.RTT) THEN zx_qs = qsats(t(i,k))/(pplay(i,k)) ELSE zx_qs = qsatl(t(i,k))/(pplay(i,k)) ENDIF rh_targ(i,k) = q(i,k)/zx_qs ENDDO ENDDO print *, 't_targ',t_targ print *, 'rh_targ',rh_targ ! ! RETURN END Subroutine Nudge_UV_init (paprs,pplay,u,v,u_targ,v_targ) ! ======================================================== USE dimphy implicit none ! ======================================================== REAL paprs(klon,klevp1) REAL pplay(klon,klev) ! ! Variables d'etat REAL u(klon,klev) REAL v(klon,klev) ! ! Profiles cible REAL u_targ(klon,klev) REAL v_targ(klon,klev) ! INTEGER k,i ! DO k = 1,klev DO i = 1,klon u_targ(i,k) = u(i,k) v_targ(i,k) = v(i,k) ENDDO ENDDO print *, 'u_targ',u_targ print *, 'v_targ',v_targ ! ! RETURN END Subroutine Nudge_RHT (dtime,paprs,pplay,t_targ,rh_targ,t,q, & & d_t,d_q) ! ======================================================== USE dimphy implicit none ! ======================================================== REAL dtime REAL paprs(klon,klevp1) REAL pplay(klon,klev) ! ! Variables d'etat REAL t(klon,klev) REAL q(klon,klev) ! ! Tendances REAL d_t(klon,klev) REAL d_q(klon,klev) ! ! Profiles cible REAL t_targ(klon,klev) REAL rh_targ(klon,klev) ! ! Temps de relaxation REAL tau !c DATA tau /3600./ !! DATA tau /5400./ DATA tau /1800./ ! INTEGER k,i REAL zx_qs, rh, tnew, d_rh ! Declaration des constantes et des fonctions thermodynamiques ! include "YOMCST.h" include "YOETHF.h" ! ! ---------------------------------------- ! Statement functions include "FCTTRE.h" ! ---------------------------------------- ! print *,'dtime, tau ',dtime,tau print *, 't_targ',t_targ print *, 'rh_targ',rh_targ print *,'temp ',t print *,'hum ',q DO k = 1,klev DO i = 1,klon !! IF (paprs(i,1)-pplay(i,k) .GT. 10000.) THEN IF (t(i,k).LT.RTT) THEN zx_qs = qsats(t(i,k))/(pplay(i,k)) ELSE zx_qs = qsatl(t(i,k))/(pplay(i,k)) ENDIF rh = q(i,k)/zx_qs ! d_t(i,k) = d_t(i,k) + 1./tau*(t_targ(i,k)-t(i,k)) d_rh = 1./tau*(rh_targ(i,k)-rh) ! tnew = t(i,k)+d_t(i,k) IF (tnew.LT.RTT) THEN zx_qs = qsats(tnew)/(pplay(i,k)) ELSE zx_qs = qsatl(tnew)/(pplay(i,k)) ENDIF d_q(i,k) = d_q(i,k) + d_rh*zx_qs ! print *,' k,d_t,rh,d_rh,d_q ', & k,d_t(i,k),rh,d_rh,d_q(i,k) !! ENDIF ! ENDDO ENDDO ! RETURN END Subroutine Nudge_UV (dtime,paprs,pplay,u_targ,v_targ,u,v, & & d_u,d_v) ! ======================================================== USE dimphy implicit none ! ======================================================== REAL dtime REAL paprs(klon,klevp1) REAL pplay(klon,klev) ! ! Variables d'etat REAL u(klon,klev) REAL v(klon,klev) ! ! Tendances REAL d_u(klon,klev) REAL d_v(klon,klev) ! ! Profiles cible REAL u_targ(klon,klev) REAL v_targ(klon,klev) ! ! Temps de relaxation REAL tau !c DATA tau /3600./ DATA tau /5400./ ! INTEGER k,i ! print *,'dtime, tau ',dtime,tau print *, 'u_targ',u_targ print *, 'v_targ',v_targ print *,'zonal velocity ',u print *,'meridional velocity ',v DO k = 1,klev DO i = 1,klon IF (paprs(i,1)-pplay(i,k) .GT. 10000.) THEN ! d_u(i,k) = d_u(i,k) + 1./tau*(u_targ(i,k)-u(i,k)) d_v(i,k) = d_v(i,k) + 1./tau*(v_targ(i,k)-v(i,k)) ! print *,' k,u,d_u,v,d_v ', & k,u(i,k),d_u(i,k),v(i,k),d_v(i,k) ENDIF ! ENDDO ENDDO ! RETURN END 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  integer, parameter :: ilesfile=1
3239  integer :: ierr,k,itrac,nt1,nt2
3240 
3241  if(.not.(llesread)) return
3242 
3243  open (ilesfile,file='prof.inp.001',status='old',iostat=ierr)
3244  if (ierr /= 0) stop 'ERROR:Prof.inp does not exist'
3245  read (ilesfile,*) kmax
3246  do k=1,kmax
3247  read (ilesfile,*) height(k),thlprof(k),qtprof (k), &
3248  & uprof (k),vprof (k),e12prof(k)
3249  enddo
3250  close(ilesfile)
3251 
3252  open(ilesfile,file='lscale.inp.001',status='old',iostat=ierr)
3253  if (ierr /= 0) stop 'ERROR:Lscale.inp does not exist'
3254  read (ilesfile,*) kmax2
3255  if (kmax .ne. kmax2) then
3256  print *, 'fichiers prof.inp et lscale.inp incompatibles :'
3257  print *, 'nbre de niveaux : ',kmax,' et ',kmax2
3258  stop 'lecture profiles'
3259  endif
3260  do k=1,kmax
3261  read (ilesfile,*) height(k),ugprof(k),vgprof(k),wfls(k), &
3262  & dqtdxls(k),dqtdyls(k),dqtdtls(k),thlpcar(k)
3263  end do
3264  close(ilesfile)
3265 
3266  open(ilesfile,file='trac.inp.001',status='old',iostat=ierr)
3267  if (ierr /= 0) then
3268  print*,'WARNING : trac.inp does not exist'
3269  else
3270  read (ilesfile,*) kmax2,nt1,nt2
3271  if (nt2>ntrac) then
3272  stop'Augmenter le nombre de traceurs dans traceur.def'
3273  endif
3274  if (kmax .ne. kmax2) then
3275  print *, 'fichiers prof.inp et lscale.inp incompatibles :'
3276  print *, 'nbre de niveaux : ',kmax,' et ',kmax2
3277  stop 'lecture profiles'
3278  endif
3279  do k=1,kmax
3280  read (ilesfile,*) height(k),(tracer(k,itrac),itrac=nt1,nt2)
3281  end do
3282  close(ilesfile)
3283  endif
3284 
3285  return
3286  end
3287 !======================================================================
3288  subroutine readprofile_sandu(nlev_max,kmax,height,pprof,tprof, &
3289  & thlprof,qprof,uprof,vprof,wprof,omega,o3mmr)
3290 !======================================================================
3291  implicit none
3292 
3293  integer nlev_max,kmax
3294  logical :: llesread = .true.
3295 
3296  real height(nlev_max),pprof(nlev_max),tprof(nlev_max)
3297  real thlprof(nlev_max)
3298  real qprof(nlev_max),uprof(nlev_max),vprof(nlev_max)
3299  real wprof(nlev_max),omega(nlev_max),o3mmr(nlev_max)
3300 
3301  integer, parameter :: ilesfile=1
3302  integer :: k,ierr
3303 
3304  if(.not.(llesread)) return
3305 
3306  open (ilesfile,file='prof.inp.001',status='old',iostat=ierr)
3307  if (ierr /= 0) stop 'ERROR:Prof.inp does not exist'
3308  read (ilesfile,*) kmax
3309  do k=1,kmax
3310  read (ilesfile,*) height(k),pprof(k), tprof(k),thlprof(k), &
3311  & qprof (k),uprof(k), vprof(k), wprof(k), &
3312  & omega (k),o3mmr(k)
3313  enddo
3314  close(ilesfile)
3315 
3316  return
3317  end
3318 
3319 !======================================================================
3320  subroutine readprofile_astex(nlev_max,kmax,height,pprof,tprof, &
3321  & thlprof,qvprof,qlprof,qtprof,uprof,vprof,wprof,tkeprof,o3mmr)
3322 !======================================================================
3323  implicit none
3324 
3325  integer nlev_max,kmax
3326  logical :: llesread = .true.
3327 
3328  real height(nlev_max),pprof(nlev_max),tprof(nlev_max), &
3329  & thlprof(nlev_max),qlprof(nlev_max),qtprof(nlev_max), &
3330  & qvprof(nlev_max),uprof(nlev_max),vprof(nlev_max), &
3331  & wprof(nlev_max),tkeprof(nlev_max),o3mmr(nlev_max)
3332 
3333  integer, parameter :: ilesfile=1
3334  integer :: ierr,k
3335 
3336  if(.not.(llesread)) return
3337 
3338  open (ilesfile,file='prof.inp.001',status='old',iostat=ierr)
3339  if (ierr /= 0) stop 'ERROR:Prof.inp does not exist'
3340  read (ilesfile,*) kmax
3341  do k=1,kmax
3342  read (ilesfile,*) height(k),pprof(k), tprof(k),thlprof(k), &
3343  & qvprof (k),qlprof (k),qtprof (k), &
3344  & uprof(k), vprof(k), wprof(k),tkeprof(k),o3mmr(k)
3345  enddo
3346  close(ilesfile)
3347 
3348  return
3349  end
3350 
3351 
3352 
3353 !======================================================================
3354  subroutine readprofile_armcu(nlev_max,kmax,height,pprof,uprof, &
3355  & vprof,thetaprof,tprof,qvprof,rvprof,aprof,bprof)
3356 !======================================================================
3357  implicit none
3358 
3359  integer nlev_max,kmax
3360  logical :: llesread = .true.
3361 
3362  real height(nlev_max),pprof(nlev_max),tprof(nlev_max)
3363  real thetaprof(nlev_max),rvprof(nlev_max)
3364  real qvprof(nlev_max),uprof(nlev_max),vprof(nlev_max)
3365  real aprof(nlev_max+1),bprof(nlev_max+1)
3366 
3367  integer, parameter :: ilesfile=1
3368  integer, parameter :: ifile=2
3369  integer :: ierr,jtot,k
3370 
3371  if(.not.(llesread)) return
3372 
3373 ! Read profiles at full levels
3374  IF(nlev_max.EQ.19) THEN
3375  open (ilesfile,file='prof.inp.19',status='old',iostat=ierr)
3376  print *,'On ouvre prof.inp.19'
3377  ELSE
3378  open (ilesfile,file='prof.inp.40',status='old',iostat=ierr)
3379  print *,'On ouvre prof.inp.40'
3380  ENDIF
3381  if (ierr /= 0) stop 'ERROR:Prof.inp does not exist'
3382  read (ilesfile,*) kmax
3383  do k=1,kmax
3384  read (ilesfile,*) height(k) ,pprof(k), uprof(k), vprof(k), &
3385  & thetaprof(k) ,tprof(k), qvprof(k),rvprof(k)
3386  enddo
3387  close(ilesfile)
3388 
3389 ! Vertical coordinates half levels for eta-coordinates (plev = alpha + beta * psurf)
3390  IF(nlev_max.EQ.19) THEN
3391  open (ifile,file='proh.inp.19',status='old',iostat=ierr)
3392  print *,'On ouvre proh.inp.19'
3393  if (ierr /= 0) stop 'ERROR:Proh.inp.19 does not exist'
3394  ELSE
3395  open (ifile,file='proh.inp.40',status='old',iostat=ierr)
3396  print *,'On ouvre proh.inp.40'
3397  if (ierr /= 0) stop 'ERROR:Proh.inp.40 does not exist'
3398  ENDIF
3399  read (ifile,*) kmax
3400  do k=1,kmax
3401  read (ifile,*) jtot,aprof(k),bprof(k)
3402  enddo
3403  close(ifile)
3404 
3405  return
3406  end
3407 
3408 !=====================================================================
3409  subroutine read_fire(fich_fire,nlevel,ntime &
3410  & ,zz,thl,qt,u,v,tke &
3411  & ,ug,vg,wls,dqtdx,dqtdy,dqtdt,thl_rad)
3412 
3413 !program reading forcings of the FIRE case study
3414 
3415 
3416  implicit none
3417 
3418 #include "netcdf.inc"
3419 
3420  integer ntime,nlevel
3421  character*80 :: fich_fire
3422  real*8 zz(nlevel)
3423 
3424  real*8 thl(nlevel)
3425  real*8 qt(nlevel),u(nlevel)
3426  real*8 v(nlevel),tke(nlevel)
3427  real*8 ug(nlevel,ntime),vg(nlevel,ntime),wls(nlevel,ntime)
3428  real*8 dqtdx(nlevel,ntime),dqtdy(nlevel,ntime)
3429  real*8 dqtdt(nlevel,ntime),thl_rad(nlevel,ntime)
3430 
3431  integer nid, ierr
3432  integer nbvar3d
3433  parameter(nbvar3d=30)
3434  integer var3didin(nbvar3d)
3435 
3436  ierr = NF_OPEN(fich_fire,NF_NOWRITE,nid)
3437  if (ierr.NE.NF_NOERR) then
3438  write(*,*) 'ERROR: Pb opening forcings nc file '
3439  write(*,*) NF_STRERROR(ierr)
3440  stop ""
3441  endif
3442 
3443 
3444  ierr=NF_INQ_VARID(nid,"zz",var3didin(1))
3445  if(ierr/=NF_NOERR) then
3446  write(*,*) NF_STRERROR(ierr)
3447  stop 'lev'
3448  endif
3449 
3450 
3451  ierr=NF_INQ_VARID(nid,"thetal",var3didin(2))
3452  if(ierr/=NF_NOERR) then
3453  write(*,*) NF_STRERROR(ierr)
3454  stop 'temp'
3455  endif
3456 
3457  ierr=NF_INQ_VARID(nid,"qt",var3didin(3))
3458  if(ierr/=NF_NOERR) then
3459  write(*,*) NF_STRERROR(ierr)
3460  stop 'qv'
3461  endif
3462 
3463  ierr=NF_INQ_VARID(nid,"u",var3didin(4))
3464  if(ierr/=NF_NOERR) then
3465  write(*,*) NF_STRERROR(ierr)
3466  stop 'u'
3467  endif
3468 
3469  ierr=NF_INQ_VARID(nid,"v",var3didin(5))
3470  if(ierr/=NF_NOERR) then
3471  write(*,*) NF_STRERROR(ierr)
3472  stop 'v'
3473  endif
3474 
3475  ierr=NF_INQ_VARID(nid,"tke",var3didin(6))
3476  if(ierr/=NF_NOERR) then
3477  write(*,*) NF_STRERROR(ierr)
3478  stop 'tke'
3479  endif
3480 
3481  ierr=NF_INQ_VARID(nid,"ugeo",var3didin(7))
3482  if(ierr/=NF_NOERR) then
3483  write(*,*) NF_STRERROR(ierr)
3484  stop 'ug'
3485  endif
3486 
3487  ierr=NF_INQ_VARID(nid,"vgeo",var3didin(8))
3488  if(ierr/=NF_NOERR) then
3489  write(*,*) NF_STRERROR(ierr)
3490  stop 'vg'
3491  endif
3492 
3493  ierr=NF_INQ_VARID(nid,"wls",var3didin(9))
3494  if(ierr/=NF_NOERR) then
3495  write(*,*) NF_STRERROR(ierr)
3496  stop 'wls'
3497  endif
3498 
3499  ierr=NF_INQ_VARID(nid,"dqtdx",var3didin(10))
3500  if(ierr/=NF_NOERR) then
3501  write(*,*) NF_STRERROR(ierr)
3502  stop 'dqtdx'
3503  endif
3504 
3505  ierr=NF_INQ_VARID(nid,"dqtdy",var3didin(11))
3506  if(ierr/=NF_NOERR) then
3507  write(*,*) NF_STRERROR(ierr)
3508  stop 'dqtdy'
3509  endif
3510 
3511  ierr=NF_INQ_VARID(nid,"dqtdt",var3didin(12))
3512  if(ierr/=NF_NOERR) then
3513  write(*,*) NF_STRERROR(ierr)
3514  stop 'dqtdt'
3515  endif
3516 
3517  ierr=NF_INQ_VARID(nid,"thl_rad",var3didin(13))
3518  if(ierr/=NF_NOERR) then
3519  write(*,*) NF_STRERROR(ierr)
3520  stop 'thl_rad'
3521  endif
3522 !dimensions lecture
3523 ! call catchaxis(nid,ntime,nlevel,time,z,ierr)
3524 
3525 #ifdef NC_DOUBLE
3526  ierr = NF_GET_VAR_DOUBLE(nid,var3didin(1),zz)
3527 #else
3528  ierr = NF_GET_VAR_REAL(nid,var3didin(1),zz)
3529 #endif
3530  if(ierr/=NF_NOERR) then
3531  write(*,*) NF_STRERROR(ierr)
3532  stop "getvarup"
3533  endif
3534 ! write(*,*)'lecture z ok',zz
3535 
3536 #ifdef NC_DOUBLE
3537  ierr = NF_GET_VAR_DOUBLE(nid,var3didin(2),thl)
3538 #else
3539  ierr = NF_GET_VAR_REAL(nid,var3didin(2),thl)
3540 #endif
3541  if(ierr/=NF_NOERR) then
3542  write(*,*) NF_STRERROR(ierr)
3543  stop "getvarup"
3544  endif
3545 ! write(*,*)'lecture thl ok',thl
3546 
3547 #ifdef NC_DOUBLE
3548  ierr = NF_GET_VAR_DOUBLE(nid,var3didin(3),qt)
3549 #else
3550  ierr = NF_GET_VAR_REAL(nid,var3didin(3),qt)
3551 #endif
3552  if(ierr/=NF_NOERR) then
3553  write(*,*) NF_STRERROR(ierr)
3554  stop "getvarup"
3555  endif
3556 ! write(*,*)'lecture qt ok',qt
3557 
3558 #ifdef NC_DOUBLE
3559  ierr = NF_GET_VAR_DOUBLE(nid,var3didin(4),u)
3560 #else
3561  ierr = NF_GET_VAR_REAL(nid,var3didin(4),u)
3562 #endif
3563  if(ierr/=NF_NOERR) then
3564  write(*,*) NF_STRERROR(ierr)
3565  stop "getvarup"
3566  endif
3567 ! write(*,*)'lecture u ok',u
3568 
3569 #ifdef NC_DOUBLE
3570  ierr = NF_GET_VAR_DOUBLE(nid,var3didin(5),v)
3571 #else
3572  ierr = NF_GET_VAR_REAL(nid,var3didin(5),v)
3573 #endif
3574  if(ierr/=NF_NOERR) then
3575  write(*,*) NF_STRERROR(ierr)
3576  stop "getvarup"
3577  endif
3578 ! write(*,*)'lecture v ok',v
3579 
3580 #ifdef NC_DOUBLE
3581  ierr = NF_GET_VAR_DOUBLE(nid,var3didin(6),tke)
3582 #else
3583  ierr = NF_GET_VAR_REAL(nid,var3didin(6),tke)
3584 #endif
3585  if(ierr/=NF_NOERR) then
3586  write(*,*) NF_STRERROR(ierr)
3587  stop "getvarup"
3588  endif
3589 ! write(*,*)'lecture tke ok',tke
3590 
3591 #ifdef NC_DOUBLE
3592  ierr = NF_GET_VAR_DOUBLE(nid,var3didin(7),ug)
3593 #else
3594  ierr = NF_GET_VAR_REAL(nid,var3didin(7),ug)
3595 #endif
3596  if(ierr/=NF_NOERR) then
3597  write(*,*) NF_STRERROR(ierr)
3598  stop "getvarup"
3599  endif
3600 ! write(*,*)'lecture ug ok',ug
3601 
3602 #ifdef NC_DOUBLE
3603  ierr = NF_GET_VAR_DOUBLE(nid,var3didin(8),vg)
3604 #else
3605  ierr = NF_GET_VAR_REAL(nid,var3didin(8),vg)
3606 #endif
3607  if(ierr/=NF_NOERR) then
3608  write(*,*) NF_STRERROR(ierr)
3609  stop "getvarup"
3610  endif
3611 ! write(*,*)'lecture vg ok',vg
3612 
3613 #ifdef NC_DOUBLE
3614  ierr = NF_GET_VAR_DOUBLE(nid,var3didin(9),wls)
3615 #else
3616  ierr = NF_GET_VAR_REAL(nid,var3didin(9),wls)
3617 #endif
3618  if(ierr/=NF_NOERR) then
3619  write(*,*) NF_STRERROR(ierr)
3620  stop "getvarup"
3621  endif
3622 ! write(*,*)'lecture wls ok',wls
3623 
3624 #ifdef NC_DOUBLE
3625  ierr = NF_GET_VAR_DOUBLE(nid,var3didin(10),dqtdx)
3626 #else
3627  ierr = NF_GET_VAR_REAL(nid,var3didin(10),dqtdx)
3628 #endif
3629  if(ierr/=NF_NOERR) then
3630  write(*,*) NF_STRERROR(ierr)
3631  stop "getvarup"
3632  endif
3633 ! write(*,*)'lecture dqtdx ok',dqtdx
3634 
3635 #ifdef NC_DOUBLE
3636  ierr = NF_GET_VAR_DOUBLE(nid,var3didin(11),dqtdy)
3637 #else
3638  ierr = NF_GET_VAR_REAL(nid,var3didin(11),dqtdy)
3639 #endif
3640  if(ierr/=NF_NOERR) then
3641  write(*,*) NF_STRERROR(ierr)
3642  stop "getvarup"
3643  endif
3644 ! write(*,*)'lecture dqtdy ok',dqtdy
3645 
3646 #ifdef NC_DOUBLE
3647  ierr = NF_GET_VAR_DOUBLE(nid,var3didin(12),dqtdt)
3648 #else
3649  ierr = NF_GET_VAR_REAL(nid,var3didin(12),dqtdt)
3650 #endif
3651  if(ierr/=NF_NOERR) then
3652  write(*,*) NF_STRERROR(ierr)
3653  stop "getvarup"
3654  endif
3655 ! write(*,*)'lecture dqtdt ok',dqtdt
3656 
3657 #ifdef NC_DOUBLE
3658  ierr = NF_GET_VAR_DOUBLE(nid,var3didin(13),thl_rad)
3659 #else
3660  ierr = NF_GET_VAR_REAL(nid,var3didin(13),thl_rad)
3661 #endif
3662  if(ierr/=NF_NOERR) then
3663  write(*,*) NF_STRERROR(ierr)
3664  stop "getvarup"
3665  endif
3666 ! write(*,*)'lecture thl_rad ok',thl_rad
3667 
3668  return
3669  end subroutine read_fire
3670 !=====================================================================
3671  subroutine read_dice(fich_dice,nlevel,ntime &
3672  & ,zz,pres,th,qv,u,v,o3 &
3673  & ,shf,lhf,lwup,swup,tg,ustar,psurf,ug,vg &
3674  & ,hadvt,hadvq,hadvu,hadvv,w,omega)
3675 
3676 !program reading initial profils and forcings of the Dice case study
3677 
3678 
3679  implicit none
3680 
3681 #include "netcdf.inc"
3682 
3683  integer ntime,nlevel
3684  integer l,k
3685  character*80 :: fich_dice
3686  real*8 time(ntime)
3687  real*8 zz(nlevel)
3688 
3689  real*8 th(nlevel),pres(nlevel)
3690  real*8 qv(nlevel),u(nlevel),v(nlevel),o3(nlevel)
3691  real*8 shf(ntime),lhf(ntime),lwup(ntime),swup(ntime),tg(ntime)
3692  real*8 ustar(ntime),psurf(ntime),ug(ntime),vg(ntime)
3693  real*8 hadvt(nlevel,ntime),hadvq(nlevel,ntime),hadvu(nlevel,ntime)
3694  real*8 hadvv(nlevel,ntime),w(nlevel,ntime),omega(nlevel,ntime)
3695 
3696  integer nid, ierr
3697  integer nbvar3d
3698  parameter(nbvar3d=30)
3699  integer var3didin(nbvar3d)
3700 
3701  ierr = NF_OPEN(fich_dice,NF_NOWRITE,nid)
3702  if (ierr.NE.NF_NOERR) then
3703  write(*,*) 'ERROR: Pb opening forcings nc file '
3704  write(*,*) NF_STRERROR(ierr)
3705  stop ""
3706  endif
3707 
3708 
3709  ierr=NF_INQ_VARID(nid,"height",var3didin(1))
3710  if(ierr/=NF_NOERR) then
3711  write(*,*) NF_STRERROR(ierr)
3712  stop 'height'
3713  endif
3714 
3715  ierr=NF_INQ_VARID(nid,"pf",var3didin(11))
3716  if(ierr/=NF_NOERR) then
3717  write(*,*) NF_STRERROR(ierr)
3718  stop 'pf'
3719  endif
3720 
3721  ierr=NF_INQ_VARID(nid,"theta",var3didin(12))
3722  if(ierr/=NF_NOERR) then
3723  write(*,*) NF_STRERROR(ierr)
3724  stop 'theta'
3725  endif
3726 
3727  ierr=NF_INQ_VARID(nid,"qv",var3didin(13))
3728  if(ierr/=NF_NOERR) then
3729  write(*,*) NF_STRERROR(ierr)
3730  stop 'qv'
3731  endif
3732 
3733  ierr=NF_INQ_VARID(nid,"u",var3didin(14))
3734  if(ierr/=NF_NOERR) then
3735  write(*,*) NF_STRERROR(ierr)
3736  stop 'u'
3737  endif
3738 
3739  ierr=NF_INQ_VARID(nid,"v",var3didin(15))
3740  if(ierr/=NF_NOERR) then
3741  write(*,*) NF_STRERROR(ierr)
3742  stop 'v'
3743  endif
3744 
3745  ierr=NF_INQ_VARID(nid,"o3mmr",var3didin(16))
3746  if(ierr/=NF_NOERR) then
3747  write(*,*) NF_STRERROR(ierr)
3748  stop 'o3'
3749  endif
3750 
3751  ierr=NF_INQ_VARID(nid,"shf",var3didin(2))
3752  if(ierr/=NF_NOERR) then
3753  write(*,*) NF_STRERROR(ierr)
3754  stop 'shf'
3755  endif
3756 
3757  ierr=NF_INQ_VARID(nid,"lhf",var3didin(3))
3758  if(ierr/=NF_NOERR) then
3759  write(*,*) NF_STRERROR(ierr)
3760  stop 'lhf'
3761  endif
3762 
3763  ierr=NF_INQ_VARID(nid,"lwup",var3didin(4))
3764  if(ierr/=NF_NOERR) then
3765  write(*,*) NF_STRERROR(ierr)
3766  stop 'lwup'
3767  endif
3768 
3769  ierr=NF_INQ_VARID(nid,"swup",var3didin(5))
3770  if(ierr/=NF_NOERR) then
3771  write(*,*) NF_STRERROR(ierr)
3772  stop 'dqtdx'
3773  endif
3774 
3775  ierr=NF_INQ_VARID(nid,"Tg",var3didin(6))
3776  if(ierr/=NF_NOERR) then
3777  write(*,*) NF_STRERROR(ierr)
3778  stop 'Tg'
3779  endif
3780 
3781  ierr=NF_INQ_VARID(nid,"ustar",var3didin(7))
3782  if(ierr/=NF_NOERR) then
3783  write(*,*) NF_STRERROR(ierr)
3784  stop 'ustar'
3785  endif
3786 
3787  ierr=NF_INQ_VARID(nid,"psurf",var3didin(8))
3788  if(ierr/=NF_NOERR) then
3789  write(*,*) NF_STRERROR(ierr)
3790  stop 'psurf'
3791  endif
3792 
3793  ierr=NF_INQ_VARID(nid,"Ug",var3didin(9))
3794  if(ierr/=NF_NOERR) then
3795  write(*,*) NF_STRERROR(ierr)
3796  stop 'Ug'
3797  endif
3798 
3799  ierr=NF_INQ_VARID(nid,"Vg",var3didin(10))
3800  if(ierr/=NF_NOERR) then
3801  write(*,*) NF_STRERROR(ierr)
3802  stop 'Vg'
3803  endif
3804 
3805  ierr=NF_INQ_VARID(nid,"hadvT",var3didin(17))
3806  if(ierr/=NF_NOERR) then
3807  write(*,*) NF_STRERROR(ierr)
3808  stop 'hadvT'
3809  endif
3810 
3811  ierr=NF_INQ_VARID(nid,"hadvq",var3didin(18))
3812  if(ierr/=NF_NOERR) then
3813  write(*,*) NF_STRERROR(ierr)
3814  stop 'hadvq'
3815  endif
3816 
3817  ierr=NF_INQ_VARID(nid,"hadvu",var3didin(19))
3818  if(ierr/=NF_NOERR) then
3819  write(*,*) NF_STRERROR(ierr)
3820  stop 'hadvu'
3821  endif
3822 
3823  ierr=NF_INQ_VARID(nid,"hadvv",var3didin(20))
3824  if(ierr/=NF_NOERR) then
3825  write(*,*) NF_STRERROR(ierr)
3826  stop 'hadvv'
3827  endif
3828 
3829  ierr=NF_INQ_VARID(nid,"w",var3didin(21))
3830  if(ierr/=NF_NOERR) then
3831  write(*,*) NF_STRERROR(ierr)
3832  stop 'w'
3833  endif
3834 
3835  ierr=NF_INQ_VARID(nid,"omega",var3didin(22))
3836  if(ierr/=NF_NOERR) then
3837  write(*,*) NF_STRERROR(ierr)
3838  stop 'omega'
3839  endif
3840 !dimensions lecture
3841 ! call catchaxis(nid,ntime,nlevel,time,z,ierr)
3842 
3843 #ifdef NC_DOUBLE
3844  ierr = NF_GET_VAR_DOUBLE(nid,var3didin(1),zz)
3845 #else
3846  ierr = NF_GET_VAR_REAL(nid,var3didin(1),zz)
3847 #endif
3848  if(ierr/=NF_NOERR) then
3849  write(*,*) NF_STRERROR(ierr)
3850  stop "getvarup"
3851  endif
3852 ! write(*,*)'lecture zz ok',zz
3853 
3854 #ifdef NC_DOUBLE
3855  ierr = NF_GET_VAR_DOUBLE(nid,var3didin(11),pres)
3856 #else
3857  ierr = NF_GET_VAR_REAL(nid,var3didin(11),pres)
3858 #endif
3859  if(ierr/=NF_NOERR) then
3860  write(*,*) NF_STRERROR(ierr)
3861  stop "getvarup"
3862  endif
3863 ! write(*,*)'lecture pres ok',pres
3864 
3865 #ifdef NC_DOUBLE
3866  ierr = NF_GET_VAR_DOUBLE(nid,var3didin(12),th)
3867 #else
3868  ierr = NF_GET_VAR_REAL(nid,var3didin(12),th)
3869 #endif
3870  if(ierr/=NF_NOERR) then
3871  write(*,*) NF_STRERROR(ierr)
3872  stop "getvarup"
3873  endif
3874 ! write(*,*)'lecture th ok',th
3875 
3876 #ifdef NC_DOUBLE
3877  ierr = NF_GET_VAR_DOUBLE(nid,var3didin(13),qv)
3878 #else
3879  ierr = NF_GET_VAR_REAL(nid,var3didin(13),qv)
3880 #endif
3881  if(ierr/=NF_NOERR) then
3882  write(*,*) NF_STRERROR(ierr)
3883  stop "getvarup"
3884  endif
3885 ! write(*,*)'lecture qv ok',qv
3886 
3887 #ifdef NC_DOUBLE
3888  ierr = NF_GET_VAR_DOUBLE(nid,var3didin(14),u)
3889 #else
3890  ierr = NF_GET_VAR_REAL(nid,var3didin(14),u)
3891 #endif
3892  if(ierr/=NF_NOERR) then
3893  write(*,*) NF_STRERROR(ierr)
3894  stop "getvarup"
3895  endif
3896 ! write(*,*)'lecture u ok',u
3897 
3898 #ifdef NC_DOUBLE
3899  ierr = NF_GET_VAR_DOUBLE(nid,var3didin(15),v)
3900 #else
3901  ierr = NF_GET_VAR_REAL(nid,var3didin(15),v)
3902 #endif
3903  if(ierr/=NF_NOERR) then
3904  write(*,*) NF_STRERROR(ierr)
3905  stop "getvarup"
3906  endif
3907 ! write(*,*)'lecture v ok',v
3908 
3909 #ifdef NC_DOUBLE
3910  ierr = NF_GET_VAR_DOUBLE(nid,var3didin(16),o3)
3911 #else
3912  ierr = NF_GET_VAR_REAL(nid,var3didin(16),o3)
3913 #endif
3914  if(ierr/=NF_NOERR) then
3915  write(*,*) NF_STRERROR(ierr)
3916  stop "getvarup"
3917  endif
3918 ! write(*,*)'lecture o3 ok',o3
3919 
3920 #ifdef NC_DOUBLE
3921  ierr = NF_GET_VAR_DOUBLE(nid,var3didin(2),shf)
3922 #else
3923  ierr = NF_GET_VAR_REAL(nid,var3didin(2),shf)
3924 #endif
3925  if(ierr/=NF_NOERR) then
3926  write(*,*) NF_STRERROR(ierr)
3927  stop "getvarup"
3928  endif
3929 ! write(*,*)'lecture shf ok',shf
3930 
3931 #ifdef NC_DOUBLE
3932  ierr = NF_GET_VAR_DOUBLE(nid,var3didin(3),lhf)
3933 #else
3934  ierr = NF_GET_VAR_REAL(nid,var3didin(3),lhf)
3935 #endif
3936  if(ierr/=NF_NOERR) then
3937  write(*,*) NF_STRERROR(ierr)
3938  stop "getvarup"
3939  endif
3940 ! write(*,*)'lecture lhf ok',lhf
3941 
3942 #ifdef NC_DOUBLE
3943  ierr = NF_GET_VAR_DOUBLE(nid,var3didin(4),lwup)
3944 #else
3945  ierr = NF_GET_VAR_REAL(nid,var3didin(4),lwup)
3946 #endif
3947  if(ierr/=NF_NOERR) then
3948  write(*,*) NF_STRERROR(ierr)
3949  stop "getvarup"
3950  endif
3951 ! write(*,*)'lecture lwup ok',lwup
3952 
3953 #ifdef NC_DOUBLE
3954  ierr = NF_GET_VAR_DOUBLE(nid,var3didin(5),swup)
3955 #else
3956  ierr = NF_GET_VAR_REAL(nid,var3didin(5),swup)
3957 #endif
3958  if(ierr/=NF_NOERR) then
3959  write(*,*) NF_STRERROR(ierr)
3960  stop "getvarup"
3961  endif
3962 ! write(*,*)'lecture swup ok',swup
3963 
3964 #ifdef NC_DOUBLE
3965  ierr = NF_GET_VAR_DOUBLE(nid,var3didin(6),tg)
3966 #else
3967  ierr = NF_GET_VAR_REAL(nid,var3didin(6),tg)
3968 #endif
3969  if(ierr/=NF_NOERR) then
3970  write(*,*) NF_STRERROR(ierr)
3971  stop "getvarup"
3972  endif
3973 ! write(*,*)'lecture tg ok',tg
3974 
3975 #ifdef NC_DOUBLE
3976  ierr = NF_GET_VAR_DOUBLE(nid,var3didin(7),ustar)
3977 #else
3978  ierr = NF_GET_VAR_REAL(nid,var3didin(7),ustar)
3979 #endif
3980  if(ierr/=NF_NOERR) then
3981  write(*,*) NF_STRERROR(ierr)
3982  stop "getvarup"
3983  endif
3984 ! write(*,*)'lecture ustar ok',ustar
3985 
3986 #ifdef NC_DOUBLE
3987  ierr = NF_GET_VAR_DOUBLE(nid,var3didin(8),psurf)
3988 #else
3989  ierr = NF_GET_VAR_REAL(nid,var3didin(8),psurf)
3990 #endif
3991  if(ierr/=NF_NOERR) then
3992  write(*,*) NF_STRERROR(ierr)
3993  stop "getvarup"
3994  endif
3995 ! write(*,*)'lecture psurf ok',psurf
3996 
3997 #ifdef NC_DOUBLE
3998  ierr = NF_GET_VAR_DOUBLE(nid,var3didin(9),ug)
3999 #else
4000  ierr = NF_GET_VAR_REAL(nid,var3didin(9),ug)
4001 #endif
4002  if(ierr/=NF_NOERR) then
4003  write(*,*) NF_STRERROR(ierr)
4004  stop "getvarup"
4005  endif
4006 ! write(*,*)'lecture ug ok',ug
4007 
4008 #ifdef NC_DOUBLE
4009  ierr = NF_GET_VAR_DOUBLE(nid,var3didin(10),vg)
4010 #else
4011  ierr = NF_GET_VAR_REAL(nid,var3didin(10),vg)
4012 #endif
4013  if(ierr/=NF_NOERR) then
4014  write(*,*) NF_STRERROR(ierr)
4015  stop "getvarup"
4016  endif
4017 ! write(*,*)'lecture vg ok',vg
4018 
4019 #ifdef NC_DOUBLE
4020  ierr = NF_GET_VAR_DOUBLE(nid,var3didin(17),hadvt)
4021 #else
4022  ierr = NF_GET_VAR_REAL(nid,var3didin(17),hadvt)
4023 #endif
4024  if(ierr/=NF_NOERR) then
4025  write(*,*) NF_STRERROR(ierr)
4026  stop "getvarup"
4027  endif
4028 ! write(*,*)'lecture hadvt ok',hadvt
4029 
4030 #ifdef NC_DOUBLE
4031  ierr = NF_GET_VAR_DOUBLE(nid,var3didin(18),hadvq)
4032 #else
4033  ierr = NF_GET_VAR_REAL(nid,var3didin(18),hadvq)
4034 #endif
4035  if(ierr/=NF_NOERR) then
4036  write(*,*) NF_STRERROR(ierr)
4037  stop "getvarup"
4038  endif
4039 ! write(*,*)'lecture hadvq ok',hadvq
4040 
4041 #ifdef NC_DOUBLE
4042  ierr = NF_GET_VAR_DOUBLE(nid,var3didin(19),hadvu)
4043 #else
4044  ierr = NF_GET_VAR_REAL(nid,var3didin(19),hadvu)
4045 #endif
4046  if(ierr/=NF_NOERR) then
4047  write(*,*) NF_STRERROR(ierr)
4048  stop "getvarup"
4049  endif
4050 ! write(*,*)'lecture hadvu ok',hadvu
4051 
4052 #ifdef NC_DOUBLE
4053  ierr = NF_GET_VAR_DOUBLE(nid,var3didin(20),hadvv)
4054 #else
4055  ierr = NF_GET_VAR_REAL(nid,var3didin(20),hadvv)
4056 #endif
4057  if(ierr/=NF_NOERR) then
4058  write(*,*) NF_STRERROR(ierr)
4059  stop "getvarup"
4060  endif
4061 ! write(*,*)'lecture hadvv ok',hadvv
4062 
4063 #ifdef NC_DOUBLE
4064  ierr = NF_GET_VAR_DOUBLE(nid,var3didin(21),w)
4065 #else
4066  ierr = NF_GET_VAR_REAL(nid,var3didin(21),w)
4067 #endif
4068  if(ierr/=NF_NOERR) then
4069  write(*,*) NF_STRERROR(ierr)
4070  stop "getvarup"
4071  endif
4072 ! write(*,*)'lecture w ok',w
4073 
4074 #ifdef NC_DOUBLE
4075  ierr = NF_GET_VAR_DOUBLE(nid,var3didin(22),omega)
4076 #else
4077  ierr = NF_GET_VAR_REAL(nid,var3didin(22),omega)
4078 #endif
4079  if(ierr/=NF_NOERR) then
4080  write(*,*) NF_STRERROR(ierr)
4081  stop "getvarup"
4082  endif
4083 ! write(*,*)'lecture omega ok',omega
4084 
4085  return
4086  end subroutine read_dice
4087 !=====================================================================
4088 ! Reads CIRC input files
4089 
4090  SUBROUTINE read_circ(nlev_circ,cf,lwp,iwp,reliq,reice,t,z,p,pm,h2o,o3,sza)
4091 
4092  parameter (ncm_1=49180)
4093 #include "YOMCST.h"
4094 
4095  real albsfc(ncm_1), albsfc_w(ncm_1)
4096  real cf(nlev_circ), icefra(nlev_circ), deice(nlev_circ), &
4097  reliq(nlev_circ), reice(nlev_circ), lwp(nlev_circ), iwp(nlev_circ)
4098  real t(nlev_circ+1), z(nlev_circ+1), dz(nlev_circ), p(nlev_circ+1)
4099  real aer_beta(nlev_circ), waer(nlev_circ), gaer(nlev_circ)
4100  real pm(nlev_circ), tm(nlev_circ), h2o(nlev_circ), o3(nlev_circ)
4101  real co2(nlev_circ), n2o(nlev_circ), co(nlev_circ), ch4(nlev_circ), &
4102  o2(nlev_circ), ccl4(nlev_circ), f11(nlev_circ), f12(nlev_circ)
4103 ! za= zenital angle
4104 ! sza= cosinus angle zenital
4105  real wavn(ncm_1), ssf(ncm_1),za,sza
4106  integer nlev
4107 
4108 
4109 ! Open the files
4110 
4111  open (11, file='Tsfc_sza_nlev_case.txt', status='old')
4112  open (12, file='level_input_case.txt', status='old')
4113  open (13, file='layer_input_case.txt', status='old')
4114  open (14, file='aerosol_input_case.txt', status='old')
4115  open (15, file='cloud_input_case.txt', status='old')
4116  open (16, file='sfcalbedo_input_case.txt', status='old')
4117 
4118 ! Read scalar information
4119  do iskip=1,5
4120  read (11, *)
4121  enddo
4122  read (11, '(i8)') nlev
4123  read (11, '(f10.2)') tsfc
4124  read (11, '(f10.2)') za
4125  read (11, '(f10.4)') sw_dn_toa
4126  sza=cos(za/180.*RPI)
4127  print *,'nlev,tsfc,sza,sw_dn_toa,RPI',nlev,tsfc,sza,sw_dn_toa,RPI
4128  close(11)
4129 
4130 ! Read level information
4131  read (12, *)
4132  do il=1,nlev
4133  read (12, 302) ilev, z(il), p(il), t(il)
4134  z(il)=z(il)*1000. ! z donne en km
4135  p(il)=p(il)*100. ! p donne en mb
4136  enddo
4137 302 format (i8, f8.3, 2f9.2)
4138  close(12)
4139 
4140 ! Read layer information (midpoint values)
4141  do iskip=1,3
4142  read (13, *)
4143  enddo
4144  do il=1,nlev-1
4145  read (13, 303) ilev,pm(il),tm(il),h2o(il),co2(il),o3(il), &
4146  n2o(il),co(il),ch4(il),o2(il),ccl4(il), &
4147  f11(il),f12(il)
4148  pm(il)=pm(il)*100.
4149  enddo
4150 303 format (i8, 2f9.2, 10(2x,e13.7))
4151  close(13)
4152 
4153 ! Read aerosol layer information
4154  do iskip=1,3
4155  read (14, *)
4156  enddo
4157  read (14, '(f10.2)') aer_alpha
4158  read (14, *)
4159  read (14, *)
4160  do il=1,nlev-1
4161  read (14, 304) ilev, aer_beta(il), waer(il), gaer(il)
4162  enddo
4163 304 format (i8, f9.5, 2f8.3)
4164  close(14)
4165 
4166 ! Read cloud information
4167  do iskip=1,3
4168  read (15, *)
4169  enddo
4170  do il=1,nlev-1
4171  read (15, 305) ilev, cf(il), lwp(il), iwp(il), reliq(il), reice(il)
4172  lwp(il)=lwp(il)/1000. ! lwp donne en g/kg
4173  iwp(il)=iwp(il)/1000. ! iwp donne en g/kg
4174  reliq(il)=reliq(il)/1000000. ! reliq donne en microns
4175  reice(il)=reice(il)/1000000. ! reice donne en microns
4176  enddo
4177 305 format (i8, f8.3, 4f9.2)
4178  close(15)
4179 
4180 ! Read surface albedo (weighted & unweighted) and spectral solar irradiance
4181  do iskip=1,6
4182  read (16, *)
4183  enddo
4184  do icm_1=1,ncm_1
4185  read (16, 306) wavn(icm_1), albsfc(icm_1), albsfc_w(icm_1), ssf(icm_1)
4186  enddo
4187 306 format(f10.1, 2f12.5, f14.8)
4188  close(16)
4189 
4190  return
4191  end subroutine read_circ
4192 !=====================================================================
4193 ! Reads RTMIP input files
4194 
4195  SUBROUTINE read_rtmip(nlev_rtmip,play,plev,t,h2o,o3)
4196 
4197 #include "YOMCST.h"
4198 
4199  real t(nlev_rtmip), pt(nlev_rtmip),pb(nlev_rtmip),h2o(nlev_rtmip), o3(nlev_rtmip)
4200  real temp(nlev_rtmip), play(nlev_rtmip),ovap(nlev_rtmip), oz(nlev_rtmip),plev(nlev_rtmip+1)
4201  integer nlev
4202 
4203 
4204 ! Open the files
4205 
4206  open (11, file='low_resolution_profile.txt', status='old')
4207 
4208 ! Read level information
4209  read (11, *)
4210  do il=1,nlev_rtmip
4211  read (11, 302) pt(il), pb(il), t(il),h2o(il),o3(il)
4212  enddo
4213  do il=1,nlev_rtmip
4214  play(il)=pt(nlev_rtmip-il+1)*100. ! p donne en mb
4215  temp(il)=t(nlev_rtmip-il+1)
4216  ovap(il)=h2o(nlev_rtmip-il+1)
4217  oz(il)=o3(nlev_rtmip-il+1)
4218  enddo
4219  do il=1,39
4220  plev(il)=play(il)+(play(il+1)-play(il))/2.
4221  print *,'il p t ovap oz=',il,plev(il),temp(il),ovap(il),oz(il)
4222  enddo
4223  plev(41)=101300.
4224 302 format (e16.10,3x,e16.10,3x,e16.10,3x,e12.6,3x,e12.6)
4225  close(12)
4226 
4227  return
4228  end subroutine read_rtmip
4229 !=====================================================================
4230 
4231 ! Subroutines for nudging
4232 
4233  Subroutine Nudge_RHT_init (paprs,pplay,t,q,t_targ,rh_targ)
4234 ! ========================================================
4235  USE dimphy
4236 
4237  implicit none
4238 
4239 ! ========================================================
4240  REAL paprs(klon,klevp1)
4241  REAL pplay(klon,klev)
4242 !
4243 ! Variables d'etat
4244  REAL t(klon,klev)
4245  REAL q(klon,klev)
4246 !
4247 ! Profiles cible
4248  REAL t_targ(klon,klev)
4249  REAL rh_targ(klon,klev)
4250 !
4251  INTEGER k,i
4252  REAL zx_qs
4253 
4254 ! Declaration des constantes et des fonctions thermodynamiques
4255 !
4256 include "YOMCST.h"
4257 include "YOETHF.h"
4258 !
4259 ! ----------------------------------------
4260 ! Statement functions
4261 include "FCTTRE.h"
4262 ! ----------------------------------------
4263 !
4264  DO k = 1,klev
4265  DO i = 1,klon
4266  t_targ(i,k) = t(i,k)
4267  IF (t(i,k).LT.RTT) THEN
4268  zx_qs = qsats(t(i,k))/(pplay(i,k))
4269  ELSE
4270  zx_qs = qsatl(t(i,k))/(pplay(i,k))
4271  ENDIF
4272  rh_targ(i,k) = q(i,k)/zx_qs
4273  ENDDO
4274  ENDDO
4275  print *, 't_targ',t_targ
4276  print *, 'rh_targ',rh_targ
4277 !
4278 !
4279  RETURN
4280  END
4281 
4282  Subroutine Nudge_UV_init (paprs,pplay,u,v,u_targ,v_targ)
4283 ! ========================================================
4284  USE dimphy
4285 
4286  implicit none
4287 
4288 ! ========================================================
4289  REAL paprs(klon,klevp1)
4290  REAL pplay(klon,klev)
4291 !
4292 ! Variables d'etat
4293  REAL u(klon,klev)
4294  REAL v(klon,klev)
4295 !
4296 ! Profiles cible
4297  REAL u_targ(klon,klev)
4298  REAL v_targ(klon,klev)
4299 !
4300  INTEGER k,i
4301 !
4302  DO k = 1,klev
4303  DO i = 1,klon
4304  u_targ(i,k) = u(i,k)
4305  v_targ(i,k) = v(i,k)
4306  ENDDO
4307  ENDDO
4308  print *, 'u_targ',u_targ
4309  print *, 'v_targ',v_targ
4310 !
4311 !
4312  RETURN
4313  END
4314 
4315  Subroutine Nudge_RHT (dtime,paprs,pplay,t_targ,rh_targ,t,q, &
4316  & d_t,d_q)
4317 ! ========================================================
4318  USE dimphy
4319 
4320  implicit none
4321 
4322 ! ========================================================
4323  REAL dtime
4324  REAL paprs(klon,klevp1)
4325  REAL pplay(klon,klev)
4326 !
4327 ! Variables d'etat
4328  REAL t(klon,klev)
4329  REAL q(klon,klev)
4330 !
4331 ! Tendances
4332  REAL d_t(klon,klev)
4333  REAL d_q(klon,klev)
4334 !
4335 ! Profiles cible
4336  REAL t_targ(klon,klev)
4337  REAL rh_targ(klon,klev)
4338 !
4339 ! Temps de relaxation
4340  REAL tau
4341 !c DATA tau /3600./
4342 !! DATA tau /5400./
4343  DATA tau /1800./
4344 !
4345  INTEGER k,i
4346  REAL zx_qs, rh, tnew, d_rh
4347 
4348 ! Declaration des constantes et des fonctions thermodynamiques
4349 !
4350 include "YOMCST.h"
4351 include "YOETHF.h"
4352 !
4353 ! ----------------------------------------
4354 ! Statement functions
4355 include "FCTTRE.h"
4356 ! ----------------------------------------
4357 !
4358  print *,'dtime, tau ',dtime,tau
4359  print *, 't_targ',t_targ
4360  print *, 'rh_targ',rh_targ
4361  print *,'temp ',t
4362  print *,'hum ',q
4363  DO k = 1,klev
4364  DO i = 1,klon
4365 !! IF (paprs(i,1)-pplay(i,k) .GT. 10000.) THEN
4366  IF (t(i,k).LT.RTT) THEN
4367  zx_qs = qsats(t(i,k))/(pplay(i,k))
4368  ELSE
4369  zx_qs = qsatl(t(i,k))/(pplay(i,k))
4370  ENDIF
4371  rh = q(i,k)/zx_qs
4372 !
4373  d_t(i,k) = d_t(i,k) + 1./tau*(t_targ(i,k)-t(i,k))
4374  d_rh = 1./tau*(rh_targ(i,k)-rh)
4375 !
4376  tnew = t(i,k)+d_t(i,k)
4377  IF (tnew.LT.RTT) THEN
4378  zx_qs = qsats(tnew)/(pplay(i,k))
4379  ELSE
4380  zx_qs = qsatl(tnew)/(pplay(i,k))
4381  ENDIF
4382  d_q(i,k) = d_q(i,k) + d_rh*zx_qs
4383 !
4384  print *,' k,d_t,rh,d_rh,d_q ', &
4385  k,d_t(i,k),rh,d_rh,d_q(i,k)
4386 !! ENDIF
4387 !
4388  ENDDO
4389  ENDDO
4390 !
4391  RETURN
4392  END
4393 
4394  Subroutine Nudge_UV (dtime,paprs,pplay,u_targ,v_targ,u,v, &
4395  & d_u,d_v)
4396 ! ========================================================
4397  USE dimphy
4398 
4399  implicit none
4400 
4401 ! ========================================================
4402  REAL dtime
4403  REAL paprs(klon,klevp1)
4404  REAL pplay(klon,klev)
4405 !
4406 ! Variables d'etat
4407  REAL u(klon,klev)
4408  REAL v(klon,klev)
4409 !
4410 ! Tendances
4411  REAL d_u(klon,klev)
4412  REAL d_v(klon,klev)
4413 !
4414 ! Profiles cible
4415  REAL u_targ(klon,klev)
4416  REAL v_targ(klon,klev)
4417 !
4418 ! Temps de relaxation
4419  REAL tau
4420 !c DATA tau /3600./
4421  DATA tau /5400./
4422 !
4423  INTEGER k,i
4424 
4425 !
4426  print *,'dtime, tau ',dtime,tau
4427  print *, 'u_targ',u_targ
4428  print *, 'v_targ',v_targ
4429  print *,'zonal velocity ',u
4430  print *,'meridional velocity ',v
4431  DO k = 1,klev
4432  DO i = 1,klon
4433  IF (paprs(i,1)-pplay(i,k) .GT. 10000.) THEN
4434 !
4435  d_u(i,k) = d_u(i,k) + 1./tau*(u_targ(i,k)-u(i,k))
4436  d_v(i,k) = d_v(i,k) + 1./tau*(v_targ(i,k)-v(i,k))
4437 !
4438  print *,' k,u,d_u,v,d_v ', &
4439  k,u(i,k),d_u(i,k),v(i,k),d_v(i,k)
4440  ENDIF
4441 !
4442  ENDDO
4443  ENDDO
4444 !
4445  RETURN
4446  END
4447 
level
!$Id NSTRA real GKLIFT real GVSEC REAL GWD_RANDO_RUWMAX!Maximum Eliassen Palm flux at launch level
Definition: YOEGWD.h:12
if
!$Id Turb_fcg_gcssold if(prt_level.ge.1) then print *
klon
INTERFACE SUBROUTINE RRTM_ECRT_140GP klon
Definition: rrtm_ecrt_140gp.intfb.h:3
dzoomy
!$Header!c!c!c include serre h!c REAL dzoomy
Definition: serre.h:8
itau_dyn
!$Id && itau_dyn
Definition: temps.h:15
day_end
!$Id day_end
Definition: temps.h:15
caldyn_mod::p
real, dimension(:,:), pointer, save p
Definition: caldyn_mod.F90:6
$Id
c c $Id
Definition: ini_bilKP_ave.h:11
grossismx
!$Header!c!c!c include serre h!c REAL && grossismx
Definition: serre.h:8
nudge_tsoil
logical nudge_tsoil integer isoil_nudge real tau_soil_nudge common tsoilnudge nudge_tsoil
Definition: tsoilnudge.h:3
isoil_nudge
logical nudge_tsoil integer isoil_nudge real tau_soil_nudge common tsoilnudge isoil_nudge
Definition: tsoilnudge.h:3
call_calfis_mod::q
real, dimension(:,:,:), pointer, save q
Definition: call_calfis_mod.F90:10
g
!$Id mode_top_bound COMMON comconstr g
Definition: comconst.h:7
phys_state_var_mod::heat
real, dimension(:,:), allocatable, save heat
Definition: phys_state_var_mod.F90:300
ocean_forced_mod::read_tsurf1d
subroutine read_tsurf1d(knon, sst_out)
Definition: ocean_forced_mod.F90:277
day1
!$Id Turb_fcg_gcssold get_uvd hqturb_gcssold endif!large scale llm day1
Definition: 1D_interp_cases.h:45
ts_fcg_gcssold
!$Id ts_fcg_gcssold
Definition: 1D_interp_cases.h:9
REAL
!$Id ok_orolf LOGICAL ok_limitvrai LOGICAL ok_all_xml INTEGER iflag_ener_conserv REAL solaire REAL(kind=8) RCO2
preff
!$Id preff
Definition: comvert.h:8
l
c c $Id c c calculs statistiques distribution nuage ftion du regime dynamique c c Ce calcul doit etre fait a partir de valeurs mensuelles CALL nbregdyn DO kmaxm1 DO l
Definition: calcul_REGDYN.h:13
geometry_mod::longitude
real, dimension(:), allocatable, save longitude
Definition: geometry_mod.F90:5
R
!$Id RNAVO!A1 Astronomical constants REAL ROMEGA!A1 bis Constantes concernant l orbite de la R_incl!A1 Geoide REAL R1SA!A1 Radiation!REAL RI0 REAL RSIGMA!A1 Thermodynamic gas phase REAL R
Definition: YOMCST.h:11
plev_prof
!$Id Turb_fcg_gcssold get_uvd hqturb_gcssold endif!large scale llm day day1 day day1 *dt_toga endif!time annee_ref dt_toga u_toga vq_toga plev_prof
Definition: 1D_interp_cases.h:45
llmp1
!$Header llmp1
Definition: paramet.h:14
dtime
!$Header!integer nvarmx dtime
Definition: gradsdef.h:20
kappa
!$Id mode_top_bound COMMON comconstr kappa
Definition: comconst.h:7
tend_q
!$Id tend_q
Definition: compar1d.h:5
delta
!$Id!Parameters for nlm real spfac!IM cf epmax real ptcrit real omtrain real dttrig real alpha real delta real betad COMMON cv30param nlm spfac &!IM cf ptcrit omtrain dttrig delta
Definition: cv30param.h:5
imp_fcg_gcssold
!$Id && imp_fcg_gcssold
Definition: 1D_interp_cases.h:9
iostart::length
integer, parameter length
Definition: iostart.F90:8
control_mod::dayref
integer, save dayref
Definition: control_mod.F90:26
Tsoil_nudge
logical nudge_tsoil integer isoil_nudge real Tsoil_nudge
Definition: tsoilnudge.h:3
iostart::open_startphy
subroutine, public open_startphy(filename)
Definition: iostart.F90:32
clon
!$Header!c!c!c include serre h!c REAL clon
Definition: serre.h:8
FALSE
!$Id ***************************************!ECRITURE DU FALSE
Definition: write_histrac.h:9
nudge_t
!$Id tend_rayo nudge_t
Definition: compar1d.h:5
im
!$Id mode_top_bound COMMON comconstr omeg dissip_zref ihf INTEGER im
Definition: comconst.h:7
jm
!$Id jm
Definition: comconst.h:7
beta
!$Id!Parameters for nlm real spfac!IM cf epmax real ptcrit real omtrain real dttrig real beta
Definition: cv30param.h:5
tau
!$Id!Parameters for nlm real spfac!IM cf epmax real ptcrit real omtrain real dttrig real tau
Definition: cv30param.h:5
dt
!$Id calend INTEGER itaufin INTEGER itau_phy INTEGER day_ref REAL dt
Definition: temps.h:15
phys_state_var_mod::cool
real, dimension(:,:), allocatable, save cool
Definition: phys_state_var_mod.F90:304
abort_gcm
subroutine abort_gcm(modname, message, ierr)
Definition: abort_gcm.F:7
llm
!$Id Turb_fcg_gcssold get_uvd hqturb_gcssold endif!large scale llm day day1 day day1 *dt_toga endif!time annee_ref dt_toga u_toga vq_toga w_prof vq_prof llm day day1 day day1 *dt_dice endif!time annee_ref dt_dice swup_dice vg_dice omega_dice tg_prof vg_profd w_profd omega_profd!do llm!print llm l llm
Definition: 1D_interp_cases.h:103
INTEGER
INTERFACE SUBROUTINE JPRB INTEGER(KIND=JPIM)
dimphy::klev
integer, save klev
Definition: dimphy.F90:7
airephy
!$Id!c c c Common de passage de la geometrie de la dynamique a la physique real airephy(klon)
pi
!$Id mode_top_bound COMMON comconstr && pi
Definition: comconst.h:7
forcing
!$Id Turb_fcg_gcssold get_uvd hqturb_gcssold endif!large scale forcing
Definition: 1D_interp_cases.h:16
fichier
!$Header!integer nvarmx s s s fichier
Definition: gradsdef.h:20
k
c c $Id c c calculs statistiques distribution nuage ftion du regime dynamique c c Ce calcul doit etre fait a partir de valeurs mensuelles CALL nbregdyn DO k
Definition: calcul_REGDYN.h:12
i
!$Id klon initialisation mois suivants day_rain itap ENDIF!Calcul fin de nday_rain calcul nday_rain itap DO i
Definition: calcul_divers.h:24
imposed
!surface temperature imposed
Definition: 1DUTILS.h:81
infotrac::nqtot
integer, save nqtot
Definition: infotrac.F90:6
leapfrog_mod::w
real, dimension(:,:), pointer, save w
Definition: leapfrog_mod.F90:15
pression
subroutine pression(ngrid, ap, bp, ps, p)
Definition: pression.F90:2
Tp_fcg_gcssold
!$Id && Tp_fcg_gcssold
Definition: 1D_interp_cases.h:9
day_ini
!$Id && day_ini
Definition: temps.h:15
pressure
!$Id vert_prof_dissip LOGICAL lstardis INTEGER niterh integer vert_prof_dissip!vertical profile of horizontal dissipation!Allowed function of pressure
Definition: comdissnew.h:13
nivsigs
!$Id nivsigs(llm)
false
!$Id itapm1 ENDIF!IM on interpole les champs sur les niveaux STD de pression!IM a chaque pas de temps de la physique c!positionnement de l argument logique a false c!pour ne pas recalculer deux fois la meme chose!c!a cet effet un appel a plevel_new a ete deplace c!a la fin de la serie d appels c!la boucle DO nlevSTD a ete internalisee c!dans d ou la creation de cette routine c c!CALL false
Definition: calcul_STDlev.h:26
etot0
!$Id etot0
Definition: ener.h:11
ioipsl_getincom
Definition: ioipsl_getincom.F90:6
rad
!$Id mode_top_bound COMMON comconstr rad
Definition: comconst.h:7
tg
!$Id tg
Definition: flux_arp.h:11
mod_1d_cases_read::plev_prof_cas
real, dimension(:), allocatable plev_prof_cas
Definition: mod_1D_cases_read.F90:55
dpres
!$Id dpres(llm)
dynredem_mod::modname
character(len=256), save, public modname
Definition: dynredem_mod.F90:9
conf_gcm.F90
ioipsl_getincom::getin
Definition: ioipsl_getincom.F90:22
day_ref
!$Id day_ref
Definition: temps.h:15
RPI
!$Id RNAVO!A1 Astronomical constants REAL ROMEGA!A1 bis Constantes concernant l orbite de la R_incl!A1 Geoide REAL R1SA!A1 Radiation!REAL RI0 REAL RSIGMA!A1 Thermodynamic gas phase REAL RCVV REAL RETV Thermodynamic solid phases REAL RCS!A1 Thermodynamic transition of phase REAL RATM!A1 Curve of saturation REAL RGAMS REAL RGAMD!COMMON YOMCST RPI
Definition: YOMCST.h:11
ts_fcg
!$Id ts_fcg
Definition: 1Dconv.h:4
fxyhypb
!$Id fxyhypb
Definition: logic.h:10
lmdz1d
program lmdz1d
Definition: lmdz1d.F90:8
physiq
subroutine physiq(nlon, nlev, debut, lafin, jD_cur, jH_cur, pdtphys, paprs, pplay, pphi, pphis, presnivs, u, v, rot, t, qx, flxmass_w, d_u, d_v, d_t, d_qx, d_ps, dudyn)
Definition: physiq.F90:11
call_calfis_mod::phi
real, dimension(:,:), pointer, save phi
Definition: call_calfis_mod.F90:17
cpp
!$Id mode_top_bound COMMON comconstr cpp
Definition: comconst.h:7
nid
nid INTEGER nid
Definition: iotd.h:16
call_calfis_mod::vcov
real, dimension(:,:), pointer, save vcov
Definition: call_calfis_mod.F90:5
day_ini_toga
!Declarations specifiques au cas Toga character *::fich_toga!integer nlev_prof nt_toga day_ini_toga
Definition: 1D_decl_cases.h:9
j
!$Id Turb_fcg!implicit none!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc!cette routine permet d obtenir hq et ainsi de!pouvoir calculer la convergence et le cisaillement dans la physiq!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc INTEGER klev REAL j
Definition: 1Dconv.h:27
ok_flux_surf
!$Id ok_flux_surf
Definition: flux_arp.h:11
gradiv2_mod::div
real, dimension(:,:), pointer, save div
Definition: gradiv2_mod.F90:5
mod_phys_lmdz_para
Definition: mod_phys_lmdz_para.F90:4
pplay
!$Id itapm1 ENDIF!IM on interpole les champs sur les niveaux STD de pression!IM a chaque pas de temps de la physique c!positionnement de l argument logique a false c!pour ne pas recalculer deux fois la meme chose!c!a cet effet un appel a plevel_new a ete deplace c!a la fin de la serie d appels c!la boucle DO nlevSTD a ete internalisee c!dans d ou la creation de cette routine c c!CALL pplay
Definition: calcul_STDlev.h:26
iostart::close_restartphy
subroutine, public close_restartphy
Definition: iostart.F90:341
height
!$Id nt_fire height
Definition: 1D_read_forc_cases.h:21
day
!$Id Turb_fcg_gcssold get_uvd hqturb_gcssold endif!large scale llm &&& day
Definition: 1D_interp_cases.h:45
imp_fcg
!$Id imp_fcg
Definition: 1Dconv.h:4
prt_level
!FH On elimine toutes les clefs physiques dans la dynamique prt_level
Definition: write_paramLMDZ_dyn.h:9
iostart::close_startphy
subroutine, public close_startphy
Definition: iostart.F90:50
daysec
!$Id mode_top_bound COMMON comconstr daysec
Definition: comconst.h:7
paprs
!$Id itapm1 ENDIF!IM on interpole les champs sur les niveaux STD de pression!IM a chaque pas de temps de la physique c!positionnement de l argument logique a false c!pour ne pas recalculer deux fois la meme chose!c!a cet effet un appel a plevel_new a ete deplace c!a la fin de la serie d appels c!la boucle DO nlevSTD a ete internalisee c!dans d ou la creation de cette routine c c!CALL ulevSTD CALL &zphi philevSTD CALL &zx_rh rhlevSTD!DO klev DO klon klev DO klon klev DO klon klev DO klon klev DO klon klev DO klon klev DO klon klev DO klon klev DO klon klev DO klon du jour ou toutes les read_climoz CALL &zphi geo500!IM on interpole a chaque pas de temps le paprs
Definition: calcul_STDlev.h:144
IF
!$Id we still need to ch2 ch3 ch4 *Initilisation parametres du lmdz1D lunout IF(lunout/=5.and.lunout/=6) THEN ENDIF!Config Key
ztot0
!$Id ztot0
Definition: ener.h:11
phis
!$Id ***************************************!ECRITURE DU phis
Definition: write_histrac.h:9
u
!$Id Turb_fcg_gcssold get_uvd hqturb_gcssold endif!large scale llm day day1 day day1 *dt_toga endif!time annee_ref dt_toga u_toga vq_toga w_prof vq_prof llm u(l)
nudge_w
!$Id tend_rayo nudge_w
Definition: compar1d.h:5
x
!$Header!c c INCLUDE fxyprim h c c c Fonctions in line c c REAL fyprim REAL rj c c il faut la calculer avant d appeler ces fonctions c c c Fonctions a changer selon x(x) et y(y) choisis.c-----------------------------------------------------------------c c.....ici
disvert
subroutine disvert()
Definition: disvert.F90:4
pa
!$Id && pa
Definition: comvert.h:8
write_field_phy
Definition: write_field_phy.F90:4
timer_filtre::time
real time
Definition: timer_filtre.F90:4
true
!$Id itapm1 ENDIF!IM on interpole les champs sur les niveaux STD de pression!IM a chaque pas de temps de la physique c!positionnement de l argument logique a false c!pour ne pas recalculer deux fois la meme chose!c!a cet effet un appel a plevel_new a ete deplace c!a la fin de la serie d appels c!la boucle DO nlevSTD a ete internalisee c!dans d ou la creation de cette routine c c!CALL ulevSTD CALL &zphi philevSTD CALL &zx_rh rhlevSTD!DO klev DO klon klev DO klon klev DO klon klev DO klon klev DO klon klev DO klon klev DO klon klev DO klon klev DO klon klev DO klon du jour ou toutes les read_climoz CALL true
Definition: calcul_STDlev.h:139
limit_read_mod::rugos
real, dimension(:), allocatable, save, private rugos
Definition: limit_read_mod.F90:18
rlat
c c $Id c nbregdyn DO klon c rlat(i) c ENDIF!lon c ENDIF!lat ENDIF!pctsrf ENDDO!klon ENDDO!nbregdyn cIM 190504 ENDIF!ok_regdyn cIM somme de toutes les nhistoW BEG IF(debut) THEN DO nreg
vertical_layers_mod::ap
real, dimension(:), allocatable, save ap
Definition: vertical_layers_mod.F90:7
dzoomx
!$Header!c!c!c include serre h!c REAL dzoomx
Definition: serre.h:8
lllm
!$Id lllm
Definition: comconst.h:7
sur
c c zjulian c cym CALL iim cym klev iim cym jjmp1 cym On stoke le fichier bilKP instantanne sur
Definition: ini_bilKP_ins.h:41
alpha
!$Id!Parameters for nlm real spfac!IM cf epmax real ptcrit real omtrain real dttrig real alpha real delta real betad COMMON cv30param nlm spfac &!IM cf ptcrit omtrain dttrig alpha
Definition: cv30param.h:5
phys_state_var_mod::ustar
real, dimension(:,:), allocatable, save ustar
Definition: phys_state_var_mod.F90:370
print_control_mod
Definition: print_control_mod.F90:2
albedo
Definition: albedo.F90:2
grossismy
!$Header!c!c!c include serre h!c REAL grossismy
Definition: serre.h:8
dtvr
!$Id mode_top_bound COMMON comconstr dtvr
Definition: comconst.h:7
call_calfis_mod::du
real, dimension(:,:), pointer, save du
Definition: call_calfis_mod.F90:18
tpm_fft::t
type(fft_type), pointer t
Definition: tpm_fft.F90:14
taux
!$Header!c!c!c include serre h!c REAL taux
Definition: serre.h:8
tend_t
!$Id tend_t
Definition: compar1d.h:5
ONLY
!$Id we still need to ONLY
Definition: 1DUTILS.h:34
mod_phy____grd::sigma
real(kind=real8), dimension(:), allocatable, save sigma
Definition: Mod_PHY____grd.f90:59
iim
c c zjulian c cym CALL iim cym klev iim
Definition: ini_bilKP_ave.h:24
tend_w
!$Id tend_w
Definition: compar1d.h:5
phys_state_var_mod::agesno
real, dimension(:,:), allocatable, save agesno
Definition: phys_state_var_mod.F90:28
phys_state_var_mod::swup
real, dimension(:,:), allocatable, save swup
Definition: phys_state_var_mod.F90:103
tend_v
!$Id tend_v
Definition: compar1d.h:5
mod_grid_phy_lmdz
Definition: mod_grid_phy_lmdz.F90:4
radiation_ar4_param::at
real(kind=8), dimension(8, 3), parameter at
Definition: radiation_ar4_param.F90:223
stot0
!$Id stot0
Definition: ener.h:11
clat
!$Header!c!c!c include serre h!c REAL clat
Definition: serre.h:8
use
!$Id we still need to use(a local version of) getin use ioipsl_getincom USE print_control_mod
u_mod
do llm!au dessus on relaxe vers profil init!on fait l hypothese que dans ce il n y a plus d eau liq au dessus!donc la relaxation en thetal et qt devient relaxation en tempe et qv u_mod(l)!if(l.ge.llm700) then relax_q(l
itaufin
!$Id itaufin
Definition: temps.h:15
play
!$Id Turb_fcg!implicit none!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc!cette routine permet d obtenir hq et ainsi de!pouvoir calculer la convergence et le cisaillement dans la physiq!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc INTEGER klev REAL in CHARACTER file_fordat COMMON com1_phys_gcss play
Definition: 1Dconv.h:27
RTT
!$Id RNAVO!A1 Astronomical constants REAL ROMEGA!A1 bis Constantes concernant l orbite de la R_incl!A1 Geoide REAL R1SA!A1 Radiation!REAL RI0 REAL RSIGMA!A1 Thermodynamic gas phase REAL RCVV REAL RETV Thermodynamic solid phases REAL RCS!A1 Thermodynamic transition of phase REAL RTT
Definition: YOMCST.h:11
def
!$Id pressure_exner real ap!hybrid pressure contribution at interlayers real bp!hybrid sigma contribution at interlayer real based on!preff and scaleheight integer disvert_type!type of vertical!automatic!using z2sig def(or 'esasig.def) file logical pressure_exner!compute pressure inside layers using Exner function
dynredem_mod::put_var
subroutine, public put_var(ncid, var, title, did, v, units)
Definition: dynredem_mod.F90:198
control_mod::nday
integer, save nday
Definition: control_mod.F90:14
radiation_ar4_param::d
real(kind=8), dimension(2, 3), parameter d
Definition: radiation_ar4_param.F90:211
de
do llm!au dessus de
Definition: 1D_nudge_sandu_astex.h:1
klevp1
!$Id itapm1 ENDIF!IM on interpole les champs sur les niveaux STD de pression!IM a chaque pas de temps de la physique c!positionnement de l argument logique a false c!pour ne pas recalculer deux fois la meme chose!c!a cet effet un appel a plevel_new a ete deplace c!a la fin de la serie d appels c!la boucle DO nlevSTD a ete internalisee c!dans d ou la creation de cette routine c c!CALL ulevSTD CALL &zphi philevSTD CALL &zx_rh rhlevSTD!DO klev DO klon klev DO klon klev DO klon klev DO klon klev DO klon klev DO klon klev DO klon klev DO klon klev DO klon klev DO klon du jour ou toutes les read_climoz CALL &zphi geo500!IM on interpole a chaque pas de temps le klevp1
Definition: calcul_STDlev.h:144
phys_state_var_mod
Definition: phys_state_var_mod.F90:4
call_calfis_mod::ucov
real, dimension(:,:), pointer, save ucov
Definition: call_calfis_mod.F90:4
PARAMETER
!$Id!Valid and equivalent for either free source form or fixed source form INTEGER i1_fin INTEGER i2_fin!cc PARAMETER(i1_deb=21, i1_fin=40)!cc PARAMETER(i2_deb
flux_arp.h
vertical_layers_mod::bp
real, dimension(:), allocatable, save bp
Definition: vertical_layers_mod.F90:9
phys_local_var_mod::omega
real, dimension(:,:), allocatable, save omega
Definition: phys_local_var_mod.F90:342
Tp_ini_gcssold
!$Id Tp_fcg_gcssold logical::Tp_ini_gcssold logical::xTurb_fcg_gcssold common fcg_gcssold && Tp_ini_gcssold
Definition: fcg_gcssold.h:4
ptot0
!$Id ptot0
Definition: ener.h:11
dimphy
Definition: dimphy.F90:1
phys_state_var_mod::lwup
real, dimension(:,:), allocatable, save lwup
Definition: phys_state_var_mod.F90:111
fcg_gcssold.h
real
!IM for NMC files!real nfiles!real nfiles!real nfiles!real nlevSTD3 real nlevSTD3 real nlevSTD3 real nlevSTD3 real nlevSTD8 real nlevSTD8 real nlevSTD8 real nlevSTD8 real
Definition: declare_STDlev.h:27
phys_local_var_mod::theta
real, dimension(:,:), allocatable, save theta
Definition: phys_local_var_mod.F90:349
press_coefoz_m::plev
real, dimension(:), pointer, save plev
Definition: press_coefoz_m.F90:6
fcg_racmo.h
nudge_v
!$Id tend_rayo nudge_v
Definition: compar1d.h:5
control_mod
Definition: control_mod.F90:5
nivsig
!$Id nivsig(llm+1)
gr_fi_dyn
subroutine gr_fi_dyn(nfield, ngrid, im, jm, pfi, pdyn)
Definition: gr_fi_dyn.F:5
iophy
Definition: iophy.F90:4
rlon
c c $Id c nbregdyn DO klon c rlon(i)
cas
do llm!au dessus on relaxe vers profil init!on fait l hypothese que dans ce cas
Definition: 1D_nudge_sandu_astex.h:1
iostart::open_restartphy
subroutine, public open_restartphy(filename)
Definition: iostart.F90:312
vertical_layers_mod::presnivs
real, dimension(:), allocatable, save presnivs
Definition: vertical_layers_mod.F90:11
annee_ref
!$Id annee_ref
Definition: temps.h:15
lunout
!$Header!gestion des impressions de sorties et de débogage la sortie standard prt_level COMMON comprint lunout
Definition: iniprint.h:7
PB
INTERFACE SUBROUTINE LWB && PB
Definition: lwb.intfb.h:3
geometry_mod::latitude
real, dimension(:), allocatable, save latitude
Definition: geometry_mod.F90:8
soil
subroutine soil(ptimestep, indice, knon, snow, ptsrf, ptsoil, pcapcal, pfluxgrd)
Definition: soil.F90:6