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