source: LMDZ6/trunk/libf/phylmd/dyn1d/m_1dconv_mod_h.f90 @ 5450

Last change on this file since 5450 was 5368, checked in by abarral, 4 weeks ago

(WIP) Turn implicit into explicit declarations
Turn 1dconv.h into a module to that end

  • Property svn:executable set to *
  • Property svn:keywords set to Id
File size: 30.4 KB
Line 
1MODULE m_1dconv_mod_h
2  IMPLICIT NONE; PRIVATE
3  PUBLIC get_uvd, copie, get_uvd2
4
5  REAL play(100)  !pression en Pa au milieu de chaque couche GCM
6  INTEGER JM(100) !pression en Pa au milieu de chaque couche GCM
7  REAL coef1(100) !coefficient d interpolation
8  REAL coef2(100) !coefficient d interpolation
9  INTEGER klev
10
11  INTEGER nblvlm !nombre de niveau de pression du mesoNH
12  REAL playm(100)  !pression en Pa au milieu de chaque couche Meso-NH
13  REAL hplaym(100) !pression en hPa milieux des couches Meso-NH
14
15CONTAINS
16
17        subroutine get_uvd(itap,dtime,file_forctl,file_fordat,                  &
18     &       ht,hq,hw,hu,hv,hthturb,hqturb,                                     &
19     &       Ts,imp_fcg,ts_fcg,Tp_fcg,Turb_fcg)                                 
20!
21        USE yomcst_mod_h
22implicit none
23
24!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
25! cette routine permet d obtenir u_convg,v_convg,ht,hq et ainsi de
26! pouvoir calculer la convergence et le cisaillement dans la physiq
27!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
28      integer i,j,k,ll,in
29
30      CHARACTER*80 file_forctl,file_fordat
31
32!======================================================================
33! methode: on va chercher les donnees du mesoNH de meteo france, on y
34!          a acces a tout pas detemps grace a la routine rdgrads qui
35!          est une boucle lisant dans ces fichiers.
36!          Puis on interpole ces donnes sur les 11 niveaux du gcm et
37!          et sur les pas de temps de ce meme gcm
38!----------------------------------------------------------------------
39! input:
40!       pasmax     :nombre de pas de temps maximum du mesoNH
41!       dt         :pas de temps du meso_NH (en secondes)
42!----------------------------------------------------------------------
43      integer pasmax,dt
44      save pasmax,dt
45!----------------------------------------------------------------------
46! arguments:
47!           itap   :compteur de la physique(le nombre de ces pas est
48!                   fixe dans la subroutine calcul_ini_gcm de interpo
49!                   -lation
50!           dtime  :pas detemps du gcm (en secondes)
51!           ht     :convergence horizontale de temperature(K/s)
52!           hq     :    "         "       d humidite (kg/kg/s)
53!           hw     :vitesse verticale moyenne (m/s**2)
54!           hu     :convergence horizontale d impulsion le long de x
55!                  (kg/(m^2 s^2)
56!           hv     : idem le long de y.
57!           Ts     : Temperature de surface (K)
58!           imp_fcg: var. logical .eq. T si forcage en impulsion
59!           ts_fcg: var. logical .eq. T si forcage en Ts present dans fichier
60!           Tp_fcg: var. logical .eq. T si forcage donne en Temp potentielle
61!           Turb_fcg: var. logical .eq. T si forcage turbulent present dans fichier
62!----------------------------------------------------------------------
63        integer itap
64        real dtime
65        real ht(:)
66        real hq(:)
67        real hu(:)
68        real hv(:)
69        real hw(:)
70        real hthturb(:)
71        real hqturb(:)
72        real Ts, Ts_subr
73        logical imp_fcg
74        logical ts_fcg
75        logical Tp_fcg
76        logical Turb_fcg
77!----------------------------------------------------------------------
78! Variables internes de get_uvd (note : l interpolation temporelle
79! est faite entre les pas de temps before et after, sur les variables
80! definies sur la grille du SCM; on atteint exactement les valeurs Meso
81! aux milieux des pas de temps Meso)
82!     time0     :date initiale en secondes
83!     time      :temps associe a chaque pas du SCM
84!     pas       :numero du pas du meso_NH (on lit en pas : le premier pas
85!                 des donnees est duplique)
86!     pasprev   :numero du pas de lecture precedent
87!     htaft     :advection horizontale de temp. au pas de temps after
88!     hqaft     :    "         "      d humidite        "
89!     hwaft     :vitesse verticalle moyenne  au pas de temps after
90!     huaft,hvaft :advection horizontale d impulsion au pas de temps after
91!     tsaft     : surface temperature 'after time step'
92!     htbef     :idem htaft, mais pour le pas de temps before
93!     hqbef     :voir hqaft
94!     hwbef     :voir hwaft
95!     hubef,hvbef : idem huaft,hvaft, mais pour before
96!     tsbef     : surface temperature 'before time step'
97!----------------------------------------------------------------------
98        integer time0,pas,pasprev
99        save time0,pas,pasprev
100        real time
101        real htaft(100),hqaft(100),hwaft(100),huaft(100),hvaft(100)
102        real hthturbaft(100),hqturbaft(100)
103        real Tsaft
104        save htaft,hqaft,hwaft,huaft,hvaft,hthturbaft,hqturbaft
105        real htbef(100),hqbef(100),hwbef(100),hubef(100),hvbef(100)
106        real hthturbbef(100),hqturbbef(100)
107        real Tsbef
108        save htbef,hqbef,hwbef,hubef,hvbef,hthturbbef,hqturbbef
109!
110        real timeaft,timebef
111        save timeaft,timebef
112        integer temps
113        character*4 string
114!----------------------------------------------------------------------
115! variables arguments de la subroutine rdgrads
116!---------------------------------------------------------------------
117        integer icompt,icomp1 !compteurs de rdgrads
118        real z(100)         ! altitude (grille Meso)
119        real ht_mes(100)    !convergence horizontale de temperature
120                            !-(grille Meso)
121        real hq_mes(100)    !convergence horizontale d humidite
122                            !(grille Meso)
123        real hw_mes(100)    !vitesse verticale moyenne
124                            !(grille Meso)
125        real hu_mes(100),hv_mes(100)    !convergence horizontale d impulsion
126                                        !(grille Meso)
127        real hthturb_mes(100) !tendance horizontale de T_pot, due aux
128                              !flux turbulents
129        real hqturb_mes(100) !tendance horizontale d humidite, due aux
130                              !flux turbulents
131!
132!---------------------------------------------------------------------
133! variable argument de la subroutine copie
134!---------------------------------------------------------------------
135! SB        real pplay(100)    !pression en milieu de couche du gcm
136! SB                            !argument de la physique
137!---------------------------------------------------------------------
138! variables destinees a la lecture du pas de temps du fichier de donnees
139!---------------------------------------------------------------------
140       character*80 aaa,atemps,apasmax
141       integer nch,imn,ipa
142!---------------------------------------------------------------------
143               print*,'le pas itap est:',itap
144!*** on determine le pas du meso_NH correspondant au nouvel itap ***
145!*** pour aller chercher les champs dans rdgrads                 ***
146!
147        time=time0+itap*dtime
148!c        temps=int(time/dt+1)
149!c        pas=min(temps,pasmax)
150        temps = 1 + int((dt + 2*time)/(2*dt))
151        pas=min(temps,pasmax-1)
152             print*,'le pas Meso est:',pas
153!
154!
155!===================================================================
156!
157!*** on remplit les champs before avec les champs after du pas   ***
158!*** precedent en format gcm                                     ***
159        if(pas.gt.pasprev)then
160          do i=1,klev
161             htbef(i)=htaft(i)
162             hqbef(i)=hqaft(i)
163             hwbef(i)=hwaft(i)
164             hubef(i)=huaft(i)
165             hvbef(i)=hvaft(i)
166             hThTurbbef(i)=hThTurbaft(i)
167             hqTurbbef(i)=hqTurbaft(i)
168          enddo
169          tsbef = tsaft
170          timebef=pasprev*dt
171          timeaft=timebef+dt
172          icomp1 = nblvlm*4
173          IF (ts_fcg) icomp1 = icomp1 + 1
174          IF (imp_fcg) icomp1 = icomp1 + nblvlm*2
175          IF (Turb_fcg) icomp1 = icomp1 + nblvlm*2
176          icompt = icomp1*pas
177         print *, 'imp_fcg,ts_fcg,Turb_fcg,pas,nblvlm,icompt'
178         print *, imp_fcg,ts_fcg,Turb_fcg,pas,nblvlm,icompt
179                       print*,'le pas pas est:',pas
180!*** on va chercher les nouveaux champs after dans toga.dat     ***
181!*** champs en format meso_NH                                   ***
182          open(99,FILE=file_fordat,FORM='UNFORMATTED',                        &
183     &             ACCESS='DIRECT',RECL=8)
184          call rdgrads(99,icompt,nblvlm,z,ht_mes,hq_mes,hw_mes                &
185     &                  ,hu_mes,hv_mes,hthturb_mes,hqturb_mes                 &
186     &                  ,ts_fcg,ts_subr,imp_fcg,Turb_fcg)
187!
188
189               if(Tp_fcg) then
190!     (le forcage est donne en temperature potentielle)
191         do i = 1,nblvlm
192           ht_mes(i) = ht_mes(i)*(hplaym(i)/1000.)**rkappa
193         enddo
194               endif ! Tp_fcg
195        if(Turb_fcg) then
196         do i = 1,nblvlm
197           hThTurb_mes(i) = hThTurb_mes(i)*(hplaym(i)/1000.)**rkappa
198         enddo
199        endif  ! Turb_fcg
200!
201               print*,'ht_mes ',(ht_mes(i),i=1,nblvlm)
202               print*,'hq_mes ',(hq_mes(i),i=1,nblvlm)
203               print*,'hw_mes ',(hw_mes(i),i=1,nblvlm)
204                  if(imp_fcg) then
205               print*,'hu_mes ',(hu_mes(i),i=1,nblvlm)
206               print*,'hv_mes ',(hv_mes(i),i=1,nblvlm)
207                  endif
208                  if(Turb_fcg) then
209               print*,'hThTurb_mes ',(hThTurb_mes(i),i=1,nblvlm)
210               print*,'hqTurb_mes ',(hqTurb_mes(i),i=1,nblvlm)
211                  endif
212          IF (ts_fcg) print*,'ts_subr', ts_subr
213!*** on interpole les champs meso_NH sur les niveaux de pression***
214!*** gcm . on obtient le nouveau champ after                    ***
215            do k=1,klev
216             if (JM(k) .eq. 0) then
217         htaft(k)=              ht_mes(jm(k)+1)
218         hqaft(k)=              hq_mes(jm(k)+1)
219         hwaft(k)=              hw_mes(jm(k)+1)
220               if(imp_fcg) then
221           huaft(k)=              hu_mes(jm(k)+1)
222           hvaft(k)=              hv_mes(jm(k)+1)
223               endif ! imp_fcg
224               if(Turb_fcg) then
225           hThTurbaft(k)=         hThTurb_mes(jm(k)+1)
226           hqTurbaft(k)=          hqTurb_mes(jm(k)+1)
227               endif ! Turb_fcg
228             else ! JM(k) .eq. 0
229           htaft(k)=coef1(k)*ht_mes(jm(k))+coef2(k)*ht_mes(jm(k)+1)
230           hqaft(k)=coef1(k)*hq_mes(jm(k))+coef2(k)*hq_mes(jm(k)+1)
231           hwaft(k)=coef1(k)*hw_mes(jm(k))+coef2(k)*hw_mes(jm(k)+1)
232               if(imp_fcg) then
233           huaft(k)=coef1(k)*hu_mes(jm(k))+coef2(k)*hu_mes(jm(k)+1)
234           hvaft(k)=coef1(k)*hv_mes(jm(k))+coef2(k)*hv_mes(jm(k)+1)
235               endif ! imp_fcg
236               if(Turb_fcg) then
237           hThTurbaft(k)=coef1(k)*hThTurb_mes(jm(k))                            &
238     &               +coef2(k)*hThTurb_mes(jm(k)+1)
239           hqTurbaft(k) =coef1(k)*hqTurb_mes(jm(k))                             &
240     &               +coef2(k)*hqTurb_mes(jm(k)+1)
241               endif ! Turb_fcg
242             endif ! JM(k) .eq. 0
243            enddo
244            tsaft = ts_subr
245            pasprev=pas
246         else ! pas.gt.pasprev
247            print*,'timebef est:',timebef
248         endif ! pas.gt.pasprev    fin du bloc relatif au passage au pas
249                                  !de temps (meso) suivant
250!*** si on atteint le pas max des donnees experimentales ,on     ***
251!*** on conserve les derniers champs calcules                    ***
252      if(temps.ge.pasmax)then
253          do ll=1,klev
254               ht(ll)=htaft(ll)
255               hq(ll)=hqaft(ll)
256               hw(ll)=hwaft(ll)
257               hu(ll)=huaft(ll)
258               hv(ll)=hvaft(ll)
259               hThTurb(ll)=hThTurbaft(ll)
260               hqTurb(ll)=hqTurbaft(ll)
261          enddo
262          ts_subr = tsaft
263      else ! temps.ge.pasmax
264!*** on interpole sur les pas de temps de 10mn du gcm a partir   ***
265!** des pas de temps de 1h du meso_NH                            ***
266         do j=1,klev
267         ht(j)=((timeaft-time)*htbef(j)+(time-timebef)*htaft(j))/dt
268         hq(j)=((timeaft-time)*hqbef(j)+(time-timebef)*hqaft(j))/dt
269         hw(j)=((timeaft-time)*hwbef(j)+(time-timebef)*hwaft(j))/dt
270             if(imp_fcg) then
271         hu(j)=((timeaft-time)*hubef(j)+(time-timebef)*huaft(j))/dt
272         hv(j)=((timeaft-time)*hvbef(j)+(time-timebef)*hvaft(j))/dt
273             endif ! imp_fcg
274             if(Turb_fcg) then
275         hThTurb(j)=((timeaft-time)*hThTurbbef(j)                           &
276     &           +(time-timebef)*hThTurbaft(j))/dt
277         hqTurb(j)= ((timeaft-time)*hqTurbbef(j)                            &
278     &           +(time-timebef)*hqTurbaft(j))/dt
279             endif ! Turb_fcg
280         enddo
281         ts_subr = ((timeaft-time)*tsbef + (time-timebef)*tsaft)/dt
282       endif ! temps.ge.pasmax
283!
284        print *,' time,timebef,timeaft',time,timebef,timeaft
285        print *,' ht,htbef,htaft,hthturb,hthturbbef,hthturbaft'
286        do j= 1,klev
287           print *, j,ht(j),htbef(j),htaft(j),                              &
288     &             hthturb(j),hthturbbef(j),hthturbaft(j)
289        enddo
290        print *,' hq,hqbef,hqaft,hqturb,hqturbbef,hqturbaft'
291        do j= 1,klev
292           print *, j,hq(j),hqbef(j),hqaft(j),                              &
293     &             hqturb(j),hqturbbef(j),hqturbaft(j)
294        enddo
295!
296!-------------------------------------------------------------------
297!
298         IF (Ts_fcg) Ts = Ts_subr
299         return
300!
301!-----------------------------------------------------------------------
302! on sort les champs de "convergence" pour l instant initial 'in'
303! ceci se passe au pas temps itap=0 de la physique
304!-----------------------------------------------------------------------
305        entry get_uvd2(itap,dtime,file_forctl,file_fordat,                  &
306     &           ht,hq,hw,hu,hv,hThTurb,hqTurb,ts,                          &
307     &           imp_fcg,ts_fcg,Tp_fcg,Turb_fcg)
308             print*,'le pas itap est:',itap
309!
310!===================================================================
311!
312       write(*,'(a)') 'OPEN '//file_forctl
313       open(97,FILE=file_forctl,FORM='FORMATTED')
314!
315!------------------
316      do i=1,1000
317      read(97,1000,end=999) string
318 1000 format (a4)
319      if (string .eq. 'TDEF') go to 50
320      enddo
321 50   backspace(97)
322!-------------------------------------------------------------------
323!   *** on lit le pas de temps dans le fichier de donnees ***
324!   *** "forcing.ctl" et pasmax                           ***
325!-------------------------------------------------------------------
326      read(97,2000) aaa
327 2000  format (a80)
328         print*,'aaa est',aaa
329      aaa=spaces(aaa,1)
330         print*,'aaa',aaa
331      call getsch(aaa,' ',' ',5,atemps,nch)
332         print*,'atemps est',atemps
333        atemps=atemps(1:nch-2)
334         print*,'atemps',atemps
335        read(atemps,*) imn
336        dt=imn*60
337         print*,'le pas de temps dt',dt
338      call getsch(aaa,' ',' ',2,apasmax,nch)
339        apasmax=apasmax(1:nch)
340        read(apasmax,*) ipa
341        pasmax=ipa
342         print*,'pasmax est',pasmax
343      CLOSE(97)
344!------------------------------------------------------------------
345! *** on lit le pas de temps initial de la simulation ***
346!------------------------------------------------------------------
347                  in=itap
348!c                  pasprev=in
349!c                  time0=dt*(pasprev-1)
350                  pasprev=in-1
351                  time0=dt*pasprev
352!
353          close(98)
354!
355      write(*,'(a)') 'OPEN '//file_fordat
356      open(99,FILE=file_fordat,FORM='UNFORMATTED',                          &
357     &          ACCESS='DIRECT',RECL=8)
358          icomp1 = nblvlm*4
359          IF (ts_fcg) icomp1 = icomp1 + 1
360          IF (imp_fcg) icomp1 = icomp1 + nblvlm*2
361          IF (Turb_fcg) icomp1 = icomp1 + nblvlm*2
362          icompt = icomp1*(in-1)
363          call rdgrads(99,icompt,nblvlm,z,ht_mes,hq_mes,hw_mes              &
364     &                   ,hu_mes,hv_mes,hthturb_mes,hqturb_mes              &
365     &                   ,ts_fcg,ts_subr,imp_fcg,Turb_fcg)
366          print *, 'get_uvd : rdgrads ->'
367          print *, tp_fcg
368!
369! following commented out because we have temperature already in ARM case
370!   (otherwise this is the potential temperature )
371!------------------------------------------------------------------------
372               if(Tp_fcg) then
373          do i = 1,nblvlm
374            ht_mes(i) = ht_mes(i)*(hplaym(i)/1000.)**rkappa
375          enddo
376               endif ! Tp_fcg
377           print*,'ht_mes ',(ht_mes(i),i=1,nblvlm)
378           print*,'hq_mes ',(hq_mes(i),i=1,nblvlm)
379           print*,'hw_mes ',(hw_mes(i),i=1,nblvlm)
380              if(imp_fcg) then
381           print*,'hu_mes ',(hu_mes(i),i=1,nblvlm)
382           print*,'hv_mes ',(hv_mes(i),i=1,nblvlm)
383           print*,'t',ts_subr
384              endif ! imp_fcg
385              if(Turb_fcg) then
386                 print*,'hThTurb_mes ',(hThTurb_mes(i),i=1,nblvlm)
387                 print*,'hqTurb ',     (hqTurb_mes(i),i=1,nblvlm)
388              endif ! Turb_fcg
389!----------------------------------------------------------------------
390! on a obtenu des champs initiaux sur les niveaux du meso_NH
391! on interpole sur les niveaux du gcm(niveau pression bien sur!)
392!-----------------------------------------------------------------------
393            do k=1,klev
394             if (JM(k) .eq. 0) then
395!FKC bug? ne faut il pas convertir tsol en tendance ????
396!RT bug taken care of by removing the stuff
397           htaft(k)=              ht_mes(jm(k)+1)
398           hqaft(k)=              hq_mes(jm(k)+1)
399           hwaft(k)=              hw_mes(jm(k)+1)
400               if(imp_fcg) then
401           huaft(k)=              hu_mes(jm(k)+1)
402           hvaft(k)=              hv_mes(jm(k)+1)
403               endif ! imp_fcg
404               if(Turb_fcg) then
405           hThTurbaft(k)=         hThTurb_mes(jm(k)+1)
406           hqTurbaft(k)=          hqTurb_mes(jm(k)+1)
407               endif ! Turb_fcg
408             else ! JM(k) .eq. 0
409           htaft(k)=coef1(k)*ht_mes(jm(k))+coef2(k)*ht_mes(jm(k)+1)
410           hqaft(k)=coef1(k)*hq_mes(jm(k))+coef2(k)*hq_mes(jm(k)+1)
411           hwaft(k)=coef1(k)*hw_mes(jm(k))+coef2(k)*hw_mes(jm(k)+1)
412               if(imp_fcg) then
413           huaft(k)=coef1(k)*hu_mes(jm(k))+coef2(k)*hu_mes(jm(k)+1)
414           hvaft(k)=coef1(k)*hv_mes(jm(k))+coef2(k)*hv_mes(jm(k)+1)
415               endif ! imp_fcg
416               if(Turb_fcg) then
417           hThTurbaft(k)=coef1(k)*hThTurb_mes(jm(k))                        &
418     &                  +coef2(k)*hThTurb_mes(jm(k)+1)
419           hqTurbaft(k) =coef1(k)*hqTurb_mes(jm(k))                         &
420     &                  +coef2(k)*hqTurb_mes(jm(k)+1)
421               endif ! Turb_fcg
422             endif ! JM(k) .eq. 0
423            enddo
424            tsaft = ts_subr
425! valeurs initiales des champs de convergence
426          do k=1,klev
427             ht(k)=htaft(k)
428             hq(k)=hqaft(k)
429             hw(k)=hwaft(k)
430                if(imp_fcg) then
431             hu(k)=huaft(k)
432             hv(k)=hvaft(k)
433                endif ! imp_fcg
434                if(Turb_fcg) then
435             hThTurb(k)=hThTurbaft(k)
436             hqTurb(k) =hqTurbaft(k)
437                endif ! Turb_fcg
438          enddo
439          ts_subr = tsaft
440          close(99)
441          close(98)
442!
443!-------------------------------------------------------------------
444!
445!
446 100      IF (Ts_fcg) Ts = Ts_subr
447        return
448!
449999     continue
450        stop 'erreur lecture, file forcing.ctl'
451        end
452
453      SUBROUTINE advect_tvl(dtime,zt,zq,vu_f,vv_f,t_f,q_f                   &
454     &                     ,d_t_adv,d_q_adv)
455      use dimphy
456      USE dimensions_mod, ONLY: iim, jjm, llm, ndm
457implicit none
458
459
460!cccc      INCLUDE "dimphy.h"
461
462      integer k
463      real dtime, fact, du, dv, cx, cy, alx, aly
464      real zt(klev), zq(klev,3)
465      real vu_f(klev), vv_f(klev), t_f(klev), q_f(klev,3)
466
467      real d_t_adv(klev), d_q_adv(klev,3)
468
469! Velocity of moving cell
470      data cx,cy /12., -2./
471
472! Dimensions of moving cell
473      data alx,aly /100000.,150000./
474
475      do k = 1, klev
476            du = abs(vu_f(k)-cx)/alx
477            dv = abs(vv_f(k)-cy)/aly
478            fact = dtime *(du+dv-du*dv*dtime)
479            d_t_adv(k) = fact * (t_f(k)-zt(k))
480            d_q_adv(k,1) = fact * (q_f(k,1)-zq(k,1))
481            d_q_adv(k,2) = fact * (q_f(k,2)-zq(k,2))
482            d_q_adv(k,3) = fact * (q_f(k,3)-zq(k,3))
483      enddo
484
485      return
486      end
487
488      SUBROUTINE copie(klevgcm,playgcm,psolgcm,file_forctl)
489      implicit none
490
491!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
492! cette routine remplit les variables du module
493!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
494
495      integer k,klevgcm
496      real playgcm(klevgcm) ! pression en milieu de couche du gcm
497      real psolgcm
498      character*80 file_forctl
499
500      klev = klevgcm
501
502!---------------------------------------------------------------------
503! pression au milieu des couches du gcm dans la physiq
504! (SB: remplace le call conv_lipress_gcm(playgcm) )
505!---------------------------------------------------------------------
506
507       do k = 1, klev
508        play(k) = playgcm(k)
509        print*,'la pression gcm est:',play(k)
510       enddo
511
512!----------------------------------------------------------------------
513! lecture du descripteur des donnees Meso-NH (forcing.ctl):
514!  -> nb niveaux du meso.NH (nblvlm) + pressions meso.NH
515! (on remplit le COMMON com2_phys_gcss)
516!----------------------------------------------------------------------
517
518      call mesolupbis(file_forctl)
519
520      print*,'la valeur de nblvlm est:',nblvlm
521
522!----------------------------------------------------------------------
523! etude de la correspondance entre les niveaux meso.NH et GCM;
524! calcul des coefficients d interpolation coef1 et coef2
525! (on remplit le COMMON com1_phys_gcss)
526!----------------------------------------------------------------------
527
528      call corresbis(psolgcm)
529
530!---------------------------------------------------------
531! TEST sur le remplissage de com1_phys_gcss et com2_phys_gcss:
532!---------------------------------------------------------
533 
534      write(*,*) ' '
535      write(*,*) 'TESTS com1_phys_gcss et com2_phys_gcss dans copie.F'
536      write(*,*) '--------------------------------------'
537      write(*,*) 'GCM: nb niveaux:',klev,' et pression, coeffs:'
538      do k = 1, klev
539      write(*,*) play(k), coef1(k), coef2(k)
540      enddo
541      write(*,*) 'MESO-NH: nb niveaux:',nblvlm,' et pression:'
542      do k = 1, nblvlm
543      write(*,*) playm(k), hplaym(k)
544      enddo
545      write(*,*) ' '
546
547      end
548      SUBROUTINE mesolupbis(file_forctl)
549      implicit none
550!
551!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
552!
553! Lecture descripteur des donnees MESO-NH (forcing.ctl):
554! -------------------------------------------------------
555!
556!     Cette subroutine lit dans le fichier de controle "essai.ctl"
557!     et affiche le nombre de niveaux du Meso-NH ainsi que les valeurs
558!     des pressions en milieu de couche du Meso-NH (en Pa puis en hPa).
559!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
560!
561
562      INTEGER i,lu,mlz,mlzh
563 
564      character*80 file_forctl
565
566      character*4 a
567      character*80 aaa,anblvl
568      integer nch
569
570      lu=9
571      open(lu,file=file_forctl,form='formatted')
572!
573      do i=1,1000
574      read(lu,1000,end=999) a
575      if (a .eq. 'ZDEF') go to 100
576      enddo
577!
578 100  backspace(lu)
579      print*,'  DESCRIPTION DES 2 MODELES : '
580      print*,' '
581!
582      read(lu,2000) aaa
583 2000  format (a80)
584       aaa=spaces(aaa,1)
585       call getsch(aaa,' ',' ',2,anblvl,nch)
586         read(anblvl,*) nblvlm
587
588!
589      print*,'nbre de niveaux de pression Meso-NH :',nblvlm
590      print*,' '
591      print*,'pression en Pa de chaque couche du meso-NH :'
592!
593      read(lu,*) (playm(mlz),mlz=1,nblvlm)
594!      Si la pression est en HPa, la multiplier par 100
595      if (playm(1) .lt. 10000.) then
596        do mlz = 1,nblvlm
597         playm(mlz) = playm(mlz)*100.
598        enddo
599      endif
600      print*,(playm(mlz),mlz=1,nblvlm)
601!
602 1000 format (a4)
603 1001 format(5x,i2)
604!
605      print*,' '
606      do mlzh=1,nblvlm
607      hplaym(mlzh)=playm(mlzh)/100.
608      enddo
609!
610      print*,'pression en hPa de chaque couche du meso-NH: '
611      print*,(hplaym(mlzh),mlzh=1,nblvlm)
612!
613      close (lu)
614      return
615!
616 999  stop 'erreur lecture des niveaux pression des donnees'
617      end
618
619      SUBROUTINE rdgrads(itape,icount,nl,z,ht,hq,hw,hu,hv,hthtur,hqtur,     &
620     &  ts_fcg,ts,imp_fcg,Turb_fcg)
621      IMPLICIT none
622      INTEGER itape,icount,icomp, nl
623      real z(nl),ht(nl),hq(nl),hw(nl),hu(nl),hv(nl)
624      real hthtur(nl),hqtur(nl)
625      real ts
626!
627      INTEGER k
628!
629      LOGICAL imp_fcg,ts_fcg,Turb_fcg
630!
631      icomp = icount
632!
633!
634         do k=1,nl
635            icomp=icomp+1
636            read(itape,rec=icomp)z(k)
637            print *,'icomp,k,z(k) ',icomp,k,z(k)
638         enddo
639         do k=1,nl
640            icomp=icomp+1
641            read(itape,rec=icomp)hT(k)
642             print*, hT(k), k
643         enddo
644         do k=1,nl
645            icomp=icomp+1
646            read(itape,rec=icomp)hQ(k)
647         enddo
648!
649             if(turb_fcg) then
650         do k=1,nl
651            icomp=icomp+1
652           read(itape,rec=icomp)hThTur(k)
653         enddo
654         do k=1,nl
655            icomp=icomp+1
656           read(itape,rec=icomp)hqTur(k)
657         enddo
658             endif
659         print *,' apres lecture hthtur, hqtur'
660!
661          if(imp_fcg) then
662
663         do k=1,nl
664            icomp=icomp+1
665           read(itape,rec=icomp)hu(k)
666         enddo
667         do k=1,nl
668            icomp=icomp+1
669            read(itape,rec=icomp)hv(k)
670         enddo
671
672          endif
673!
674         do k=1,nl
675            icomp=icomp+1
676            read(itape,rec=icomp)hw(k)
677         enddo
678!
679              if(ts_fcg) then
680         icomp=icomp+1
681         read(itape,rec=icomp)ts
682              endif
683!
684      print *,' rdgrads ->'
685
686      RETURN
687      END
688
689      SUBROUTINE corresbis(psol)
690      implicit none
691
692!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
693! Cette subroutine calcule et affiche les valeurs des coefficients
694! d interpolation qui serviront dans la formule d interpolation elle-
695! meme.
696!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
697
698      REAL psol
699      REAL val
700      INTEGER k, mlz
701
702
703      do k=1,klev
704         val=play(k)
705       if (val .gt. playm(1)) then
706          mlz = 0
707          JM(1) = mlz
708          coef1(1)=(playm(mlz+1)-val)/(playm(mlz+1)-psol)
709          coef2(1)=(val-psol)/(playm(mlz+1)-psol)
710       else if (val .gt. playm(nblvlm)) then
711         do mlz=1,nblvlm
712          if (     val .le. playm(mlz).and. val .gt. playm(mlz+1))then
713           JM(k)=mlz
714           coef1(k)=(playm(mlz+1)-val)/(playm(mlz+1)-playm(mlz))
715           coef2(k)=(val-playm(mlz))/(playm(mlz+1)-playm(mlz))
716          endif
717         enddo
718       else
719         JM(k) = nblvlm-1
720         coef1(k) = 0.
721         coef2(k) = 0.
722       endif
723      enddo
724!
725!c      if (play(klev) .le. playm(nblvlm)) then
726!c         mlz=nblvlm-1
727!c         JM(klev)=mlz
728!c         coef1(klev)=(playm(mlz+1)-val)
729!c     *            /(playm(mlz+1)-playm(mlz))
730!c         coef2(klev)=(val-playm(mlz))
731!c     *            /(playm(mlz+1)-playm(mlz))
732!c      endif
733!
734      print*,' '
735      print*,'         INTERPOLATION  : '
736      print*,' '
737      print*,'correspondance de 9 niveaux du GCM sur les 53 du meso-NH:'
738      print*,(JM(k),k=1,klev)
739      print*,'correspondance de 9 niveaux du GCM sur les 53 du meso-NH:'
740      print*,(JM(k),k=1,klev)
741      print*,' '
742      print*,'vals du premier coef d"interpolation pour les 9 niveaux: '
743      print*,(coef1(k),k=1,klev)
744      print*,' '
745      print*,'valeurs du deuxieme coef d"interpolation pour les 9 niveaux:'
746      print*,(coef2(k),k=1,klev)
747!
748      return
749      end
750      SUBROUTINE GETSCH(STR,DEL,TRM,NTH,SST,NCH)
751!***************************************************************
752!*                                                             *
753!*                                                             *
754!* GETSCH                                                      *
755!*                                                             *
756!*                                                             *
757!* modified by :                                               *
758!***************************************************************
759!*   Return in SST the character string found between the NTH-1 and NTH
760!*   occurence of the delimiter 'DEL' but before the terminator 'TRM' in
761!*   the input string 'STR'. If TRM=DEL then STR is considered unlimited.
762!*   NCH=Length of the string returned in SST or =-1 if NTH is <1 or if
763!*   NTH is greater than the number of delimiters in STR.
764      IMPLICIT INTEGER (A-Z)
765      CHARACTER STR*(*),DEL*1,TRM*1,SST*(*)
766      NCH=-1
767      SST=' '
768      IF(NTH.GT.0) THEN
769        IF(TRM.EQ.DEL) THEN
770          LENGTH=LEN(STR)
771        ELSE
772          LENGTH=INDEX(STR,TRM)-1
773          IF(LENGTH.LT.0) LENGTH=LEN(STR)
774        ENDIF
775!*     Find beginning and end of the NTH DEL-limited substring in STR
776        END=-1
777        DO 1,N=1,NTH
778        IF(END.EQ.LENGTH) RETURN
779        BEG=END+2
780        END=BEG+INDEX(STR(BEG:LENGTH),DEL)-2
781        IF(END.EQ.BEG-2) END=LENGTH
782!*        PRINT *,'NTH,LENGTH,N,BEG,END=',NTH,LENGTH,N,BEG,END
783    1   CONTINUE
784        NCH=END-BEG+1
785        IF(NCH.GT.0) SST=STR(BEG:END)
786      ENDIF
787      END
788
789      CHARACTER*(80) FUNCTION SPACES(STR,NSPACE)
790!
791! CERN PROGLIB# M433    SPACES          .VERSION KERNFOR  4.14  860211
792! ORIG.  6/05/86 M.GOOSSENS/DD
793!
794!-    The function value SPACES returns the character string STR with
795!-    leading blanks removed and each occurence of one or more blanks
796!-    replaced by NSPACE blanks inside the string STR
797!
798      CHARACTER*(80), INTENT(OUT) :: str
799      INTEGER :: nspace
800      INTEGER :: iblank, inonbl, ispace, lenspa, i, lens
801!
802      LENSPA = LEN(SPACES)
803      SPACES = ' '
804      IF (NSPACE.LT.0) NSPACE = 0
805      IBLANK = 1
806      ISPACE = 1
807  100 INONBL = INDEXC(STR(IBLANK:),' ')
808      IF (INONBL.EQ.0) THEN
809          SPACES(ISPACE:) = STR(IBLANK:)
810                                                    GO TO 999
811      ENDIF
812      INONBL = INONBL + IBLANK - 1
813      IBLANK = INDEX(STR(INONBL:),' ')
814      IF (IBLANK.EQ.0) THEN
815          SPACES(ISPACE:) = STR(INONBL:)
816                                                    GO TO 999
817      ENDIF
818      IBLANK = IBLANK + INONBL - 1
819      SPACES(ISPACE:) = STR(INONBL:IBLANK-1)
820      ISPACE = ISPACE + IBLANK - INONBL + NSPACE
821      IF (ISPACE.LE.LENSPA)                         GO TO 100
822  999 END
823
824      INTEGER FUNCTION INDEXC(STR,SSTR)
825!
826! CERN PROGLIB# M433    INDEXC          .VERSION KERNFOR  4.14  860211
827! ORIG. 26/03/86 M.GOOSSENS/DD
828!
829!-    Find the leftmost position where substring SSTR does not match
830!-    string STR scanning forward
831!
832      CHARACTER*(*), INTENT(IN) :: str, sstr
833      INTEGER :: lens, lenss, i
834!
835      LENS   = LEN(STR)
836      LENSS  = LEN(SSTR)
837!
838      DO 10 I=1,LENS-LENSS+1
839          IF (STR(I:I+LENSS-1).NE.SSTR) THEN
840              INDEXC = I
841                                         GO TO 999
842          ENDIF
843   10 CONTINUE
844      INDEXC = 0
845!
846  999 END
847
848END MODULE m_1dconv_mod_h
Note: See TracBrowser for help on using the repository browser.