source: LMDZ5/trunk/libf/phy1d/1DUTILS.h_with_writelim @ 1749

Last change on this file since 1749 was 1702, checked in by Laurent Fairhead, 12 years ago

To avoid conflict (incompatible and redundant open statements) with
conf_gcm FH

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