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