source: LMDZ5/branches/LMDZ5-DOFOCO/libf/phy1d/1DUTILS.h_with_writelim @ 5228

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

In 1DUTILS.h disvert renamed disvert0.
There are now 2 routines to get a vertical grid: disvert0 or disvert (../dyn3d/disvert.F90).
The new flag ok_old_disvert allows the use of either disvert0 (ok_old_disvert=y) or disvert (ok_old_disvert=n)

In lmdz1d.F: rlat_rad and rlon_rad for coordinates in radians, cfdt ans sfdt concern the geostrophic wind

File size: 80.7 KB
Line 
1!
2! $Id: conf_unicol.F 1279 2010-08-04 17:20:56Z lahellec $
3!
4c
5c
6      SUBROUTINE conf_unicol( tapedef )
7c
8#ifdef CPP_IOIPSL
9      use IOIPSL
10#else
11! if not using IOIPSL, we still need to use (a local version of) getin
12      use ioipsl_getincom
13#endif
14      IMPLICIT NONE
15c-----------------------------------------------------------------------
16c     Auteurs :   A. Lahellec  .
17c
18c     Arguments :
19c
20c     tapedef   :
21
22       INTEGER tapedef
23c
24c   Declarations :
25c   --------------
26
27#include "compar1d.h"
28#include "flux_arp.h"
29#include "fcg_gcssold.h"
30#include "fcg_racmo.h"
31#include "iniprint.h"
32c
33c
34c   local:
35c   ------
36
37c      CHARACTER ch1*72,ch2*72,ch3*72,ch4*12
38     
39c
40c  -------------------------------------------------------------------
41c
42c      .........    Initilisation parametres du lmdz1D      ..........
43c
44c---------------------------------------------------------------------
45c   initialisations:
46c   ----------------
47
48!Config  Key  = lunout
49!Config  Desc = unite de fichier pour les impressions
50!Config  Def  = 6
51!Config  Help = unite de fichier pour les impressions
52!Config         (defaut sortie standard = 6)
53      lunout=6
54!      CALL getin('lunout', lunout)
55      IF (lunout /= 5 .and. lunout /= 6) THEN
56        OPEN(lunout,FILE='lmdz.out')
57      ENDIF
58
59!Config  Key  = prt_level
60!Config  Desc = niveau d'impressions de débogage
61!Config  Def  = 0
62!Config  Help = Niveau d'impression pour le débogage
63!Config         (0 = minimum d'impression)
64c      prt_level = 0
65c      CALL getin('prt_level',prt_level)
66
67c-----------------------------------------------------------------------
68c  Parametres de controle du run:
69c-----------------------------------------------------------------------
70
71!Config  Key  = restart
72!Config  Desc = on repart des startphy et start1dyn
73!Config  Def  = false
74!Config  Help = les fichiers restart doivent etre renomme en start
75       restart = .FALSE.
76       CALL getin('restart',restart)
77
78!Config  Key  = forcing_type
79!Config  Desc = defines the way the SCM is forced:
80!Config  Def  = 0
81!!Config  Help = 0 ==> forcing_les = .true.
82!             initial profiles from file prof.inp.001
83!             no forcing by LS convergence ;
84!             surface temperature imposed ;
85!             radiative cooling may be imposed (iflag_radia=0 in physiq.def)
86!         = 1 ==> forcing_radconv = .true.
87!             idem forcing_type = 0, but the imposed radiative cooling
88!             is set to 0 (hence, if iflag_radia=0 in physiq.def,
89!             then there is no radiative cooling at all)
90!         = 2 ==> forcing_toga = .true.
91!             initial profiles from TOGA-COARE IFA files
92!             LS convergence and SST imposed from TOGA-COARE IFA files
93!         = 3 ==> forcing_GCM2SCM = .true.
94!             initial profiles from the GCM output
95!             LS convergence imposed from the GCM output
96!         = 4 ==> forcing_twpi = .true.
97!             initial profiles from TWPICE nc files
98!             LS convergence and SST imposed from TWPICE nc files
99!         = 5 ==> forcing_rico = .true.
100!             initial profiles from RICO idealized
101!             LS convergence imposed from  RICO (cst)
102!         = 40 ==> forcing_GCSSold = .true.
103!             initial profile from GCSS file
104!             LS convergence imposed from GCSS file
105!
106       forcing_type = 0
107       CALL getin('forcing_type',forcing_type)
108         imp_fcg_gcssold   = .false.
109         ts_fcg_gcssold    = .false. 
110         Tp_fcg_gcssold    = .false. 
111         Tp_ini_gcssold    = .false. 
112         xTurb_fcg_gcssold = .false.
113        IF (forcing_type .eq.40) THEN
114          CALL getin('imp_fcg',imp_fcg_gcssold)
115          CALL getin('ts_fcg',ts_fcg_gcssold)
116          CALL getin('tp_fcg',Tp_fcg_gcssold)
117          CALL getin('tp_ini',Tp_ini_gcssold)
118          CALL getin('turb_fcg',xTurb_fcg_gcssold)
119        ENDIF
120
121!Config  Key  = ok_flux_surf
122!Config  Desc = forcage ou non par les flux de surface
123!Config  Def  = false
124!Config  Help = forcage ou non par les flux de surface
125       ok_flux_surf = .FALSE.
126       CALL getin('ok_flux_surf',ok_flux_surf)
127
128!Config  Key  = ok_old_disvert
129!Config  Desc = utilisation de l ancien programme disvert0 (dans 1DUTILS.h)
130!Config  Def  = false
131!Config  Help = utilisation de l ancien programme disvert0 (dans 1DUTILS.h)
132       ok_old_disvert = .FALSE.
133       CALL getin('ok_old_disvert',ok_old_disvert)
134
135!Config  Key  = time_ini
136!Config  Desc = meaningless in this  case
137!Config  Def  = 0.
138!Config  Help =
139       tsurf = 0.
140       CALL getin('time_ini',time_ini)
141
142!Config  Key  = rlat et rlon
143!Config  Desc = latitude et longitude
144!Config  Def  = 0.0  0.0
145!Config  Help = fixe la position de la colonne
146       xlat = 0.
147       xlon = 0.
148       CALL getin('rlat',xlat)
149       CALL getin('rlon',xlon)
150
151!Config  Key  = airephy
152!Config  Desc = Grid cell area
153!Config  Def  = 1.e11
154!Config  Help =
155       airefi = 1.e11
156       CALL getin('airephy',airefi)
157
158!Config  Key  = nat_surf
159!Config  Desc = surface type
160!Config  Def  = 0 (ocean)
161!Config  Help = 0=ocean,1=land,2=glacier,3=banquise
162       nat_surf = 0.
163       CALL getin('nat_surf',nat_surf)
164
165!Config  Key  = tsurf
166!Config  Desc = surface temperature
167!Config  Def  = 290.
168!Config  Help = not used if type_ts_forcing=1 in lmdz1d.F
169       tsurf = 290.
170       CALL getin('tsurf',tsurf)
171
172!Config  Key  = psurf
173!Config  Desc = surface pressure
174!Config  Def  = 102400.
175!Config  Help =
176       psurf = 102400.
177       CALL getin('psurf',psurf)
178
179!Config  Key  = zsurf
180!Config  Desc = surface altitude
181!Config  Def  = 0.
182!Config  Help =
183       zsurf = 0.
184       CALL getin('zsurf',zsurf)
185
186!Config  Key  = rugos
187!Config  Desc = coefficient de frottement
188!Config  Def  = 0.0001
189!Config  Help = calcul du Cdrag
190       rugos = 0.0001
191       CALL getin('rugos',rugos)
192
193!Config  Key  = wtsurf et wqsurf
194!Config  Desc = ???
195!Config  Def  = 0.0 0.0
196!Config  Help =
197       wtsurf = 0.0
198       wqsurf = 0.0
199       CALL getin('wtsurf',wtsurf)
200       CALL getin('wqsurf',wqsurf)
201
202!Config  Key  = albedo
203!Config  Desc = albedo
204!Config  Def  = 0.09
205!Config  Help =
206       albedo = 0.09
207       CALL getin('albedo',albedo)
208
209!Config  Key  = agesno
210!Config  Desc = age de la neige
211!Config  Def  = 30.0
212!Config  Help =
213       xagesno = 30.0
214       CALL getin('agesno',xagesno)
215
216!Config  Key  = restart_runoff
217!Config  Desc = age de la neige
218!Config  Def  = 30.0
219!Config  Help =
220       restart_runoff = 0.0
221       CALL getin('restart_runoff',restart_runoff)
222
223!Config  Key  = qsolinp
224!Config  Desc = initial bucket water content (kg/m2) when land (5std)
225!Config  Def  = 30.0
226!Config  Help =
227       qsolinp = 1.
228       CALL getin('qsolinp',qsolinp)
229
230!Config  Key  = zpicinp
231!Config  Desc = denivellation orographie
232!Config  Def  = 300.
233!Config  Help =  input brise
234       zpicinp = 300.
235       CALL getin('zpicinp',zpicinp)
236
237
238      write(lunout,*)' +++++++++++++++++++++++++++++++++++++++'
239      write(lunout,*)' Configuration des parametres du gcm1D: '
240      write(lunout,*)' +++++++++++++++++++++++++++++++++++++++'
241      write(lunout,*)' restart = ', restart
242      write(lunout,*)' forcing_type = ', forcing_type
243      write(lunout,*)' time_ini = ', time_ini
244      write(lunout,*)' rlat = ', xlat
245      write(lunout,*)' rlon = ', xlon
246      write(lunout,*)' airephy = ', airefi
247      write(lunout,*)' nat_surf = ', nat_surf
248      write(lunout,*)' tsurf = ', tsurf
249      write(lunout,*)' psurf = ', psurf
250      write(lunout,*)' zsurf = ', zsurf
251      write(lunout,*)' rugos = ', rugos
252      write(lunout,*)' wtsurf = ', wtsurf
253      write(lunout,*)' wqsurf = ', wqsurf
254      write(lunout,*)' albedo = ', albedo
255      write(lunout,*)' xagesno = ', xagesno
256      write(lunout,*)' restart_runoff = ', restart_runoff
257      write(lunout,*)' qsolinp = ', qsolinp
258      write(lunout,*)' zpicinp = ', zpicinp
259      IF (forcing_type .eq.40) THEN
260        write(lunout,*) '--- Forcing type GCSS Old --- with:'
261        write(lunout,*)'imp_fcg',imp_fcg_gcssold
262        write(lunout,*)'ts_fcg',ts_fcg_gcssold
263        write(lunout,*)'tp_fcg',Tp_fcg_gcssold
264        write(lunout,*)'tp_ini',Tp_ini_gcssold
265        write(lunout,*)'xturb_fcg',xTurb_fcg_gcssold
266      ENDIF
267
268      write(lunout,*)' +++++++++++++++++++++++++++++++++++++++'
269      write(lunout,*)
270c
271      RETURN
272      END
273!
274! $Id: dyn1deta0.F 1279 2010/07/30 A Lahellec$
275!
276c
277      SUBROUTINE dyn1deta0(fichnom,plev,play,phi,phis,presnivs,
278     &                          ucov,vcov,temp,q,omega2)
279      USE dimphy
280      USE mod_grid_phy_lmdz
281      USE mod_phys_lmdz_para
282      USE iophy
283      USE phys_state_var_mod
284      USE iostart
285      USE write_field_phy
286      USE infotrac
287      use control_mod
288
289      IMPLICIT NONE
290c=======================================================
291c Ecriture du fichier de redemarrage sous format NetCDF
292c=======================================================
293c   Declarations:
294c   -------------
295#include "dimensions.h"
296#include "comconst.h"
297#include "temps.h"
298!!#include "control.h"
299#include "logic.h"
300#include "netcdf.inc"
301
302c   Arguments:
303c   ----------
304      CHARACTER*(*) fichnom
305cAl1 plev tronque pour .nc mais plev(klev+1):=0
306      real :: plev(klon,klev),play (klon,klev),phi(klon,klev)
307      real :: presnivs(klon,klev)
308      real :: ucov(klon,klev),vcov(klon,klev),temp(klon,klev)
309      real :: q(klon,klev,nqtot),omega2(klon,klev)
310c      real :: ug(klev),vg(klev),fcoriolis
311      real :: phis(klon)
312
313c   Variables locales pour NetCDF:
314c   ------------------------------
315      INTEGER nid, nvarid
316      INTEGER idim_s
317      INTEGER ierr, ierr_file
318      INTEGER iq
319      INTEGER length
320      PARAMETER (length = 100)
321      REAL tab_cntrl(length) ! tableau des parametres du run
322      character*4 nmq(nqtot)
323      character*12 modname
324      character*80 abort_message
325      LOGICAL found
326c
327      INTEGER nb
328
329      modname = 'dyn1deta0 : '
330      nmq(1)="vap"
331      nmq(2)="cond"
332      do iq=3,nqtot
333        write(nmq(iq),'("tra",i1)') iq-2
334      enddo
335      print*,'in dyn1deta0 ',fichnom,klon,klev,nqtot
336      CALL open_startphy(fichnom)
337      print*,'after open startphy ',fichnom,nmq
338
339c
340c Lecture des parametres de controle:
341c
342      CALL get_var("controle",tab_cntrl)
343       
344
345      im         = tab_cntrl(1)
346      jm         = tab_cntrl(2)
347      lllm       = tab_cntrl(3)
348      day_ref    = tab_cntrl(4)
349      annee_ref  = tab_cntrl(5)
350c      rad        = tab_cntrl(6)
351c      omeg       = tab_cntrl(7)
352c      g          = tab_cntrl(8)
353c      cpp        = tab_cntrl(9)
354c      kappa      = tab_cntrl(10)
355c      daysec     = tab_cntrl(11)
356c      dtvr       = tab_cntrl(12)
357c      etot0      = tab_cntrl(13)
358c      ptot0      = tab_cntrl(14)
359c      ztot0      = tab_cntrl(15)
360c      stot0      = tab_cntrl(16)
361c      ang0       = tab_cntrl(17)
362c      pa         = tab_cntrl(18)
363c      preff      = tab_cntrl(19)
364c
365c      clon       = tab_cntrl(20)
366c      clat       = tab_cntrl(21)
367c      grossismx  = tab_cntrl(22)
368c      grossismy  = tab_cntrl(23)
369c
370      IF ( tab_cntrl(24).EQ.1. )  THEN
371        fxyhypb  = . TRUE .
372c        dzoomx   = tab_cntrl(25)
373c        dzoomy   = tab_cntrl(26)
374c        taux     = tab_cntrl(28)
375c        tauy     = tab_cntrl(29)
376      ELSE
377        fxyhypb = . FALSE .
378        ysinus  = . FALSE .
379        IF( tab_cntrl(27).EQ.1. ) ysinus = . TRUE.
380      ENDIF
381
382      day_ini = tab_cntrl(30)
383      itau_dyn = tab_cntrl(31)
384c   .................................................................
385c
386c
387c      PRINT*,'rad,omeg,g,cpp,kappa',rad,omeg,g,cpp,kappa
388cAl1
389       Print*,'day_ref,annee_ref,day_ini,itau_dyn',
390     &              day_ref,annee_ref,day_ini,itau_dyn
391
392c  Lecture des champs
393c
394      plev(1,klev+1)=0.
395      CALL get_field("plev",plev,found)
396      IF (.NOT. found) PRINT*, modname//'Le champ <Plev> est absent'
397      CALL get_field("play",play,found)
398      IF (.NOT. found) PRINT*, modname//'Le champ <Play> est absent'
399      CALL get_field("phi",phi,found)
400      IF (.NOT. found) PRINT*, modname//'Le champ <Phi> est absent'
401      CALL get_field("phis",phis,found)
402      IF (.NOT. found) PRINT*, modname//'Le champ <Phis> est absent'
403      CALL get_field("presnivs",presnivs,found)
404      IF (.NOT. found) PRINT*, modname//'Le champ <Presnivs> est absent'
405      CALL get_field("ucov",ucov,found)
406      IF (.NOT. found) PRINT*, modname//'Le champ <ucov> est absent'
407      CALL get_field("vcov",vcov,found)
408      IF (.NOT. found) PRINT*, modname//'Le champ <vcov> est absent'
409      CALL get_field("temp",temp,found)
410      IF (.NOT. found) PRINT*, modname//'Le champ <temp> est absent'
411      CALL get_field("omega2",omega2,found)
412      IF (.NOT. found) PRINT*, modname//'Le champ <omega2> est absent'
413
414      Do iq=1,nqtot
415        CALL get_field("q"//nmq(iq),q(:,:,iq),found)
416        IF (.NOT. found)
417     &  PRINT*, modname//'Le champ <q'//nmq//'> est absent'
418      EndDo
419
420      CALL close_startphy
421      print*,' close startphy'
422     .      ,fichnom,play(1,1),play(1,klev),temp(1,klev)
423c
424      RETURN
425      END
426!
427! $Id: dyn1dredem.F 1279 2010/07/29 A Lahellec$
428!
429c
430      SUBROUTINE dyn1dredem(fichnom,plev,play,phi,phis,presnivs,
431     &                          ucov,vcov,temp,q,omega2)
432      USE dimphy
433      USE mod_grid_phy_lmdz
434      USE mod_phys_lmdz_para
435      USE phys_state_var_mod
436      USE iostart
437      USE infotrac
438      use control_mod
439
440      IMPLICIT NONE
441c=======================================================
442c Ecriture du fichier de redemarrage sous format NetCDF
443c=======================================================
444c   Declarations:
445c   -------------
446#include "dimensions.h"
447#include "comconst.h"
448#include "temps.h"
449!!#include "control.h"
450#include "logic.h"
451#include "netcdf.inc"
452
453c   Arguments:
454c   ----------
455      CHARACTER*(*) fichnom
456      REAL time
457cAl1 plev tronque pour .nc mais plev(klev+1):=0
458      real :: plev(klon,klev),play (klon,klev),phi(klon,klev)
459      real :: presnivs(klon,klev)
460      real :: ucov(klon,klev),vcov(klon,klev),temp(klon,klev)
461      real :: q(klon,klev,nqtot)
462      real :: omega2(klon,klev),rho(klon,klev+1)
463c      real :: ug(klev),vg(klev),fcoriolis
464      real :: phis(klon)
465
466c   Variables locales pour NetCDF:
467c   ------------------------------
468      INTEGER nid, nvarid
469      INTEGER idim_s
470      INTEGER ierr, ierr_file
471      INTEGER iq,l
472      INTEGER length
473      PARAMETER (length = 100)
474      REAL tab_cntrl(length) ! tableau des parametres du run
475      character*4 nmq(nqtot)
476      character*20 modname
477      character*80 abort_message
478c
479      INTEGER nb
480      SAVE nb
481      DATA nb / 0 /
482
483      REAL zan0,zjulian,hours
484      INTEGER yyears0,jjour0, mmois0
485      character*30 unites
486
487cDbg
488      CALL open_restartphy(fichnom)
489      print*,'redm1 ',fichnom,klon,klev,nqtot
490      nmq(1)="vap"
491      nmq(2)="cond"
492      nmq(3)="tra1"
493      nmq(4)="tra2"
494
495      modname = 'dyn1dredem'
496      ierr = NF_OPEN(fichnom, NF_WRITE, nid)
497      IF (ierr .NE. NF_NOERR) THEN
498         PRINT*, "Pb. d ouverture "//fichnom
499         CALL abort
500      ENDIF
501
502      DO l=1,length
503       tab_cntrl(l) = 0.
504      ENDDO
505       tab_cntrl(1)  = FLOAT(iim)
506       tab_cntrl(2)  = FLOAT(jjm)
507       tab_cntrl(3)  = FLOAT(llm)
508       tab_cntrl(4)  = FLOAT(day_ref)
509       tab_cntrl(5)  = FLOAT(annee_ref)
510       tab_cntrl(6)  = rad
511       tab_cntrl(7)  = omeg
512       tab_cntrl(8)  = g
513       tab_cntrl(9)  = cpp
514       tab_cntrl(10) = kappa
515       tab_cntrl(11) = daysec
516       tab_cntrl(12) = dtvr
517c       tab_cntrl(13) = etot0
518c       tab_cntrl(14) = ptot0
519c       tab_cntrl(15) = ztot0
520c       tab_cntrl(16) = stot0
521c       tab_cntrl(17) = ang0
522c       tab_cntrl(18) = pa
523c       tab_cntrl(19) = preff
524c
525c    .....    parametres  pour le zoom      ......   
526
527c       tab_cntrl(20)  = clon
528c       tab_cntrl(21)  = clat
529c       tab_cntrl(22)  = grossismx
530c       tab_cntrl(23)  = grossismy
531c
532      IF ( fxyhypb )   THEN
533       tab_cntrl(24) = 1.
534c       tab_cntrl(25) = dzoomx
535c       tab_cntrl(26) = dzoomy
536       tab_cntrl(27) = 0.
537c       tab_cntrl(28) = taux
538c       tab_cntrl(29) = tauy
539      ELSE
540       tab_cntrl(24) = 0.
541c       tab_cntrl(25) = dzoomx
542c       tab_cntrl(26) = dzoomy
543       tab_cntrl(27) = 0.
544       tab_cntrl(28) = 0.
545       tab_cntrl(29) = 0.
546       IF( ysinus )  tab_cntrl(27) = 1.
547      ENDIF
548CAl1 iday_end -> day_end
549       tab_cntrl(30) = FLOAT(day_end)
550       tab_cntrl(31) = FLOAT(itau_dyn + itaufin)
551c
552      CALL put_var("controle","Param. de controle Dyn1D",tab_cntrl)
553c
554
555c  Ecriture/extension de la coordonnee temps
556
557      nb = nb + 1
558
559c  Ecriture des champs
560c
561      CALL put_field("plev","p interfaces sauf la nulle",plev)
562      CALL put_field("play","",play)
563      CALL put_field("phi","geopotentielle",phi)
564      CALL put_field("phis","geopotentiell de surface",phis)
565      CALL put_field("presnivs","",presnivs)
566      CALL put_field("ucov","",ucov)
567      CALL put_field("vcov","",vcov)
568      CALL put_field("temp","",temp)
569      CALL put_field("omega2","",omega2)
570
571      Do iq=1,nqtot
572        CALL put_field("q"//nmq(iq),"eau vap ou condens et traceurs",
573     .                                                      q(:,:,iq))
574      EndDo
575      CALL close_restartphy
576
577c
578      RETURN
579      END
580      SUBROUTINE gr_fi_dyn(nfield,ngrid,im,jm,pfi,pdyn)
581      IMPLICIT NONE
582!=======================================================================
583!   passage d'un champ de la grille scalaire a la grille physique
584!=======================================================================
585 
586!-----------------------------------------------------------------------
587!   declarations:
588!   -------------
589 
590      INTEGER im,jm,ngrid,nfield
591      REAL pdyn(im,jm,nfield)
592      REAL pfi(ngrid,nfield)
593 
594      INTEGER i,j,ifield,ig
595 
596!-----------------------------------------------------------------------
597!   calcul:
598!   -------
599 
600      DO ifield=1,nfield
601!   traitement des poles
602         DO i=1,im
603            pdyn(i,1,ifield)=pfi(1,ifield)
604            pdyn(i,jm,ifield)=pfi(ngrid,ifield)
605         ENDDO
606 
607!   traitement des point normaux
608         DO j=2,jm-1
609            ig=2+(j-2)*(im-1)
610            CALL SCOPY(im-1,pfi(ig,ifield),1,pdyn(1,j,ifield),1)
611            pdyn(im,j,ifield)=pdyn(1,j,ifield)
612         ENDDO
613      ENDDO
614 
615      RETURN
616      END
617 
618 
619
620      SUBROUTINE abort_gcm(modname, message, ierr)
621 
622      USE IOIPSL
623!
624! Stops the simulation cleanly, closing files and printing various
625! comments
626!
627!  Input: modname = name of calling program
628!         message = stuff to print
629!         ierr    = severity of situation ( = 0 normal )
630 
631      character*20 modname
632      integer ierr
633      character*80 message
634 
635      write(*,*) 'in abort_gcm'
636      call histclo
637!     call histclo(2)
638!     call histclo(3)
639!     call histclo(4)
640!     call histclo(5)
641      write(*,*) 'out of histclo'
642      write(*,*) 'Stopping in ', modname
643      write(*,*) 'Reason = ',message
644      call getin_dump
645!
646      if (ierr .eq. 0) then
647        write(*,*) 'Everything is cool'
648      else
649        write(*,*) 'Houston, we have a problem ', ierr
650      endif
651      STOP
652      END
653      REAL FUNCTION fq_sat(kelvin, millibar)
654!
655      IMPLICIT none
656!======================================================================
657! Autheur(s): Z.X. Li (LMD/CNRS)
658! Objet: calculer la vapeur d'eau saturante (formule Centre Euro.)
659!======================================================================
660! Arguments:
661! kelvin---input-R: temperature en Kelvin
662! millibar--input-R: pression en mb
663!
664! fq_sat----output-R: vapeur d'eau saturante en kg/kg
665!======================================================================
666!
667      REAL kelvin, millibar
668!
669      REAL r2es
670      PARAMETER (r2es=611.14 *18.0153/28.9644)
671!
672      REAL r3les, r3ies, r3es
673      PARAMETER (R3LES=17.269)
674      PARAMETER (R3IES=21.875)
675!
676      REAL r4les, r4ies, r4es
677      PARAMETER (R4LES=35.86)
678      PARAMETER (R4IES=7.66)
679!
680      REAL rtt
681      PARAMETER (rtt=273.16)
682!
683      REAL retv
684      PARAMETER (retv=28.9644/18.0153 - 1.0)
685!
686      REAL zqsat
687      REAL temp, pres
688!     ------------------------------------------------------------------
689!
690!
691      temp = kelvin
692      pres = millibar * 100.0
693!      write(*,*)'kelvin,millibar=',kelvin,millibar
694!      write(*,*)'temp,pres=',temp,pres
695!
696      IF (temp .LE. rtt) THEN
697         r3es = r3ies
698         r4es = r4ies
699      ELSE
700         r3es = r3les
701         r4es = r4les
702      ENDIF
703!
704      zqsat=r2es/pres * EXP ( r3es*(temp-rtt) / (temp-r4es) )
705      zqsat=MIN(0.5,ZQSAT)
706      zqsat=zqsat/(1.-retv  *zqsat)
707!
708      fq_sat = zqsat
709!
710      RETURN
711      END
712      subroutine scopy(n,sx,incx,sy,incy)
713!
714      IMPLICIT NONE
715!
716      integer n,incx,incy,ix,iy,i
717      real sx((n-1)*incx+1),sy((n-1)*incy+1)
718!
719      iy=1
720      ix=1
721      do 10 i=1,n
722         sy(iy)=sx(ix)
723         ix=ix+incx
724         iy=iy+incy
72510    continue
726!
727      return
728      end
729      subroutine wrgradsfi(if,nl,field,name,titlevar)
730      implicit none
731 
732!   Declarations
733 
734#include "gradsdef.h"
735 
736!   arguments
737      integer if,nl
738      real field(imx*jmx*lmx)
739      character*10 name,file
740      character*10 titlevar
741 
742!   local
743 
744      integer im,jm,lm,i,j,l,iv,iii,iji,iif,ijf
745 
746      logical writectl
747 
748 
749      writectl=.false.
750 
751!     print*,if,iid(if),jid(if),ifd(if),jfd(if)
752      iii=iid(if)
753      iji=jid(if)
754      iif=ifd(if)
755      ijf=jfd(if)
756      im=iif-iii+1
757      jm=ijf-iji+1
758      lm=lmd(if)
759
760
761!     print*,'im,jm,lm,name,firsttime(if)'
762!     print*,im,jm,lm,name,firsttime(if)
763 
764      if(firsttime(if)) then
765         if(name.eq.var(1,if)) then
766            firsttime(if)=.false.
767            ivar(if)=1
768         print*,'fin de l initialiation de l ecriture du fichier'
769         print*,file
770           print*,'fichier no: ',if
771           print*,'unit ',unit(if)
772           print*,'nvar  ',nvar(if)
773           print*,'vars ',(var(iv,if),iv=1,nvar(if))
774         else
775            ivar(if)=ivar(if)+1
776            nvar(if)=ivar(if)
777            var(ivar(if),if)=name
778            tvar(ivar(if),if)=trim(titlevar)
779            nld(ivar(if),if)=nl
780            print*,'initialisation ecriture de ',var(ivar(if),if)
781            print*,'if ivar(if) nld ',if,ivar(if),nld(ivar(if),if)
782         endif
783         writectl=.true.
784         itime(if)=1
785      else
786         ivar(if)=mod(ivar(if),nvar(if))+1
787         if (ivar(if).eq.nvar(if)) then
788            writectl=.true.
789            itime(if)=itime(if)+1
790         endif
791 
792         if(var(ivar(if),if).ne.name) then
793           print*,'Il faut stoker la meme succession de champs a chaque'
794           print*,'pas de temps'
795           print*,'fichier no: ',if
796           print*,'unit ',unit(if)
797           print*,'nvar  ',nvar(if)
798           print*,'vars ',(var(iv,if),iv=1,nvar(if))
799 
800           stop
801         endif
802      endif
803 
804!     print*,'ivar(if),nvar(if),var(ivar(if),if),writectl'
805!     print*,ivar(if),nvar(if),var(ivar(if),if),writectl
806      do l=1,nl
807         irec(if)=irec(if)+1
808!        print*,'Ecrit rec=',irec(if),iii,iif,iji,ijf,
809!    s (l-1)*imd(if)*jmd(if)+(iji-1)*imd(if)+iii
810!    s ,(l-1)*imd(if)*jmd(if)+(ijf-1)*imd(if)+iif
811         write(unit(if)+1,rec=irec(if))
812     s   ((field((l-1)*imd(if)*jmd(if)+(j-1)*imd(if)+i)
813     s   ,i=iii,iif),j=iji,ijf)
814      enddo
815      if (writectl) then
816 
817      file=fichier(if)
818!   WARNING! on reecrase le fichier .ctl a chaque ecriture
819      open(unit(if),file=trim(file)//'.ctl',
820     &         form='formatted',status='unknown')
821      write(unit(if),'(a5,1x,a40)')
822     &       'DSET ','^'//trim(file)//'.dat'
823 
824      write(unit(if),'(a12)') 'UNDEF 1.0E30'
825      write(unit(if),'(a5,1x,a40)') 'TITLE ',title(if)
826      call formcoord(unit(if),im,xd(iii,if),1.,.false.,'XDEF')
827      call formcoord(unit(if),jm,yd(iji,if),1.,.true.,'YDEF')
828      call formcoord(unit(if),lm,zd(1,if),1.,.false.,'ZDEF')
829      write(unit(if),'(a4,i10,a30)')
830     &       'TDEF ',itime(if),' LINEAR 07AUG1998 30MN '
831      write(unit(if),'(a4,2x,i5)') 'VARS',nvar(if)
832      do iv=1,nvar(if)
833!        print*,'if,var(iv,if),nld(iv,if),nld(iv,if)-1/nld(iv,if)'
834!        print*,if,var(iv,if),nld(iv,if),nld(iv,if)-1/nld(iv,if)
835         write(unit(if),1000) var(iv,if),nld(iv,if)-1/nld(iv,if)
836     &     ,99,tvar(iv,if)
837      enddo
838      write(unit(if),'(a7)') 'ENDVARS'
839!
8401000  format(a5,3x,i4,i3,1x,a39)
841 
842      close(unit(if))
843 
844      endif ! writectl
845 
846      return
847 
848      END
849 
850      subroutine inigrads(if,im
851     s  ,x,fx,xmin,xmax,jm,y,ymin,ymax,fy,lm,z,fz
852     s  ,dt,file,titlel)
853 
854 
855      implicit none
856 
857      integer if,im,jm,lm,i,j,l
858      real x(im),y(jm),z(lm),fx,fy,fz,dt
859      real xmin,xmax,ymin,ymax
860      integer nf
861 
862      character file*10,titlel*40
863 
864#include "gradsdef.h"
865 
866      data unit/24,32,34,36,38,40,42,44,46,48/
867      data nf/0/
868 
869      if (if.le.nf) stop'verifier les appels a inigrads'
870 
871      print*,'Entree dans inigrads'
872 
873      nf=if
874      title(if)=titlel
875      ivar(if)=0
876 
877      fichier(if)=trim(file)
878 
879      firsttime(if)=.true.
880      dtime(if)=dt
881 
882      iid(if)=1
883      ifd(if)=im
884      imd(if)=im
885      do i=1,im
886         xd(i,if)=x(i)*fx
887         if(xd(i,if).lt.xmin) iid(if)=i+1
888         if(xd(i,if).le.xmax) ifd(if)=i
889      enddo
890      print*,'On stoke du point ',iid(if),'  a ',ifd(if),' en x'
891 
892      jid(if)=1
893      jfd(if)=jm
894      jmd(if)=jm
895      do j=1,jm
896         yd(j,if)=y(j)*fy
897         if(yd(j,if).gt.ymax) jid(if)=j+1
898         if(yd(j,if).ge.ymin) jfd(if)=j
899      enddo
900      print*,'On stoke du point ',jid(if),'  a ',jfd(if),' en y'
901
902      print*,'Open de dat'
903      print*,'file=',file
904      print*,'fichier(if)=',fichier(if)
905 
906      print*,4*(ifd(if)-iid(if))*(jfd(if)-jid(if))
907      print*,trim(file)//'.dat'
908 
909      OPEN (unit(if)+1,FILE=trim(file)//'.dat',
910     s   FORM='UNFORMATTED',
911     s   ACCESS='DIRECT'
912     s  ,RECL=4*(ifd(if)-iid(if)+1)*(jfd(if)-jid(if)+1))
913 
914      print*,'Open de dat ok'
915 
916      lmd(if)=lm
917      do l=1,lm
918         zd(l,if)=z(l)*fz
919      enddo
920 
921      irec(if)=0
922!CR
923!      print*,if,imd(if),jmd(if),lmd(if)
924!      print*,'if,imd(if),jmd(if),lmd(if)'
925 
926      return
927      end
928      SUBROUTINE gr_dyn_fi(nfield,im,jm,ngrid,pdyn,pfi)
929      IMPLICIT NONE
930!=======================================================================
931!   passage d'un champ de la grille scalaire a la grille physique
932!=======================================================================
933 
934!-----------------------------------------------------------------------
935!   declarations:
936!   -------------
937 
938      INTEGER im,jm,ngrid,nfield
939      REAL pdyn(im,jm,nfield)
940      REAL pfi(ngrid,nfield)
941 
942      INTEGER j,ifield,ig
943 
944!-----------------------------------------------------------------------
945!   calcul:
946!   -------
947 
948      IF(ngrid.NE.2+(jm-2)*(im-1).AND.ngrid.NE.1)
949     s    STOP 'probleme de dim'
950!   traitement des poles
951      CALL SCOPY(nfield,pdyn,im*jm,pfi,ngrid)
952      CALL SCOPY(nfield,pdyn(1,jm,1),im*jm,pfi(ngrid,1),ngrid)
953 
954!   traitement des point normaux
955      DO ifield=1,nfield
956         DO j=2,jm-1
957            ig=2+(j-2)*(im-1)
958            CALL SCOPY(im-1,pdyn(1,j,ifield),1,pfi(ig,ifield),1)
959         ENDDO
960      ENDDO
961 
962      RETURN
963      END
964 
965!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
966      subroutine writelim
967     s   (phy_nat,phy_alb,phy_sst,phy_bil,phy_rug,phy_ice,
968     s    phy_fter,phy_foce,phy_flic,phy_fsic)
969
970      use dimphy
971      implicit none
972!
973#include "dimensions.h"
974#include "netcdf.inc"
975 
976      integer klon1d
977      integer k
978      parameter (klon1d=1)
979      REAL phy_nat(klon1d,360)
980      REAL phy_alb(klon1d,360)
981      REAL phy_sst(klon1d,360)
982      REAL phy_bil(klon1d,360)
983      REAL phy_rug(klon1d,360)
984      REAL phy_ice(klon1d,360)
985      REAL phy_fter(klon1d,360)
986      REAL phy_foce(klon1d,360)
987      REAL phy_flic(klon1d,360)
988      REAL phy_fsic(klon1d,360)
989 
990      INTEGER ierr
991      INTEGER dimfirst(3)
992      INTEGER dimlast(3)
993!
994      INTEGER nid, ndim, ntim
995      INTEGER dims(2), debut(2), epais(2)
996      INTEGER id_tim
997      INTEGER id_NAT, id_SST, id_BILS, id_RUG, id_ALB
998      INTEGER id_FTER,id_FOCE,id_FSIC,id_FLIC
999 
1000      PRINT*, 'Ecriture du fichier limit'
1001!
1002      ierr = NF_CREATE ("limit.nc", NF_CLOBBER, nid)
1003!
1004      ierr = NF_PUT_ATT_TEXT (nid, NF_GLOBAL, "title", 30,
1005     .                       "Fichier conditions aux limites")
1006      ierr = NF_DEF_DIM (nid, "points_physiques", klon1d, ndim)
1007      ierr = NF_DEF_DIM (nid, "time", NF_UNLIMITED, ntim)
1008!
1009      dims(1) = ndim
1010      dims(2) = ntim
1011!
1012!cc      ierr = NF_DEF_VAR (nid, "TEMPS", NF_DOUBLE, 1,ntim, id_tim)
1013      ierr = NF_DEF_VAR (nid, "TEMPS", NF_FLOAT, 1,ntim, id_tim)
1014      ierr = NF_PUT_ATT_TEXT (nid, id_tim, "title", 17,
1015     .                        "Jour dans l annee")
1016!cc      ierr = NF_DEF_VAR (nid, "NAT", NF_DOUBLE, 2,dims, id_NAT)
1017      ierr = NF_DEF_VAR (nid, "NAT", NF_FLOAT, 2,dims, id_NAT)
1018      ierr = NF_PUT_ATT_TEXT (nid, id_NAT, "title", 23,
1019     .                        "Nature du sol (0,1,2,3)")
1020!cc      ierr = NF_DEF_VAR (nid, "SST", NF_DOUBLE, 2,dims, id_SST)
1021      ierr = NF_DEF_VAR (nid, "SST", NF_FLOAT, 2,dims, id_SST)
1022      ierr = NF_PUT_ATT_TEXT (nid, id_SST, "title", 35,
1023     .                        "Temperature superficielle de la mer")
1024!cc      ierr = NF_DEF_VAR (nid, "BILS", NF_DOUBLE, 2,dims, id_BILS)
1025      ierr = NF_DEF_VAR (nid, "BILS", NF_FLOAT, 2,dims, id_BILS)
1026      ierr = NF_PUT_ATT_TEXT (nid, id_BILS, "title", 32,
1027     .                        "Reference flux de chaleur au sol")
1028!cc      ierr = NF_DEF_VAR (nid, "ALB", NF_DOUBLE, 2,dims, id_ALB)
1029      ierr = NF_DEF_VAR (nid, "ALB", NF_FLOAT, 2,dims, id_ALB)
1030      ierr = NF_PUT_ATT_TEXT (nid, id_ALB, "title", 19,
1031     .                        "Albedo a la surface")
1032!cc      ierr = NF_DEF_VAR (nid, "RUG", NF_DOUBLE, 2,dims, id_RUG)
1033      ierr = NF_DEF_VAR (nid, "RUG", NF_FLOAT, 2,dims, id_RUG)
1034      ierr = NF_PUT_ATT_TEXT (nid, id_RUG, "title", 8,
1035     .                        "Rugosite")
1036
1037      ierr = NF_DEF_VAR (nid, "FTER", NF_FLOAT, 2,dims, id_FTER)
1038      ierr = NF_PUT_ATT_TEXT (nid, id_FTER, "title", 8,"Frac. Terre")
1039      ierr = NF_DEF_VAR (nid, "FOCE", NF_FLOAT, 2,dims, id_FOCE)
1040      ierr = NF_PUT_ATT_TEXT (nid, id_FOCE, "title", 8,"Frac. Terre")
1041      ierr = NF_DEF_VAR (nid, "FSIC", NF_FLOAT, 2,dims, id_FSIC)
1042      ierr = NF_PUT_ATT_TEXT (nid, id_FSIC, "title", 8,"Frac. Terre")
1043      ierr = NF_DEF_VAR (nid, "FLIC", NF_FLOAT, 2,dims, id_FLIC)
1044      ierr = NF_PUT_ATT_TEXT (nid, id_FLIC, "title", 8,"Frac. Terre")
1045!
1046      ierr = NF_ENDDEF(nid)
1047!
1048      DO k = 1, 360
1049!
1050      debut(1) = 1
1051      debut(2) = k
1052      epais(1) = klon1d
1053      epais(2) = 1
1054!
1055#ifdef NC_DOUBLE
1056      ierr = NF_PUT_VAR1_DOUBLE (nid,id_tim,k,DBLE(k))
1057      ierr = NF_PUT_VARA_DOUBLE (nid,id_NAT,debut,epais,phy_nat(1,k))
1058      ierr = NF_PUT_VARA_DOUBLE (nid,id_SST,debut,epais,phy_sst(1,k))
1059      ierr = NF_PUT_VARA_DOUBLE (nid,id_BILS,debut,epais,phy_bil(1,k))
1060      ierr = NF_PUT_VARA_DOUBLE (nid,id_ALB,debut,epais,phy_alb(1,k))
1061      ierr = NF_PUT_VARA_DOUBLE (nid,id_RUG,debut,epais,phy_rug(1,k))
1062      ierr = NF_PUT_VARA_DOUBLE (nid,id_FTER,debut,epais,phy_fter(1,k))
1063      ierr = NF_PUT_VARA_DOUBLE (nid,id_FOCE,debut,epais,phy_foce(1,k))
1064      ierr = NF_PUT_VARA_DOUBLE (nid,id_FSIC,debut,epais,phy_fsic(1,k))
1065      ierr = NF_PUT_VARA_DOUBLE (nid,id_FLIC,debut,epais,phy_flic(1,k))
1066#else
1067      ierr = NF_PUT_VAR1_REAL (nid,id_tim,k,FLOAT(k))
1068      ierr = NF_PUT_VARA_REAL (nid,id_NAT,debut,epais,phy_nat(1,k))
1069      ierr = NF_PUT_VARA_REAL (nid,id_SST,debut,epais,phy_sst(1,k))
1070      ierr = NF_PUT_VARA_REAL (nid,id_BILS,debut,epais,phy_bil(1,k))
1071      ierr = NF_PUT_VARA_REAL (nid,id_ALB,debut,epais,phy_alb(1,k))
1072      ierr = NF_PUT_VARA_REAL (nid,id_RUG,debut,epais,phy_rug(1,k))
1073      ierr = NF_PUT_VARA_REAL (nid,id_FTER,debut,epais,phy_fter(1,k))
1074      ierr = NF_PUT_VARA_REAL (nid,id_FOCE,debut,epais,phy_foce(1,k))
1075      ierr = NF_PUT_VARA_REAL (nid,id_FSIC,debut,epais,phy_fsic(1,k))
1076      ierr = NF_PUT_VARA_REAL (nid,id_FLIC,debut,epais,phy_flic(1,k))
1077
1078#endif
1079!
1080      ENDDO
1081!
1082      ierr = NF_CLOSE(nid)
1083!
1084      return
1085      end
1086 
1087!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1088      SUBROUTINE disvert0(pa,preff,ap,bp,dpres,presnivs,nivsigs,nivsig)
1089 
1090!    Auteur :  P. Le Van .
1091!
1092      IMPLICIT NONE
1093 
1094#include "dimensions.h"
1095#include "paramet.h"
1096!
1097!=======================================================================
1098!
1099!
1100!    s = sigma ** kappa   :  coordonnee  verticale
1101!    dsig(l)            : epaisseur de la couche l ds la coord.  s
1102!    sig(l)             : sigma a l'interface des couches l et l-1
1103!    ds(l)              : distance entre les couches l et l-1 en coord.s
1104!
1105!=======================================================================
1106!
1107      REAL pa,preff
1108      REAL ap(llmp1),bp(llmp1),dpres(llm),nivsigs(llm),nivsig(llmp1)
1109      REAL presnivs(llm)
1110!
1111!   declarations:
1112!   -------------
1113!
1114      REAL sig(llm+1),dsig(llm)
1115!
1116      INTEGER l
1117      REAL snorm
1118      REAL alpha,beta,gama,delta,deltaz,h
1119      INTEGER np,ierr
1120      REAL pi,x
1121 
1122!-----------------------------------------------------------------------
1123!
1124      pi=2.*ASIN(1.)
1125 
1126      OPEN(99,file='sigma.def',status='old',form='formatted',
1127     s   iostat=ierr)
1128 
1129!-----------------------------------------------------------------------
1130!   cas 1 on lit les options dans sigma.def:
1131!   ----------------------------------------
1132 
1133      IF (ierr.eq.0) THEN
1134 
1135      print*,'WARNING!!! on lit les options dans sigma.def'
1136      READ(99,*) deltaz
1137      READ(99,*) h
1138      READ(99,*) beta
1139      READ(99,*) gama
1140      READ(99,*) delta
1141      READ(99,*) np
1142      CLOSE(99)
1143      alpha=deltaz/(llm*h)
1144!
1145 
1146       DO 1  l = 1, llm
1147       dsig(l) = (alpha+(1.-alpha)*exp(-beta*(llm-l)))*
1148     $          ( (tanh(gama*l)/tanh(gama*llm))**np +
1149     $            (1.-l/FLOAT(llm))*delta )
1150   1   CONTINUE
1151 
1152       sig(1)=1.
1153       DO 101 l=1,llm-1
1154          sig(l+1)=sig(l)*(1.-dsig(l))/(1.+dsig(l))
1155101    CONTINUE
1156       sig(llm+1)=0.
1157 
1158       DO 2  l = 1, llm
1159       dsig(l) = sig(l)-sig(l+1)
1160   2   CONTINUE
1161!
1162 
1163      ELSE
1164!-----------------------------------------------------------------------
1165!   cas 2 ancienne discretisation (LMD5...):
1166!   ----------------------------------------
1167 
1168      PRINT*,'WARNING!!! Ancienne discretisation verticale'
1169 
1170      h=7.
1171      snorm  = 0.
1172      DO l = 1, llm
1173         x = 2.*asin(1.) * (FLOAT(l)-0.5) / float(llm+1)
1174         dsig(l) = 1.0 + 7.0 * SIN(x)**2
1175         snorm = snorm + dsig(l)
1176      ENDDO
1177      snorm = 1./snorm
1178      DO l = 1, llm
1179         dsig(l) = dsig(l)*snorm
1180      ENDDO
1181      sig(llm+1) = 0.
1182      DO l = llm, 1, -1
1183         sig(l) = sig(l+1) + dsig(l)
1184      ENDDO
1185 
1186      ENDIF
1187 
1188 
1189      DO l=1,llm
1190        nivsigs(l) = FLOAT(l)
1191      ENDDO
1192 
1193      DO l=1,llmp1
1194        nivsig(l)= FLOAT(l)
1195      ENDDO
1196 
1197!
1198!    ....  Calculs  de ap(l) et de bp(l)  ....
1199!    .........................................
1200!
1201!
1202!   .....  pa et preff sont lus  sur les fichiers start par lectba  .....
1203!
1204 
1205      bp(llmp1) =   0.
1206 
1207      DO l = 1, llm
1208!c
1209!cc    ap(l) = 0.
1210!cc    bp(l) = sig(l)
1211 
1212      bp(l) = EXP( 1. -1./( sig(l)*sig(l)) )
1213      ap(l) = pa * ( sig(l) - bp(l) )
1214!
1215      ENDDO
1216      ap(llmp1) = pa * ( sig(llmp1) - bp(llmp1) )
1217 
1218      PRINT *,' BP '
1219      PRINT *,  bp
1220      PRINT *,' AP '
1221      PRINT *,  ap
1222 
1223      DO l = 1, llm
1224       dpres(l) = bp(l) - bp(l+1)
1225       presnivs(l) = 0.5 *( ap(l)+bp(l)*preff + ap(l+1)+bp(l+1)*preff )
1226      ENDDO
1227 
1228      PRINT *,' PRESNIVS '
1229      PRINT *,presnivs
1230 
1231      RETURN
1232      END
1233
1234!======================================================================
1235       SUBROUTINE read_tsurf1d(knon,knindex,sst_out)
1236
1237! This subroutine specifies the surface temperature to be used in 1D simulations
1238
1239      USE dimphy, ONLY : klon
1240
1241      INTEGER, INTENT(IN)                  :: knon     ! nomber of points on compressed grid
1242      INTEGER, DIMENSION(klon), INTENT(IN) :: knindex  ! grid point number for compressed grid
1243      REAL, DIMENSION(klon), INTENT(OUT)   :: sst_out  ! tsurf used to force the single-column model
1244
1245       INTEGER :: i
1246! COMMON defined in lmdz1d.F:
1247       real ts_cur
1248       common /sst_forcing/ts_cur
1249
1250       DO i = 1, knon
1251        sst_out(i) = ts_cur
1252       ENDDO
1253
1254      END SUBROUTINE read_tsurf1d
1255
1256!===============================================================
1257      subroutine advect_vert(llm,w,dt,q,plev)
1258!===============================================================
1259!   Schema amont pour l'advection verticale en 1D
1260!   w est la vitesse verticale dp/dt en Pa/s
1261!   Traitement en volumes finis
1262!   d / dt ( zm q ) = delta_z ( omega q )
1263!   d / dt ( zm ) = delta_z ( omega )
1264!   avec zm = delta_z ( p )
1265!   si * designe la valeur au pas de temps t+dt
1266!   zm*(l) q*(l) - zm(l) q(l) = w(l+1) q(l+1) - w(l) q(l)
1267!   zm*(l) -zm(l) = w(l+1) - w(l)
1268!   avec w=omega * dt
1269!---------------------------------------------------------------
1270      implicit none
1271! arguments
1272      integer llm
1273      real w(llm+1),q(llm),plev(llm+1),dt
1274
1275! local
1276      integer l
1277      real zwq(llm+1),zm(llm+1),zw(llm+1)
1278      real qold
1279
1280!---------------------------------------------------------------
1281
1282      do l=1,llm
1283         zw(l)=dt*w(l)
1284         zm(l)=plev(l)-plev(l+1)
1285         zwq(l)=q(l)*zw(l)
1286      enddo
1287      zwq(llm+1)=0.
1288      zw(llm+1)=0.
1289 
1290      do l=1,llm
1291         qold=q(l)
1292         q(l)=(q(l)*zm(l)+zwq(l+1)-zwq(l))/(zm(l)+zw(l+1)-zw(l))
1293         print*,'ADV Q ',zm(l),zw(l),zwq(l),qold,q(l)
1294      enddo
1295
1296 
1297      return
1298      end
1299
1300!===============================================================
1301
1302
1303       SUBROUTINE advect_va(llm,omega,d_t_va,d_q_va,d_u_va,d_v_va,
1304     !                q,temp,u,v,
1305     !            play,plev)
1306!itlmd
1307!----------------------------------------------------------------------
1308!   Calcul de l'advection verticale (ascendance et subsidence) de
1309!   température et d'humidité. Hypothèse : ce qui rentre de l'extérieur
1310!   a les mêmes caractéristiques que l'air de la colonne 1D (WTG) ou
1311!   sans WTG rajouter une advection horizontale
1312!---------------------------------------------------------------------- 
1313        implicit none
1314#include "YOMCST.h"
1315!        argument
1316        integer llm
1317        real  omega(llm+1),d_t_va(llm), d_q_va(llm,3)
1318        real  d_u_va(llm), d_v_va(llm)
1319        real  q(llm,3),temp(llm)
1320        real  u(llm),v(llm)
1321        real  play(llm),plev(llm+1)
1322! interne
1323        integer l
1324        real alpha,omgdown,omgup
1325
1326      do l= 1,llm
1327       if(l.eq.1) then
1328!si omgup pour la couche 1, alors tendance nulle
1329        omgdown=max(omega(2),0.0)
1330        alpha = rkappa*temp(l)*(1.+q(l,1)*rv/rd)/(play(l)*
1331     &           (1.+q(l,1)))
1332        d_t_va(l)= alpha*(omgdown)-
1333     &              omgdown*(temp(l)-temp(l+1))
1334     &              /(play(l)-play(l+1))             
1335
1336        d_q_va(l,:)= -omgdown*(q(l,:)-q(l+1,:))
1337     &              /(play(l)-play(l+1))             
1338
1339        d_u_va(l)= -omgdown*(u(l)-u(l+1))
1340     &              /(play(l)-play(l+1))             
1341        d_v_va(l)= -omgdown*(v(l)-v(l+1))
1342     &              /(play(l)-play(l+1))             
1343
1344       
1345       elseif(l.eq.llm) then
1346        omgup=min(omega(l),0.0)
1347        alpha = rkappa*temp(l)*(1.+q(l,1)*rv/rd)/(play(l)*
1348     &           (1.+q(l,1)))
1349        d_t_va(l)= alpha*(omgup)-
1350!bug?     &              omgup*(temp(l-1)-temp(l))/(play(l-1)-plev(l))
1351     &              omgup*(temp(l-1)-temp(l))/(play(l-1)-play(l))
1352        d_q_va(l,:)= -omgup*(q(l-1,:)-q(l,:))/(play(l-1)-play(l))
1353        d_u_va(l)= -omgup*(u(l-1)-u(l))/(play(l-1)-play(l))
1354        d_v_va(l)= -omgup*(v(l-1)-v(l))/(play(l-1)-play(l))
1355       
1356       else
1357        omgup=min(omega(l),0.0)
1358        omgdown=max(omega(l+1),0.0)
1359        alpha = rkappa*temp(l)*(1.+q(l,1)*rv/rd)/(play(l)*
1360     &           (1.+q(l,1)))
1361        d_t_va(l)= alpha*(omgup+omgdown)-
1362     &              omgdown*(temp(l)-temp(l+1))
1363     &              /(play(l)-play(l+1))-             
1364!bug?     &              omgup*(temp(l-1)-temp(l))/(play(l-1)-plev(l))
1365     &              omgup*(temp(l-1)-temp(l))/(play(l-1)-play(l))
1366!      print*, '  ??? '
1367
1368        d_q_va(l,:)= -omgdown*(q(l,:)-q(l+1,:))
1369     &              /(play(l)-play(l+1))-             
1370     &              omgup*(q(l-1,:)-q(l,:))/(play(l-1)-play(l))
1371        d_u_va(l)= -omgdown*(u(l)-u(l+1))
1372     &              /(play(l)-play(l+1))-             
1373     &              omgup*(u(l-1)-u(l))/(play(l-1)-play(l))
1374        d_v_va(l)= -omgdown*(v(l)-v(l+1))
1375     &              /(play(l)-play(l+1))-             
1376     &              omgup*(v(l-1)-v(l))/(play(l-1)-play(l))
1377       
1378      endif
1379         
1380      enddo
1381!fin itlmd
1382        return
1383        end
1384!       SUBROUTINE lstendH(llm,omega,d_t_va,d_q_va,d_u_va,d_v_va,
1385       SUBROUTINE lstendH(llm,nqtot,omega,d_t_va,d_q_va,
1386     !                q,temp,u,v,play)
1387!itlmd
1388!----------------------------------------------------------------------
1389!   Calcul de l'advection verticale (ascendance et subsidence) de
1390!   température et d'humidité. Hypothèse : ce qui rentre de l'extérieur
1391!   a les mêmes caractéristiques que l'air de la colonne 1D (WTG) ou
1392!   sans WTG rajouter une advection horizontale
1393!---------------------------------------------------------------------- 
1394        implicit none
1395#include "YOMCST.h"
1396!        argument
1397        integer llm,nqtot
1398        real  omega(llm+1),d_t_va(llm), d_q_va(llm,nqtot)
1399!        real  d_u_va(llm), d_v_va(llm)
1400        real  q(llm,nqtot),temp(llm)
1401        real  u(llm),v(llm)
1402        real  play(llm)
1403        real cor(llm)
1404!        real dph(llm),dudp(llm),dvdp(llm),dqdp(llm),dtdp(llm)
1405        real dph(llm),dqdp(llm),dtdp(llm)
1406! interne
1407        integer l,k
1408        real alpha,omdn,omup
1409
1410!        dudp=0.
1411!        dvdp=0.
1412        dqdp=0.
1413        dtdp=0.
1414!        d_u_va=0.
1415!        d_v_va=0.
1416
1417      cor(:) = rkappa*temp*(1.+q(:,1)*rv/rd)/(play*(1.+q(:,1)))
1418
1419
1420      do k=2,llm-1
1421
1422       dph  (k-1) = (play(k  )- play(k-1  ))
1423!       dudp (k-1) = (u   (k  )- u   (k-1  ))/dph(k-1)
1424!       dvdp (k-1) = (v   (k  )- v   (k-1  ))/dph(k-1)
1425       dqdp (k-1) = (q   (k,1)- q   (k-1,1))/dph(k-1)
1426       dtdp (k-1) = (temp(k  )- temp(k-1  ))/dph(k-1)
1427
1428      enddo
1429
1430!      dudp (  llm  ) = dudp ( llm-1 )
1431!      dvdp (  llm  ) = dvdp ( llm-1 )
1432      dqdp (  llm  ) = dqdp ( llm-1 )
1433      dtdp (  llm  ) = dtdp ( llm-1 )
1434
1435      do k=2,llm-1
1436      omdn=max(0.0,omega(k+1))
1437      omup=min(0.0,omega( k ))
1438
1439!      d_u_va(k)  = -omdn*dudp(k)-omup*dudp(k-1)
1440!      d_v_va(k)  = -omdn*dvdp(k)-omup*dvdp(k-1)
1441      d_q_va(k,1)= -omdn*dqdp(k)-omup*dqdp(k-1)
1442      d_t_va(k)  = -omdn*dtdp(k)-omup*dtdp(k-1)
1443     :              +(omup+omdn)*cor(k)
1444      enddo
1445
1446      omdn=max(0.0,omega( 2 ))
1447      omup=min(0.0,omega(llm))
1448!      d_u_va( 1 )   = -omdn*dudp( 1 )
1449!      d_u_va(llm)   = -omup*dudp(llm)
1450!      d_v_va( 1 )   = -omdn*dvdp( 1 )
1451!      d_v_va(llm)   = -omup*dvdp(llm)
1452      d_q_va( 1 ,1) = -omdn*dqdp( 1 )
1453      d_q_va(llm,1) = -omup*dqdp(llm)
1454      d_t_va( 1 )   = -omdn*dtdp( 1 )+omdn*cor( 1 )
1455      d_t_va(llm)   = -omup*dtdp(llm)!+omup*cor(llm)
1456
1457!      if(abs(rlat(1))>10.) then
1458!     Calculate the tendency due agestrophic motions
1459!      du_age = fcoriolis*(v-vg)
1460!      dv_age = fcoriolis*(ug-u)
1461!      endif
1462
1463!       call writefield_phy('d_t_va',d_t_va,llm)
1464
1465          return
1466         end
1467
1468!======================================================================
1469      SUBROUTINE read_togacoare(fich_toga,nlev_toga,nt_toga
1470     :             ,ts_toga,plev_toga,t_toga,q_toga,u_toga,v_toga,w_toga
1471     :             ,ht_toga,vt_toga,hq_toga,vq_toga)
1472      implicit none
1473
1474c-------------------------------------------------------------------------
1475c Read TOGA-COARE forcing data
1476c-------------------------------------------------------------------------
1477
1478      integer nlev_toga,nt_toga
1479      real ts_toga(nt_toga),plev_toga(nlev_toga,nt_toga)
1480      real t_toga(nlev_toga,nt_toga),q_toga(nlev_toga,nt_toga)
1481      real u_toga(nlev_toga,nt_toga),v_toga(nlev_toga,nt_toga)
1482      real w_toga(nlev_toga,nt_toga)
1483      real ht_toga(nlev_toga,nt_toga),vt_toga(nlev_toga,nt_toga)
1484      real hq_toga(nlev_toga,nt_toga),vq_toga(nlev_toga,nt_toga)
1485      character*80 fich_toga
1486
1487      integer no,l,k,ip
1488      real riy,rim,rid,rih,bid
1489
1490      integer iy,im,id,ih
1491     
1492       real plev_min
1493
1494       plev_min = 55.  ! pas de tendance de vap. d eau au-dessus de 55 hPa
1495
1496      open(21,file=trim(fich_toga),form='formatted')
1497      read(21,'(a)')
1498      do ip = 1, nt_toga
1499      read(21,'(a)')
1500      read(21,'(a)')
1501      read(21,223) iy, im, id, ih, bid, ts_toga(ip), bid,bid,bid,bid
1502      read(21,'(a)')
1503      read(21,'(a)')
1504
1505       do k = 1, nlev_toga
1506         read(21,230) plev_toga(k,ip), t_toga(k,ip), q_toga(k,ip)
1507     :       ,u_toga(k,ip), v_toga(k,ip), w_toga(k,ip)
1508     :       ,ht_toga(k,ip), vt_toga(k,ip), hq_toga(k,ip), vq_toga(k,ip)
1509
1510! conversion in SI units:
1511         t_toga(k,ip)=t_toga(k,ip)+273.15     ! K
1512         q_toga(k,ip)=q_toga(k,ip)*0.001      ! kg/kg
1513         w_toga(k,ip)=w_toga(k,ip)*100./3600. ! Pa/s
1514! no water vapour tendency above 55 hPa
1515         if (plev_toga(k,ip) .lt. plev_min) then
1516          q_toga(k,ip) = 0.
1517          hq_toga(k,ip) = 0.
1518          vq_toga(k,ip) =0.
1519         endif
1520       enddo
1521
1522         ts_toga(ip)=ts_toga(ip)+273.15       ! K
1523       enddo
1524       close(21)
1525
1526  223 format(4i3,6f8.2)
1527  226 format(f7.1,1x,10f8.2)
1528  227 format(f7.1,1x,1p,4e11.3)
1529  230 format(6f9.3,4e11.3)
1530
1531          return
1532          end
1533
1534!=====================================================================
1535      subroutine read_twpice(fich_twpice,nlevel,ntime
1536     :     ,T_srf,plev,T,q,u,v,omega
1537     :     ,T_adv_h,T_adv_v,q_adv_h,q_adv_v)
1538
1539!program reading forcings of the TWP-ICE experiment
1540
1541!      use netcdf
1542
1543      implicit none
1544
1545#include "netcdf.inc"
1546
1547      integer ntime,nlevel
1548      integer l,k
1549      character*80 :: fich_twpice
1550      real*8 time(ntime)
1551      real*8 lat, lon, alt, phis       
1552      real*8 lev(nlevel)
1553      real*8 plev(nlevel,ntime)
1554
1555      real*8 T(nlevel,ntime)
1556      real*8 q(nlevel,ntime),u(nlevel,ntime)
1557      real*8 v(nlevel,ntime)
1558      real*8 omega(nlevel,ntime), div(nlevel,ntime)
1559      real*8 T_adv_h(nlevel,ntime)
1560      real*8 T_adv_v(nlevel,ntime), q_adv_h(nlevel,ntime)
1561      real*8 q_adv_v(nlevel,ntime)     
1562      real*8 s(nlevel,ntime), s_adv_h(nlevel,ntime)
1563      real*8 s_adv_v(nlevel,ntime)
1564      real*8 p_srf_aver(ntime), p_srf_center(ntime)
1565      real*8 T_srf(ntime)
1566
1567      integer nid, ierr
1568      integer nbvar3d
1569      parameter(nbvar3d=20)
1570      integer var3didin(nbvar3d)
1571
1572      ierr = NF_OPEN(fich_twpice,NF_NOWRITE,nid)
1573      if (ierr.NE.NF_NOERR) then
1574         write(*,*) 'ERROR: Pb opening forcings cdf file '
1575         write(*,*) NF_STRERROR(ierr)
1576         stop ""
1577      endif
1578
1579      ierr=NF_INQ_VARID(nid,"lat",var3didin(1))
1580         if(ierr/=NF_NOERR) then
1581           write(*,*) NF_STRERROR(ierr)
1582           stop 'lat'
1583         endif
1584     
1585       ierr=NF_INQ_VARID(nid,"lon",var3didin(2))
1586         if(ierr/=NF_NOERR) then
1587           write(*,*) NF_STRERROR(ierr)
1588           stop 'lon'
1589         endif
1590
1591       ierr=NF_INQ_VARID(nid,"alt",var3didin(3))
1592         if(ierr/=NF_NOERR) then
1593           write(*,*) NF_STRERROR(ierr)
1594           stop 'alt'
1595         endif
1596
1597      ierr=NF_INQ_VARID(nid,"phis",var3didin(4))
1598         if(ierr/=NF_NOERR) then
1599           write(*,*) NF_STRERROR(ierr)
1600           stop 'phis'
1601         endif
1602
1603      ierr=NF_INQ_VARID(nid,"T",var3didin(5))
1604         if(ierr/=NF_NOERR) then
1605           write(*,*) NF_STRERROR(ierr)
1606           stop 'T'
1607         endif
1608
1609      ierr=NF_INQ_VARID(nid,"q",var3didin(6))
1610         if(ierr/=NF_NOERR) then
1611           write(*,*) NF_STRERROR(ierr)
1612           stop 'q'
1613         endif
1614
1615      ierr=NF_INQ_VARID(nid,"u",var3didin(7))
1616         if(ierr/=NF_NOERR) then
1617           write(*,*) NF_STRERROR(ierr)
1618           stop 'u'
1619         endif
1620
1621      ierr=NF_INQ_VARID(nid,"v",var3didin(8))
1622         if(ierr/=NF_NOERR) then
1623           write(*,*) NF_STRERROR(ierr)
1624           stop 'v'
1625         endif
1626
1627      ierr=NF_INQ_VARID(nid,"omega",var3didin(9))
1628         if(ierr/=NF_NOERR) then
1629           write(*,*) NF_STRERROR(ierr)
1630           stop 'omega'
1631         endif
1632
1633      ierr=NF_INQ_VARID(nid,"div",var3didin(10))
1634         if(ierr/=NF_NOERR) then
1635           write(*,*) NF_STRERROR(ierr)
1636           stop 'div'
1637         endif
1638
1639      ierr=NF_INQ_VARID(nid,"T_adv_h",var3didin(11))
1640         if(ierr/=NF_NOERR) then
1641           write(*,*) NF_STRERROR(ierr)
1642           stop 'T_adv_h'
1643         endif
1644
1645      ierr=NF_INQ_VARID(nid,"T_adv_v",var3didin(12))
1646         if(ierr/=NF_NOERR) then
1647           write(*,*) NF_STRERROR(ierr)
1648           stop 'T_adv_v'
1649         endif
1650
1651      ierr=NF_INQ_VARID(nid,"q_adv_h",var3didin(13))
1652         if(ierr/=NF_NOERR) then
1653           write(*,*) NF_STRERROR(ierr)
1654           stop 'q_adv_h'
1655         endif
1656
1657      ierr=NF_INQ_VARID(nid,"q_adv_v",var3didin(14))
1658         if(ierr/=NF_NOERR) then
1659           write(*,*) NF_STRERROR(ierr)
1660           stop 'q_adv_v'
1661         endif
1662
1663      ierr=NF_INQ_VARID(nid,"s",var3didin(15))
1664         if(ierr/=NF_NOERR) then
1665           write(*,*) NF_STRERROR(ierr)
1666           stop 's'
1667         endif
1668
1669      ierr=NF_INQ_VARID(nid,"s_adv_h",var3didin(16))
1670         if(ierr/=NF_NOERR) then
1671           write(*,*) NF_STRERROR(ierr)
1672           stop 's_adv_h'
1673         endif
1674   
1675      ierr=NF_INQ_VARID(nid,"s_adv_v",var3didin(17))
1676         if(ierr/=NF_NOERR) then
1677           write(*,*) NF_STRERROR(ierr)
1678           stop 's_adv_v'
1679         endif
1680
1681      ierr=NF_INQ_VARID(nid,"p_srf_aver",var3didin(18))
1682         if(ierr/=NF_NOERR) then
1683           write(*,*) NF_STRERROR(ierr)
1684           stop 'p_srf_aver'
1685         endif
1686
1687      ierr=NF_INQ_VARID(nid,"p_srf_center",var3didin(19))
1688         if(ierr/=NF_NOERR) then
1689           write(*,*) NF_STRERROR(ierr)
1690           stop 'p_srf_center'
1691         endif
1692
1693      ierr=NF_INQ_VARID(nid,"T_srf",var3didin(20))
1694         if(ierr/=NF_NOERR) then
1695           write(*,*) NF_STRERROR(ierr)
1696           stop 'T_srf'
1697         endif
1698
1699!dimensions lecture
1700      call catchaxis(nid,ntime,nlevel,time,lev,ierr)
1701
1702!pressure
1703       do l=1,ntime
1704       do k=1,nlevel
1705          plev(k,l)=lev(k)
1706       enddo
1707       enddo
1708         
1709#ifdef NC_DOUBLE
1710         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(1),lat)
1711#else
1712         ierr = NF_GET_VAR_REAL(nid,var3didin(1),lat)
1713#endif
1714         if(ierr/=NF_NOERR) then
1715            write(*,*) NF_STRERROR(ierr)
1716            stop "getvarup"
1717         endif
1718!         write(*,*)'lecture lat ok',lat
1719
1720#ifdef NC_DOUBLE
1721         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(2),lon)
1722#else
1723         ierr = NF_GET_VAR_REAL(nid,var3didin(2),lon)
1724#endif
1725         if(ierr/=NF_NOERR) then
1726            write(*,*) NF_STRERROR(ierr)
1727            stop "getvarup"
1728         endif
1729!         write(*,*)'lecture lon ok',lon
1730 
1731#ifdef NC_DOUBLE
1732         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(3),alt)
1733#else
1734         ierr = NF_GET_VAR_REAL(nid,var3didin(3),alt)
1735#endif
1736         if(ierr/=NF_NOERR) then
1737            write(*,*) NF_STRERROR(ierr)
1738            stop "getvarup"
1739         endif
1740!          write(*,*)'lecture alt ok',alt
1741 
1742#ifdef NC_DOUBLE
1743         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(4),phis)
1744#else
1745         ierr = NF_GET_VAR_REAL(nid,var3didin(4),phis)
1746#endif
1747         if(ierr/=NF_NOERR) then
1748            write(*,*) NF_STRERROR(ierr)
1749            stop "getvarup"
1750         endif
1751!          write(*,*)'lecture phis ok',phis
1752         
1753#ifdef NC_DOUBLE
1754         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(5),T)
1755#else
1756         ierr = NF_GET_VAR_REAL(nid,var3didin(5),T)
1757#endif
1758         if(ierr/=NF_NOERR) then
1759            write(*,*) NF_STRERROR(ierr)
1760            stop "getvarup"
1761         endif
1762!         write(*,*)'lecture T ok'
1763
1764#ifdef NC_DOUBLE
1765         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(6),q)
1766#else
1767         ierr = NF_GET_VAR_REAL(nid,var3didin(6),q)
1768#endif
1769         if(ierr/=NF_NOERR) then
1770            write(*,*) NF_STRERROR(ierr)
1771            stop "getvarup"
1772         endif
1773!         write(*,*)'lecture q ok'
1774!q in kg/kg
1775       do l=1,ntime
1776       do k=1,nlevel
1777          q(k,l)=q(k,l)/1000.
1778       enddo
1779       enddo
1780#ifdef NC_DOUBLE
1781         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(7),u)
1782#else
1783         ierr = NF_GET_VAR_REAL(nid,var3didin(7),u)
1784#endif
1785         if(ierr/=NF_NOERR) then
1786            write(*,*) NF_STRERROR(ierr)
1787            stop "getvarup"
1788         endif
1789!         write(*,*)'lecture u ok'
1790
1791#ifdef NC_DOUBLE
1792         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(8),v)
1793#else
1794         ierr = NF_GET_VAR_REAL(nid,var3didin(8),v)
1795#endif
1796         if(ierr/=NF_NOERR) then
1797            write(*,*) NF_STRERROR(ierr)
1798            stop "getvarup"
1799         endif
1800!         write(*,*)'lecture v ok'
1801
1802#ifdef NC_DOUBLE
1803         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(9),omega)
1804#else
1805         ierr = NF_GET_VAR_REAL(nid,var3didin(9),omega)
1806#endif
1807         if(ierr/=NF_NOERR) then
1808            write(*,*) NF_STRERROR(ierr)
1809            stop "getvarup"
1810         endif
1811!         write(*,*)'lecture omega ok'
1812!omega in mb/hour
1813       do l=1,ntime
1814       do k=1,nlevel
1815          omega(k,l)=omega(k,l)*100./3600.
1816       enddo
1817       enddo
1818
1819#ifdef NC_DOUBLE
1820         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(10),div)
1821#else
1822         ierr = NF_GET_VAR_REAL(nid,var3didin(10),div)
1823#endif
1824         if(ierr/=NF_NOERR) then
1825            write(*,*) NF_STRERROR(ierr)
1826            stop "getvarup"
1827         endif
1828!         write(*,*)'lecture div ok'
1829
1830#ifdef NC_DOUBLE
1831         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(11),T_adv_h)
1832#else
1833         ierr = NF_GET_VAR_REAL(nid,var3didin(11),T_adv_h)
1834#endif
1835         if(ierr/=NF_NOERR) then
1836            write(*,*) NF_STRERROR(ierr)
1837            stop "getvarup"
1838         endif
1839!         write(*,*)'lecture T_adv_h ok'
1840!T adv in K/s
1841       do l=1,ntime
1842       do k=1,nlevel
1843          T_adv_h(k,l)=T_adv_h(k,l)/3600.
1844       enddo
1845       enddo
1846
1847
1848#ifdef NC_DOUBLE
1849         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(12),T_adv_v)
1850#else
1851         ierr = NF_GET_VAR_REAL(nid,var3didin(12),T_adv_v)
1852#endif
1853         if(ierr/=NF_NOERR) then
1854            write(*,*) NF_STRERROR(ierr)
1855            stop "getvarup"
1856         endif
1857!         write(*,*)'lecture T_adv_v ok'
1858!T adv in K/s
1859       do l=1,ntime
1860       do k=1,nlevel
1861          T_adv_v(k,l)=T_adv_v(k,l)/3600.
1862       enddo
1863       enddo
1864
1865#ifdef NC_DOUBLE
1866         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(13),q_adv_h)
1867#else
1868         ierr = NF_GET_VAR_REAL(nid,var3didin(13),q_adv_h)
1869#endif
1870         if(ierr/=NF_NOERR) then
1871            write(*,*) NF_STRERROR(ierr)
1872            stop "getvarup"
1873         endif
1874!         write(*,*)'lecture q_adv_h ok'
1875!q adv in kg/kg/s
1876       do l=1,ntime
1877       do k=1,nlevel
1878          q_adv_h(k,l)=q_adv_h(k,l)/1000./3600.
1879       enddo
1880       enddo
1881
1882
1883#ifdef NC_DOUBLE
1884         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(14),q_adv_v)
1885#else
1886         ierr = NF_GET_VAR_REAL(nid,var3didin(14),q_adv_v)
1887#endif
1888         if(ierr/=NF_NOERR) then
1889            write(*,*) NF_STRERROR(ierr)
1890            stop "getvarup"
1891         endif
1892!         write(*,*)'lecture q_adv_v ok'
1893!q adv in kg/kg/s
1894       do l=1,ntime
1895       do k=1,nlevel
1896          q_adv_v(k,l)=q_adv_v(k,l)/1000./3600.
1897       enddo
1898       enddo
1899
1900
1901#ifdef NC_DOUBLE
1902         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(15),s)
1903#else
1904         ierr = NF_GET_VAR_REAL(nid,var3didin(15),s)
1905#endif
1906         if(ierr/=NF_NOERR) then
1907            write(*,*) NF_STRERROR(ierr)
1908            stop "getvarup"
1909         endif
1910
1911#ifdef NC_DOUBLE
1912         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(16),s_adv_h)
1913#else
1914         ierr = NF_GET_VAR_REAL(nid,var3didin(16),s_adv_h)
1915#endif
1916         if(ierr/=NF_NOERR) then
1917            write(*,*) NF_STRERROR(ierr)
1918            stop "getvarup"
1919         endif
1920
1921#ifdef NC_DOUBLE
1922         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(17),s_adv_v)
1923#else
1924         ierr = NF_GET_VAR_REAL(nid,var3didin(17),s_adv_v)
1925#endif
1926         if(ierr/=NF_NOERR) then
1927            write(*,*) NF_STRERROR(ierr)
1928            stop "getvarup"
1929         endif
1930
1931#ifdef NC_DOUBLE
1932         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(18),p_srf_aver)
1933#else
1934         ierr = NF_GET_VAR_REAL(nid,var3didin(18),p_srf_aver)
1935#endif
1936         if(ierr/=NF_NOERR) then
1937            write(*,*) NF_STRERROR(ierr)
1938            stop "getvarup"
1939         endif
1940
1941#ifdef NC_DOUBLE
1942         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(19),p_srf_center)
1943#else
1944         ierr = NF_GET_VAR_REAL(nid,var3didin(19),p_srf_center)
1945#endif
1946         if(ierr/=NF_NOERR) then
1947            write(*,*) NF_STRERROR(ierr)
1948            stop "getvarup"
1949         endif
1950
1951#ifdef NC_DOUBLE
1952         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(20),T_srf)
1953#else
1954         ierr = NF_GET_VAR_REAL(nid,var3didin(20),T_srf)
1955#endif
1956         if(ierr/=NF_NOERR) then
1957            write(*,*) NF_STRERROR(ierr)
1958            stop "getvarup"
1959         endif
1960!         write(*,*)'lecture T_srf ok', T_srf
1961
1962         return
1963         end subroutine read_twpice
1964!=====================================================================
1965         subroutine catchaxis(nid,ttm,llm,time,lev,ierr)
1966
1967!         use netcdf
1968
1969         implicit none
1970#include "netcdf.inc"
1971         integer nid,ttm,llm
1972         real*8 time(ttm)
1973         real*8 lev(llm)
1974         integer ierr
1975
1976         integer i
1977         integer timevar,levvar
1978         integer timelen,levlen
1979         integer timedimin,levdimin
1980
1981! Control & lecture on dimensions
1982! ===============================
1983         ierr=NF_INQ_DIMID(nid,"time",timedimin)
1984         ierr=NF_INQ_VARID(nid,"time",timevar)
1985         if (ierr.NE.NF_NOERR) then
1986            write(*,*) 'ERROR: Field <time> is missing'
1987            stop "" 
1988         endif
1989         ierr=NF_INQ_DIMLEN(nid,timedimin,timelen)
1990
1991         ierr=NF_INQ_DIMID(nid,"lev",levdimin)
1992         ierr=NF_INQ_VARID(nid,"lev",levvar)
1993         if (ierr.NE.NF_NOERR) then
1994             write(*,*) 'ERROR: Field <lev> is lacking'
1995             stop ""
1996         endif
1997         ierr=NF_INQ_DIMLEN(nid,levdimin,levlen)
1998
1999         if((timelen/=ttm).or.(levlen/=llm)) then
2000            write(*,*) 'ERROR: Not the good lenght for axis'
2001            write(*,*) 'longitude: ',timelen,ttm+1
2002            write(*,*) 'latitude: ',levlen,llm
2003            stop "" 
2004         endif
2005
2006!#ifdef NC_DOUBLE
2007         ierr = NF_GET_VAR_DOUBLE(nid,timevar,time)
2008         ierr = NF_GET_VAR_DOUBLE(nid,levvar,lev)
2009!#else
2010!        ierr = NF_GET_VAR_REAL(nid,timevar,time)
2011!        ierr = NF_GET_VAR_REAL(nid,levvar,lev)
2012!#endif
2013
2014       return
2015       end
2016
2017!======================================================================
2018      SUBROUTINE read_rico(fich_rico,nlev_rico,ps_rico,play
2019     :             ,ts_rico,t_rico,q_rico,u_rico,v_rico,w_rico
2020     :             ,dth_dyn,dqh_dyn)
2021      implicit none
2022
2023c-------------------------------------------------------------------------
2024c Read RICO forcing data
2025c-------------------------------------------------------------------------
2026#include "dimensions.h"
2027
2028
2029      integer nlev_rico
2030      real ts_rico,ps_rico
2031      real t_rico(llm),q_rico(llm)
2032      real u_rico(llm),v_rico(llm)
2033      real w_rico(llm)
2034      real dth_dyn(llm)
2035      real dqh_dyn(llm)
2036     
2037
2038      real play(llm),zlay(llm)
2039     
2040
2041      real prico(nlev_rico),zrico(nlev_rico)
2042
2043      character*80 fich_rico
2044
2045      integer k,l
2046
2047     
2048      print*,fich_rico
2049      open(21,file=trim(fich_rico),form='formatted')
2050        do k=1,llm
2051      zlay(k)=0.
2052         enddo
2053     
2054        read(21,*) ps_rico,ts_rico
2055        prico(1)=ps_rico
2056        zrico(1)=0.0
2057      do l=2,nlev_rico
2058        read(21,*) k,prico(l),zrico(l)
2059      enddo
2060       close(21)
2061
2062      do k=1,llm
2063        do l=1,80
2064          if(prico(l)>play(k)) then
2065              if(play(k)>prico(l+1)) then
2066                zlay(k)=zrico(l)+(play(k)-prico(l)) *
2067     &              (zrico(l+1)-zrico(l))/(prico(l+1)-prico(l))
2068              else
2069                zlay(k)=zrico(l)+(play(k)-prico(80))*
2070     &              (zrico(81)-zrico(80))/(prico(81)-prico(80))
2071              endif
2072          endif
2073        enddo
2074        print*,k,zlay(k)
2075        ! U
2076        if(0 < zlay(k) .and. zlay(k) < 4000) then
2077          u_rico(k)=-9.9 + (-1.9 + 9.9)*zlay(k)/4000
2078        elseif(4000 < zlay(k) .and. zlay(k) < 12000) then
2079       u_rico(k)=  -1.9 + (30.0 + 1.9) /
2080     :          (12000 - 4000) * (zlay(k) - 4000)
2081        elseif(12000 < zlay(k) .and. zlay(k) < 13000) then
2082          u_rico(k)=30.0
2083        elseif(13000 < zlay(k) .and. zlay(k) < 20000) then
2084          u_rico(k)=30.0 - (30.0) /
2085     : (20000 - 13000) * (zlay(k) - 13000)
2086        else
2087          u_rico(k)=0.0
2088        endif
2089
2090!Q_v
2091        if(0 < zlay(k) .and. zlay(k) < 740) then
2092          q_rico(k)=16.0 + (13.8 - 16.0) / (740) * zlay(k)
2093        elseif(740 < zlay(k) .and. zlay(k) < 3260) then
2094          q_rico(k)=13.8 + (2.4 - 13.8) /
2095     :          (3260 - 740) * (zlay(k) - 740)
2096        elseif(3260 < zlay(k) .and. zlay(k) < 4000) then
2097          q_rico(k)=2.4 + (1.8 - 2.4) /
2098     :               (4000 - 3260) * (zlay(k) - 3260)
2099        elseif(4000 < zlay(k) .and. zlay(k) < 9000) then
2100          q_rico(k)=1.8 + (0 - 1.8) /
2101     :             (10000 - 4000) * (zlay(k) - 4000)
2102        else
2103          q_rico(k)=0.0
2104        endif
2105
2106!T
2107        if(0 < zlay(k) .and. zlay(k) < 740) then
2108          t_rico(k)=299.2 + (292.0 - 299.2) / (740) * zlay(k)
2109        elseif(740 < zlay(k) .and. zlay(k) < 4000) then
2110          t_rico(k)=292.0 + (278.0 - 292.0) /
2111     :       (4000 - 740) * (zlay(k) - 740)
2112        elseif(4000 < zlay(k) .and. zlay(k) < 15000) then
2113          t_rico(k)=278.0 + (203.0 - 278.0) /
2114     :       (15000 - 4000) * (zlay(k) - 4000)
2115        elseif(15000 < zlay(k) .and. zlay(k) < 17500) then
2116          t_rico(k)=203.0 + (194.0 - 203.0) /
2117     :       (17500 - 15000)* (zlay(k) - 15000)
2118        elseif(17500 < zlay(k) .and. zlay(k) < 20000) then
2119          t_rico(k)=194.0 + (206.0 - 194.0) /
2120     :       (20000 - 17500)* (zlay(k) - 17500)
2121        elseif(20000 < zlay(k) .and. zlay(k) < 60000) then
2122          t_rico(k)=206.0 + (270.0 - 206.0) /
2123     :        (60000 - 20000)* (zlay(k) - 20000)
2124        endif
2125
2126! W
2127        if(0 < zlay(k) .and. zlay(k) < 2260 ) then
2128          w_rico(k)=- (0.005/2260) * zlay(k)
2129        elseif(2260 < zlay(k) .and. zlay(k) < 4000 ) then
2130          w_rico(k)=- 0.005
2131        elseif(4000 < zlay(k) .and. zlay(k) < 5000 ) then
2132       w_rico(k)=- 0.005 + (0.005/ (5000 - 4000)) * (zlay(k) - 4000)
2133        else
2134          w_rico(k)=0.0
2135        endif
2136
2137! dThrz+dTsw0+dTlw0
2138        if(0 < zlay(k) .and. zlay(k) < 4000) then
2139          dth_dyn(k)=- 2.51 / 86400 + (-2.18 + 2.51 )/
2140     :               (86400*4000) * zlay(k)
2141        elseif(4000 < zlay(k) .and. zlay(k) < 5000) then
2142          dth_dyn(k)=- 2.18 / 86400 + ( 2.18 ) /
2143     :           (86400*(5000 - 4000)) * (zlay(k) - 4000)
2144        else
2145          dth_dyn(k)=0.0
2146        endif
2147! dQhrz
2148        if(0 < zlay(k) .and. zlay(k) < 3000) then
2149          dqh_dyn(k)=-1.0 / 86400 + (0.345 + 1.0)/
2150     :                    (86400*3000) * (zlay(k))
2151        elseif(3000 < zlay(k) .and. zlay(k) < 4000) then
2152          dqh_dyn(k)=0.345 / 86400
2153        elseif(4000 < zlay(k) .and. zlay(k) < 5000) then
2154          dqh_dyn(k)=0.345 / 86400 +
2155     +   (-0.345)/(86400 * (5000 - 4000)) * (zlay(k)-4000)
2156        else
2157          dqh_dyn(k)=0.0
2158        endif
2159
2160!?        if(play(k)>6e4) then
2161!?          ratqs0(1,k)=ratqsbas*(plev(1)-play(k))/(plev(1)-6e4)
2162!?        elseif((play(k)>3e4).and.(play(k)<6e4)) then
2163!?          ratqs0(1,k)=ratqsbas+(ratqshaut-ratqsbas)&
2164!?                          *(6e4-play(k))/(6e4-3e4)
2165!?        else
2166!?          ratqs0(1,k)=ratqshaut
2167!?        endif
2168
2169      enddo
2170
2171      do k=1,llm
2172      q_rico(k)=q_rico(k)/1e3
2173      dqh_dyn(k)=dqh_dyn(k)/1e3
2174      v_rico(k)=-3.8
2175      enddo
2176
2177          return
2178          end
2179
2180!=====================================================================
2181c-------------------------------------------------------------------------
2182      SUBROUTINE read_armcu(fich_armcu,nlev_armcu,nt_armcu,
2183     : sens,flat,adv_theta,rad_theta,adv_qt)
2184      implicit none
2185
2186c-------------------------------------------------------------------------
2187c Read ARM_CU case forcing data
2188c-------------------------------------------------------------------------
2189
2190      integer nlev_armcu,nt_armcu
2191      real sens(nt_armcu),flat(nt_armcu)
2192      real adv_theta(nt_armcu),rad_theta(nt_armcu),adv_qt(nt_armcu)
2193      character*80 fich_armcu
2194
2195      integer no,l,k,ip
2196      real riy,rim,rid,rih,bid
2197
2198      integer iy,im,id,ih,in
2199
2200      open(21,file=trim(fich_armcu),form='formatted')
2201      read(21,'(a)')
2202      do ip = 1, nt_armcu
2203      read(21,'(a)')
2204      read(21,'(a)')
2205      read(21,223) iy, im, id, ih, in, sens(ip),flat(ip),
2206     :             adv_theta(ip),rad_theta(ip),adv_qt(ip)
2207      print *,'forcages=',iy,im,id,ih,in, sens(ip),flat(ip),
2208     :             adv_theta(ip),rad_theta(ip),adv_qt(ip)
2209      enddo
2210      close(21)
2211
2212  223 format(5i3,5f8.3)
2213  226 format(f7.1,1x,10f8.2)
2214  227 format(f7.1,1x,1p,4e11.3)
2215  230 format(6f9.3,4e11.3)
2216
2217          return
2218          end
2219
2220!=====================================================================
2221       SUBROUTINE interp_toga_vertical(play,nlev_toga,plev_prof
2222     :         ,t_prof,q_prof,u_prof,v_prof,w_prof
2223     :         ,ht_prof,vt_prof,hq_prof,vq_prof
2224     :         ,t_mod,q_mod,u_mod,v_mod,w_mod
2225     :         ,ht_mod,vt_mod,hq_mod,vq_mod,mxcalc)
2226 
2227       implicit none
2228 
2229#include "dimensions.h"
2230
2231c-------------------------------------------------------------------------
2232c Vertical interpolation of TOGA-COARE forcing data onto model levels
2233c-------------------------------------------------------------------------
2234 
2235       integer nlevmax
2236       parameter (nlevmax=41)
2237       integer nlev_toga,mxcalc
2238!       real play(llm), plev_prof(nlevmax)
2239!       real t_prof(nlevmax),q_prof(nlevmax)
2240!       real u_prof(nlevmax),v_prof(nlevmax), w_prof(nlevmax)
2241!       real ht_prof(nlevmax),vt_prof(nlevmax)
2242!       real hq_prof(nlevmax),vq_prof(nlevmax)
2243 
2244       real play(llm), plev_prof(nlev_toga)
2245       real t_prof(nlev_toga),q_prof(nlev_toga)
2246       real u_prof(nlev_toga),v_prof(nlev_toga), w_prof(nlev_toga)
2247       real ht_prof(nlev_toga),vt_prof(nlev_toga)
2248       real hq_prof(nlev_toga),vq_prof(nlev_toga)
2249 
2250       real t_mod(llm),q_mod(llm)
2251       real u_mod(llm),v_mod(llm), w_mod(llm)
2252       real ht_mod(llm),vt_mod(llm)
2253       real hq_mod(llm),vq_mod(llm)
2254 
2255       integer l,k,k1,k2,kp
2256       real aa,frac,frac1,frac2,fact
2257 
2258       do l = 1, llm
2259
2260        if (play(l).ge.plev_prof(nlev_toga)) then
2261 
2262        mxcalc=l
2263         k1=0
2264         k2=0
2265
2266         if (play(l).le.plev_prof(1)) then
2267
2268         do k = 1, nlev_toga-1
2269          if (play(l).le.plev_prof(k)
2270     :       .and. play(l).gt.plev_prof(k+1)) then
2271            k1=k
2272            k2=k+1
2273          endif
2274         enddo
2275
2276         if (k1.eq.0 .or. k2.eq.0) then
2277          write(*,*) 'PB! k1, k2 = ',k1,k2
2278          write(*,*) 'l,play(l) = ',l,play(l)/100
2279         do k = 1, nlev_toga-1
2280          write(*,*) 'k,plev_prof(k) = ',k,plev_prof(k)/100
2281         enddo
2282         endif
2283
2284         frac = (plev_prof(k2)-play(l))/(plev_prof(k2)-plev_prof(k1))
2285         t_mod(l)= t_prof(k2) - frac*(t_prof(k2)-t_prof(k1))
2286         q_mod(l)= q_prof(k2) - frac*(q_prof(k2)-q_prof(k1))
2287         u_mod(l)= u_prof(k2) - frac*(u_prof(k2)-u_prof(k1))
2288         v_mod(l)= v_prof(k2) - frac*(v_prof(k2)-v_prof(k1))
2289         w_mod(l)= w_prof(k2) - frac*(w_prof(k2)-w_prof(k1))
2290         ht_mod(l)= ht_prof(k2) - frac*(ht_prof(k2)-ht_prof(k1))
2291         vt_mod(l)= vt_prof(k2) - frac*(vt_prof(k2)-vt_prof(k1))
2292         hq_mod(l)= hq_prof(k2) - frac*(hq_prof(k2)-hq_prof(k1))
2293         vq_mod(l)= vq_prof(k2) - frac*(vq_prof(k2)-vq_prof(k1))
2294     
2295         else !play>plev_prof(1)
2296
2297         k1=1
2298         k2=2
2299         frac1 = (play(l)-plev_prof(k2))/(plev_prof(k1)-plev_prof(k2))
2300         frac2 = (play(l)-plev_prof(k1))/(plev_prof(k1)-plev_prof(k2))
2301         t_mod(l)= frac1*t_prof(k1) - frac2*t_prof(k2)
2302         q_mod(l)= frac1*q_prof(k1) - frac2*q_prof(k2)
2303         u_mod(l)= frac1*u_prof(k1) - frac2*u_prof(k2)
2304         v_mod(l)= frac1*v_prof(k1) - frac2*v_prof(k2)
2305         w_mod(l)= frac1*w_prof(k1) - frac2*w_prof(k2)
2306         ht_mod(l)= frac1*ht_prof(k1) - frac2*ht_prof(k2)
2307         vt_mod(l)= frac1*vt_prof(k1) - frac2*vt_prof(k2)
2308         hq_mod(l)= frac1*hq_prof(k1) - frac2*hq_prof(k2)
2309         vq_mod(l)= frac1*vq_prof(k1) - frac2*vq_prof(k2)
2310
2311         endif ! play.le.plev_prof(1)
2312
2313        else ! above max altitude of forcing file
2314 
2315cjyg
2316         fact=20.*(plev_prof(nlev_toga)-play(l))/plev_prof(nlev_toga) !jyg
2317         fact = max(fact,0.)                                           !jyg
2318         fact = exp(-fact)                                             !jyg
2319         t_mod(l)= t_prof(nlev_toga)                                   !jyg
2320         q_mod(l)= q_prof(nlev_toga)*fact                              !jyg
2321         u_mod(l)= u_prof(nlev_toga)*fact                              !jyg
2322         v_mod(l)= v_prof(nlev_toga)*fact                              !jyg
2323         w_mod(l)= 0.0                                                 !jyg
2324         ht_mod(l)= ht_prof(nlev_toga)                                 !jyg
2325         vt_mod(l)= vt_prof(nlev_toga)                                 !jyg
2326         hq_mod(l)= hq_prof(nlev_toga)*fact                            !jyg
2327         vq_mod(l)= vq_prof(nlev_toga)*fact                            !jyg
2328 
2329        endif ! play
2330 
2331       enddo ! l
2332
2333!       do l = 1,llm
2334!       print *,'t_mod(l),q_mod(l),ht_mod(l),hq_mod(l) ',
2335!     $        l,t_mod(l),q_mod(l),ht_mod(l),hq_mod(l)
2336!       enddo
2337 
2338          return
2339          end
2340 
2341!======================================================================
2342        SUBROUTINE interp_toga_time(day,day1,annee_ref
2343     i             ,year_ini_toga,day_ini_toga,nt_toga,dt_toga,nlev_toga
2344     i             ,ts_toga,plev_toga,t_toga,q_toga,u_toga,v_toga,w_toga
2345     i             ,ht_toga,vt_toga,hq_toga,vq_toga
2346     o             ,ts_prof,plev_prof,t_prof,q_prof,u_prof,v_prof,w_prof
2347     o             ,ht_prof,vt_prof,hq_prof,vq_prof)
2348        implicit none
2349
2350!---------------------------------------------------------------------------------------
2351! Time interpolation of a 2D field to the timestep corresponding to day
2352!
2353! day: current julian day (e.g. 717538.2)
2354! day1: first day of the simulation
2355! nt_toga: total nb of data in the forcing (e.g. 480 for TOGA-COARE)
2356! dt_toga: total time interval (in sec) between 2 forcing data (e.g. 6h for TOGA-COARE)
2357!---------------------------------------------------------------------------------------
2358
2359#include "compar1d.h"
2360
2361! inputs:
2362        integer annee_ref
2363        integer nt_toga,nlev_toga
2364        integer year_ini_toga
2365        real day, day1,day_ini_toga,dt_toga
2366        real ts_toga(nt_toga)
2367        real plev_toga(nlev_toga,nt_toga),t_toga(nlev_toga,nt_toga)
2368        real q_toga(nlev_toga,nt_toga),u_toga(nlev_toga,nt_toga)
2369        real v_toga(nlev_toga,nt_toga),w_toga(nlev_toga,nt_toga)
2370        real ht_toga(nlev_toga,nt_toga),vt_toga(nlev_toga,nt_toga)
2371        real hq_toga(nlev_toga,nt_toga),vq_toga(nlev_toga,nt_toga)
2372! outputs:
2373        real ts_prof
2374        real plev_prof(nlev_toga),t_prof(nlev_toga)
2375        real q_prof(nlev_toga),u_prof(nlev_toga)
2376        real v_prof(nlev_toga),w_prof(nlev_toga)
2377        real ht_prof(nlev_toga),vt_prof(nlev_toga)
2378        real hq_prof(nlev_toga),vq_prof(nlev_toga)
2379! local:
2380        integer it_toga1, it_toga2,k
2381        real timeit,time_toga1,time_toga2,frac
2382
2383
2384        if (forcing_type.eq.2) then
2385! Check that initial day of the simulation consistent with TOGA-COARE period:
2386       if (annee_ref.ne.1992 .and. annee_ref.ne.1993) then
2387        print*,'Pour TOGA-COARE, annee_ref doit etre 1992 ou 1993'
2388        print*,'Changer annee_ref dans run.def'
2389        stop
2390       endif
2391       if (annee_ref.eq.1992 .and. day1.lt.day_ini_toga) then
2392        print*,'TOGA-COARE a débuté le 1er Nov 1992 (jour julien=306)'
2393        print*,'Changer dayref dans run.def'
2394        stop
2395       endif
2396       if (annee_ref.eq.1993 .and. day1.gt.day_ini_toga+119) then
2397        print*,'TOGA-COARE a fini le 28 Feb 1993 (jour julien=59)'
2398        print*,'Changer dayref ou nday dans run.def'
2399        stop
2400       endif
2401
2402       else if (forcing_type.eq.4) then
2403
2404! Check that initial day of the simulation consistent with TWP-ICE period:
2405       if (annee_ref.ne.2006) then
2406        print*,'Pour TWP-ICE, annee_ref doit etre 2006'
2407        print*,'Changer annee_ref dans run.def'
2408        stop
2409       endif
2410       if (annee_ref.eq.2006 .and. day1.lt.day_ini_toga) then
2411        print*,'TWP-ICE a debute le 17 Jan 2006 (jour julien=17)'
2412        print*,'Changer dayref dans run.def'
2413        stop
2414       endif
2415       if (annee_ref.eq.2006 .and. day1.gt.day_ini_toga+26) then
2416        print*,'TWP-ICE a fini le 12 Feb 2006 (jour julien=43)'
2417        print*,'Changer dayref ou nday dans run.def'
2418        stop
2419       endif
2420
2421       endif
2422
2423! Determine timestep relative to the 1st day of TOGA-COARE:
2424!       timeit=(day-day1)*86400.
2425!       if (annee_ref.eq.1992) then
2426!        timeit=(day-day_ini_toga)*86400.
2427!       else
2428!        timeit=(day+61.-1.)*86400. ! 61 days between Nov01 and Dec31 1992
2429!       endif
2430      timeit=(day-day_ini_toga)*86400
2431
2432! Determine the closest observation times:
2433       it_toga1=INT(timeit/dt_toga)+1
2434       it_toga2=it_toga1 + 1
2435       time_toga1=(it_toga1-1)*dt_toga
2436       time_toga2=(it_toga2-1)*dt_toga
2437
2438       if (it_toga1 .ge. nt_toga) then
2439        write(*,*) 'PB-stop: day, it_toga1, it_toga2, timeit: '
2440     :        ,day,it_toga1,it_toga2,timeit/86400.
2441        stop
2442       endif
2443
2444! time interpolation:
2445       frac=(time_toga2-timeit)/(time_toga2-time_toga1)
2446       frac=max(frac,0.0)
2447
2448       ts_prof = ts_toga(it_toga2)
2449     :          -frac*(ts_toga(it_toga2)-ts_toga(it_toga1))
2450
2451!        print*,
2452!     :'day,annee_ref,day_ini_toga,timeit,it_toga1,it_toga2,SST:',
2453!     :day,annee_ref,day_ini_toga,timeit/86400.,it_toga1,it_toga2,ts_prof
2454
2455       do k=1,nlev_toga
2456        plev_prof(k) = 100.*(plev_toga(k,it_toga2)
2457     :          -frac*(plev_toga(k,it_toga2)-plev_toga(k,it_toga1)))
2458        t_prof(k) = t_toga(k,it_toga2)
2459     :          -frac*(t_toga(k,it_toga2)-t_toga(k,it_toga1))
2460        q_prof(k) = q_toga(k,it_toga2)
2461     :          -frac*(q_toga(k,it_toga2)-q_toga(k,it_toga1))
2462        u_prof(k) = u_toga(k,it_toga2)
2463     :          -frac*(u_toga(k,it_toga2)-u_toga(k,it_toga1))
2464        v_prof(k) = v_toga(k,it_toga2)
2465     :          -frac*(v_toga(k,it_toga2)-v_toga(k,it_toga1))
2466        w_prof(k) = w_toga(k,it_toga2)
2467     :          -frac*(w_toga(k,it_toga2)-w_toga(k,it_toga1))
2468        ht_prof(k) = ht_toga(k,it_toga2)
2469     :          -frac*(ht_toga(k,it_toga2)-ht_toga(k,it_toga1))
2470        vt_prof(k) = vt_toga(k,it_toga2)
2471     :          -frac*(vt_toga(k,it_toga2)-vt_toga(k,it_toga1))
2472        hq_prof(k) = hq_toga(k,it_toga2)
2473     :          -frac*(hq_toga(k,it_toga2)-hq_toga(k,it_toga1))
2474        vq_prof(k) = vq_toga(k,it_toga2)
2475     :          -frac*(vq_toga(k,it_toga2)-vq_toga(k,it_toga1))
2476        enddo
2477
2478        return
2479        END
2480
2481!======================================================================
2482        SUBROUTINE interp_armcu_time(day,day1,annee_ref
2483     i             ,year_ini_armcu,day_ini_armcu,nt_armcu,dt_armcu
2484     i             ,nlev_armcu,fs_armcu,fl_armcu,at_armcu,rt_armcu
2485     i             ,aqt_armcu,fs_prof,fl_prof,at_prof,rt_prof,aqt_prof)
2486        implicit none
2487
2488!---------------------------------------------------------------------------------------
2489! Time interpolation of a 2D field to the timestep corresponding to day
2490!
2491! day: current julian day (e.g. 717538.2)
2492! day1: first day of the simulation
2493! nt_armcu: total nb of data in the forcing (e.g. 31 for armcu)
2494! dt_armcu: total time interval (in sec) between 2 forcing data (e.g. 1/2h for armcu)
2495! fs= sensible flux
2496! fl= latent flux
2497! at,rt,aqt= advective and radiative tendencies
2498!---------------------------------------------------------------------------------------
2499
2500! inputs:
2501        integer annee_ref
2502        integer nt_armcu,nlev_armcu
2503        integer year_ini_armcu
2504        real day, day1,day_ini_armcu,dt_armcu
2505        real fs_armcu(nt_armcu),fl_armcu(nt_armcu),at_armcu(nt_armcu)
2506        real rt_armcu(nt_armcu),aqt_armcu(nt_armcu)
2507! outputs:
2508        real fs_prof,fl_prof,at_prof,rt_prof,aqt_prof
2509! local:
2510        integer it_armcu1, it_armcu2,k
2511        real timeit,time_armcu1,time_armcu2,frac
2512
2513! Check that initial day of the simulation consistent with ARMCU period:
2514       if (annee_ref.ne.1997 ) then
2515        print*,'Pour ARMCU, annee_ref doit etre 1997 '
2516        print*,'Changer annee_ref dans run.def'
2517        stop
2518       endif
2519
2520      timeit=(day-day_ini_armcu)*86400
2521
2522! Determine the closest observation times:
2523       it_armcu1=INT(timeit/dt_armcu)+1
2524       it_armcu2=it_armcu1 + 1
2525       time_armcu1=(it_armcu1-1)*dt_armcu
2526       time_armcu2=(it_armcu2-1)*dt_armcu
2527       print *,'timeit day day_ini_armcu',timeit,day,day_ini_armcu
2528       print *,'it_armcu1,it_armcu2,time_armcu1,time_armcu2',
2529     .          it_armcu1,it_armcu2,time_armcu1,time_armcu2
2530
2531       if (it_armcu1 .ge. nt_armcu) then
2532        write(*,*) 'PB-stop: day, it_armcu1, it_armcu2, timeit: '
2533     :        ,day,it_armcu1,it_armcu2,timeit/86400.
2534        stop
2535       endif
2536
2537! time interpolation:
2538       frac=(time_armcu2-timeit)/(time_armcu2-time_armcu1)
2539       frac=max(frac,0.0)
2540
2541       fs_prof = fs_armcu(it_armcu2)
2542     :          -frac*(fs_armcu(it_armcu2)-fs_armcu(it_armcu1))
2543       fl_prof = fl_armcu(it_armcu2)
2544     :          -frac*(fl_armcu(it_armcu2)-fl_armcu(it_armcu1))
2545       at_prof = at_armcu(it_armcu2)
2546     :          -frac*(at_armcu(it_armcu2)-at_armcu(it_armcu1))
2547       rt_prof = rt_armcu(it_armcu2)
2548     :          -frac*(rt_armcu(it_armcu2)-rt_armcu(it_armcu1))
2549       aqt_prof = aqt_armcu(it_armcu2)
2550     :          -frac*(aqt_armcu(it_armcu2)-aqt_armcu(it_armcu1))
2551
2552         print*,
2553     :'day,annee_ref,day_ini_armcu,timeit,it_armcu1,it_armcu2,SST:',
2554     :day,annee_ref,day_ini_armcu,timeit/86400.,it_armcu1,
2555     :it_armcu2,fs_prof,fl_prof,at_prof,rt_prof,aqt_prof
2556
2557        return
2558        END
2559
2560!=====================================================================
2561      subroutine readprofiles(nlev_max,kmax,height,
2562     .           thlprof,qtprof,uprof,
2563     .           vprof,e12prof,ugprof,vgprof,
2564     .           wfls,dqtdxls,dqtdyls,dqtdtls,
2565     .           thlpcar)
2566      implicit none
2567
2568        integer nlev_max,kmax,kmax2
2569        logical :: llesread = .true.
2570
2571        real height(nlev_max),thlprof(nlev_max),qtprof(nlev_max),
2572     .       uprof(nlev_max),vprof(nlev_max),e12prof(nlev_max),
2573     .       ugprof(nlev_max),vgprof(nlev_max),wfls(nlev_max),
2574     .       dqtdxls(nlev_max),dqtdyls(nlev_max),dqtdtls(nlev_max),
2575     .           thlpcar(nlev_max)
2576
2577        integer, parameter :: ilesfile=1
2578        integer :: ierr,irad,imax,jtot,k
2579        logical :: lmoist,lcoriol,ltimedep
2580        real :: xsize,ysize
2581        real :: ustin,wsvsurf,timerad
2582        character(80) :: chmess
2583
2584        if(.not.(llesread)) return
2585
2586       open (ilesfile,file='prof.inp.001',status='old',iostat=ierr)
2587        if (ierr /= 0) stop 'ERROR:Prof.inp does not exist'
2588        read (ilesfile,*) kmax
2589        do k=1,kmax
2590          read (ilesfile,*) height(k),thlprof(k),qtprof (k),
2591     .                      uprof (k),vprof  (k),e12prof(k)
2592        enddo
2593        close(ilesfile)
2594
2595       open(ilesfile,file='lscale.inp.001',status='old',iostat=ierr)
2596        if (ierr /= 0) stop 'ERROR:Lscale.inp does not exist'
2597        read (ilesfile,*) kmax2
2598        if (kmax .ne. kmax2) then
2599          print *, 'fichiers prof.inp et lscale.inp incompatibles :'
2600          print *, 'nbre de niveaux : ',kmax,' et ',kmax2
2601          stop 'lecture profiles'
2602        endif
2603        do k=1,kmax
2604          read (ilesfile,*) height(k),ugprof(k),vgprof(k),wfls(k),
2605     .                      dqtdxls(k),dqtdyls(k),dqtdtls(k),thlpcar(k)
2606        end do
2607        close(ilesfile)
2608
2609        return
2610        end
2611
2612
2613!======================================================================
2614      subroutine readprofile_armcu(nlev_max,kmax,height,pprof,uprof,
2615     .       vprof,thetaprof,tprof,qvprof,rvprof,aprof,bprof)
2616!======================================================================
2617      implicit none
2618
2619        integer nlev_max,kmax,kmax2
2620        logical :: llesread = .true.
2621
2622        real height(nlev_max),pprof(nlev_max),tprof(nlev_max),
2623     .  thetaprof(nlev_max),rvprof(nlev_max),
2624     .  qvprof(nlev_max),uprof(nlev_max),vprof(nlev_max),
2625     .  aprof(nlev_max+1),bprof(nlev_max+1)
2626
2627        integer, parameter :: ilesfile=1
2628        integer, parameter :: ifile=2
2629        integer :: ierr,irad,imax,jtot,k
2630        logical :: lmoist,lcoriol,ltimedep
2631        real :: xsize,ysize
2632        real :: ustin,wsvsurf,timerad
2633        character(80) :: chmess
2634
2635        if(.not.(llesread)) return
2636
2637! Read profiles at full levels
2638       IF(nlev_max.EQ.19) THEN
2639       open (ilesfile,file='prof.inp.19',status='old',iostat=ierr)
2640       print *,'On ouvre prof.inp.19'
2641       ELSE
2642       open (ilesfile,file='prof.inp.40',status='old',iostat=ierr)
2643       print *,'On ouvre prof.inp.40'
2644       ENDIF
2645        if (ierr /= 0) stop 'ERROR:Prof.inp does not exist'
2646        read (ilesfile,*) kmax
2647        do k=1,kmax
2648          read (ilesfile,*) height(k)    ,pprof(k),  uprof(k), vprof(k),
2649     .                      thetaprof(k) ,tprof(k), qvprof(k),rvprof(k)
2650        enddo
2651        close(ilesfile)
2652
2653! Vertical coordinates half levels for eta-coordinates (plev = alpha + beta * psurf)
2654       IF(nlev_max.EQ.19) THEN
2655       open (ifile,file='proh.inp.19',status='old',iostat=ierr)
2656       print *,'On ouvre proh.inp.19'
2657       if (ierr /= 0) stop 'ERROR:Proh.inp.19 does not exist'
2658       ELSE
2659       open (ifile,file='proh.inp.40',status='old',iostat=ierr)
2660       print *,'On ouvre proh.inp.40'
2661       if (ierr /= 0) stop 'ERROR:Proh.inp.40 does not exist'
2662       ENDIF
2663        read (ifile,*) kmax
2664        do k=1,kmax
2665          read (ifile,*) jtot,aprof(k),bprof(k)
2666        enddo
2667        close(ifile)
2668
2669        return
2670        end
2671!===============================================================
2672      function ismin(n,sx,incx)
2673
2674      implicit none
2675      integer n,i,incx,ismin,ix
2676      real sx((n-1)*incx+1),sxmin
2677
2678      ix=1
2679      ismin=1
2680      sxmin=sx(1)
2681      do i=1,n-1
2682         ix=ix+incx
2683         if(sx(ix).lt.sxmin) then
2684             sxmin=sx(ix)
2685             ismin=i+1
2686         endif
2687      enddo
2688
2689      return
2690      end
2691
2692!===============================================================
2693      function ismax(n,sx,incx)
2694
2695      implicit none
2696      integer n,i,incx,ismax,ix
2697      real sx((n-1)*incx+1),sxmax
2698
2699      ix=1
2700      ismax=1
2701      sxmax=sx(1)
2702      do i=1,n-1
2703         ix=ix+incx
2704         if(sx(ix).gt.sxmax) then
2705             sxmax=sx(ix)
2706             ismax=i+1
2707         endif
2708      enddo
2709
2710      return
2711      end
2712
Note: See TracBrowser for help on using the repository browser.