/[lmdze]/trunk/Sources/dyn3d/Guide/guide.f
ViewVC logotype

Contents of /trunk/Sources/dyn3d/Guide/guide.f

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3 - (show annotations)
Wed Feb 27 13:16:39 2008 UTC (16 years, 2 months ago) by guez
Original Path: trunk/libf/dyn3d/guide.f90
File size: 16205 byte(s)
Initial import
1 module guide_m
2
3 ! From dyn3d/guide.F,v 1.3 2005/05/25 13:10:09
4 ! and dyn3d/guide.h,v 1.1.1.1 2004/05/19 12:53:06
5
6 real tau_min_u,tau_max_u
7 real tau_min_v,tau_max_v
8 real tau_min_T,tau_max_T
9 real tau_min_q,tau_max_q
10 real tau_min_P,tau_max_P
11 real aire_min,aire_max
12
13
14 logical guide_u,guide_v,guide_T,guide_Q,guide_P
15 real lat_min_guide,lat_max_guide
16
17 LOGICAL ncep,ini_anal
18 integer online
19
20 contains
21
22 subroutine guide(itau,ucov,vcov,teta,q,masse,ps)
23
24 use dimens_m
25 use paramet_m
26 use comconst
27 use comdissnew
28 use comvert
29 use conf_gcm_m
30 use logic
31 use comgeom
32 use serre
33 use temps
34 use tracstoke
35 use ener
36 use q_sat_m, only: q_sat
37 use exner_hyb_m, only: exner_hyb
38 use pression_m, only: pression
39 use inigrads_m, only: inigrads
40
41 IMPLICIT NONE
42
43 ! ...... Version du 10/01/98 ..........
44
45 ! avec coordonnees verticales hybrides
46 ! avec nouveaux operat. dissipation * ( gradiv2,divgrad2,nxgraro2 )
47
48 !=======================================================================
49 !
50 ! Auteur: F.Hourdin
51 ! -------
52 !
53 ! Objet:
54 ! ------
55 !
56 ! GCM LMD nouvelle grille
57 !
58 !=======================================================================
59
60 ! ... Dans inigeom , nouveaux calculs pour les elongations cu , cv
61 ! et possibilite d'appeler une fonction f(y) a derivee tangente
62 ! hyperbolique a la place de la fonction a derivee sinusoidale.
63
64 ! ... Possibilite de choisir le shema de Van-leer pour l'advection de
65 ! q , en faisant iadv = 10 dans traceur (29/04/97) .
66 !
67 !-----------------------------------------------------------------------
68 ! Declarations:
69 ! -------------
70
71 include "netcdf.inc"
72
73 ! variables dynamiques
74 REAL vcov(ip1jm,llm),ucov(ip1jmp1,llm) ! vents covariants
75 REAL teta(ip1jmp1,llm) ! temperature potentielle
76 REAL q(ip1jmp1,llm) ! temperature potentielle
77 REAL ps(ip1jmp1) ! pression au sol
78 REAL masse(ip1jmp1,llm) ! masse d'air
79
80 ! common passe pour des sorties
81 real dxdys(iip1,jjp1),dxdyu(iip1,jjp1),dxdyv(iip1,jjm)
82 common/comdxdy/dxdys,dxdyu,dxdyv
83
84 ! variables dynamiques pour les reanalyses.
85 REAL ucovrea1(ip1jmp1,llm),vcovrea1(ip1jm,llm) !vts cov reas
86 REAL tetarea1(ip1jmp1,llm) ! temp pot reales
87 REAL qrea1(ip1jmp1,llm) ! temp pot reales
88 REAL psrea1(ip1jmp1) ! ps
89 REAL ucovrea2(ip1jmp1,llm),vcovrea2(ip1jm,llm) !vts cov reas
90 REAL tetarea2(ip1jmp1,llm) ! temp pot reales
91 REAL qrea2(ip1jmp1,llm) ! temp pot reales
92 REAL masserea2(ip1jmp1,llm) ! masse
93 REAL psrea2(ip1jmp1) ! ps
94
95 real alpha_q(ip1jmp1)
96 real alpha_T(ip1jmp1),alpha_P(ip1jmp1)
97 real alpha_u(ip1jmp1),alpha_v(ip1jm)
98 real dday_step,toto,reste,itau_test
99 INTEGER step_rea,count_no_rea
100
101 !IM 180305 real aire_min,aire_max
102 integer ilon,ilat
103 real factt,ztau(ip1jmp1)
104
105 INTEGER, intent(in):: itau
106 integer ij, l
107 integer ncidpl,varidpl,nlev,status
108 integer rcod,rid
109 real ditau,tau,a
110 save nlev
111
112 ! TEST SUR QSAT
113 real p(ip1jmp1,llmp1),pk(ip1jmp1,llm),pks(ip1jmp1)
114 real pkf(ip1jmp1,llm)
115 real pres(ip1jmp1,llm)
116
117 real qsat(ip1jmp1,llm)
118 real unskap
119 real tnat(ip1jmp1,llm)
120 !cccccccccccccccc
121
122
123 LOGICAL first
124 save first
125 data first/.true./
126
127 save ucovrea1,vcovrea1,tetarea1,psrea1,qrea1
128 save ucovrea2,vcovrea2,tetarea2,masserea2,psrea2,qrea2
129
130 save alpha_T,alpha_q,alpha_u,alpha_v,alpha_P,itau_test
131 save step_rea,count_no_rea
132
133 character*10 file
134 integer igrads
135 real dtgrads
136 save igrads,dtgrads
137 data igrads,dtgrads/2,100./
138
139 print *,'Call sequence information: guide'
140
141 !-----------------------------------------------------------------------
142 ! calcul de l'humidite saturante
143 !-----------------------------------------------------------------------
144 CALL pression( ip1jmp1, ap, bp, ps, p )
145 call massdair(p,masse)
146 print*,'OK1'
147 CALL exner_hyb(ps,p,pks,pk,pkf)
148 print*,'OK2'
149 tnat(:,:)=pk(:,:)*teta(:,:)/cpp
150 print*,'OK3'
151 unskap = 1./ kappa
152 pres(:,:)=preff*(pk(:,:)/cpp)**unskap
153 print*,'OK4'
154 qsat = q_sat(tnat, pres)
155
156 !-----------------------------------------------------------------------
157
158 !-----------------------------------------------------------------------
159 ! initialisations pour la lecture des reanalyses.
160 ! alpha determine la part des injections de donnees a chaque etape
161 ! alpha=1 signifie pas d'injection
162 ! alpha=0 signifie injection totale
163 !-----------------------------------------------------------------------
164
165 print*,'ONLINE=',online
166 if(online.eq.-1) then
167 return
168 endif
169
170 if (first) then
171
172 print*,'initialisation du guide '
173 call conf_guide
174 print*,'apres conf_guide'
175
176 file='guide'
177 call inigrads(igrads &
178 ,rlonv,180./pi,-180.,180.,rlatu,-90.,90.,180./pi &
179 ,presnivs,1. &
180 ,dtgrads,file,'dyn_zon ')
181
182 print* &
183 ,'1: en-ligne, 0: hors-ligne (x=x_rea), -1: climat (x=x_gcm)'
184
185 if(online.eq.-1) return
186 if (online.eq.1) then
187
188 !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
189 ! Constantes de temps de rappel en jour
190 ! 0.1 c'est en gros 2h30.
191 ! 1e10 est une constante infinie donc en gros pas de guidage
192 !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
193 ! coordonnees du centre du zoom
194 call coordij(clon,clat,ilon,ilat)
195 ! aire de la maille au centre du zoom
196 aire_min=aire(ilon+(ilat-1)*iip1)
197 ! aire maximale de la maille
198 aire_max=0.
199 do ij=1,ip1jmp1
200 aire_max=max(aire_max,aire(ij))
201 enddo
202 ! factt = pas de temps en fraction de jour
203 factt=dtvr*iperiod/daysec
204
205 ! subroutine tau2alpha(type,im,jm,factt,taumin,taumax,alpha)
206 call tau2alpha(3,iip1,jjm ,factt,tau_min_v,tau_max_v,alpha_v)
207 call tau2alpha(2,iip1,jjp1,factt,tau_min_u,tau_max_u,alpha_u)
208 call tau2alpha(1,iip1,jjp1,factt,tau_min_T,tau_max_T,alpha_T)
209 call tau2alpha(1,iip1,jjp1,factt,tau_min_P,tau_max_P,alpha_P)
210 call tau2alpha(1,iip1,jjp1,factt,tau_min_q,tau_max_q,alpha_q)
211
212 call dump2d(iip1,jjp1,aire,'AIRE MAILLe ')
213 call dump2d(iip1,jjp1,alpha_u,'COEFF U ')
214 call dump2d(iip1,jjp1,alpha_T,'COEFF T ')
215
216 !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
217 ! Cas ou on force exactement par les variables analysees
218 else
219 alpha_T=0.
220 alpha_u=0.
221 alpha_v=0.
222 alpha_P=0.
223 ! physic=.false.
224 endif
225
226 itau_test=1001
227 step_rea=1
228 count_no_rea=0
229 ncidpl=-99
230
231 ! itau_test montre si l'importation a deja ete faite au rang itau
232 ! lecture d'un fichier netcdf pour determiner le nombre de niveaux
233 if (guide_u) then
234 if (ncidpl.eq.-99) ncidpl=NCOPN('u.nc',NCNOWRIT,rcod)
235 endif
236 !
237 if (guide_v) then
238 if (ncidpl.eq.-99) ncidpl=NCOPN('v.nc',NCNOWRIT,rcod)
239 endif
240 !
241 if (guide_T) then
242 if (ncidpl.eq.-99) ncidpl=NCOPN('T.nc',NCNOWRIT,rcod)
243 endif
244 !
245 if (guide_Q) then
246 if (ncidpl.eq.-99) ncidpl=NCOPN('hur.nc',NCNOWRIT,rcod)
247 endif
248 !
249 if (ncep) then
250 status=NF_INQ_DIMID(ncidpl,'LEVEL',rid)
251 else
252 status=NF_INQ_DIMID(ncidpl,'PRESSURE',rid)
253 endif
254 status=NF_INQ_DIMLEN(ncidpl,rid,nlev)
255 print *,'nlev', nlev
256 call ncclos(ncidpl,rcod)
257 ! Lecture du premier etat des reanalyses.
258 call read_reanalyse(1,ps &
259 ,ucovrea2,vcovrea2,tetarea2,qrea2,masserea2,psrea2,1,nlev)
260 qrea2(:,:)=max(qrea2(:,:),0.1)
261
262
263 !-----------------------------------------------------------------------
264 ! Debut de l'integration temporelle:
265 ! ----------------------------------
266
267 endif ! first
268 !
269 !-----------------------------------------------------------------------
270 !----- IMPORTATION DES VENTS,PRESSION ET TEMPERATURE REELS:
271 !-----------------------------------------------------------------------
272
273 ditau=real(itau)
274 DDAY_step=real(day_step)
275 write(*,*)'ditau,dday_step'
276 write(*,*)ditau,dday_step
277 toto=4*ditau/dday_step
278 reste=toto-aint(toto)
279 ! write(*,*)'toto,reste',toto,reste
280
281 if (reste.eq.0.) then
282 if (itau_test.eq.itau) then
283 write(*,*)'deuxieme passage de advreel a itau=',itau
284 stop
285 else
286 vcovrea1(:,:)=vcovrea2(:,:)
287 ucovrea1(:,:)=ucovrea2(:,:)
288 tetarea1(:,:)=tetarea2(:,:)
289 qrea1(:,:)=qrea2(:,:)
290
291 print*,'LECTURE REANALYSES, pas ',step_rea &
292 ,'apres ',count_no_rea,' non lectures'
293 step_rea=step_rea+1
294 itau_test=itau
295 call read_reanalyse(step_rea,ps &
296 ,ucovrea2,vcovrea2,tetarea2,qrea2,masserea2,psrea2,1,nlev)
297 qrea2(:,:)=max(qrea2(:,:),0.1)
298 factt=dtvr*iperiod/daysec
299 ztau(:)=factt/max(alpha_T(:),1.e-10)
300 call wrgrads(igrads,1,aire ,'aire ','aire ' )
301 call wrgrads(igrads,1,dxdys ,'dxdy ','dxdy ' )
302 call wrgrads(igrads,1,alpha_u,'au ','au ' )
303 call wrgrads(igrads,1,alpha_T,'at ','at ' )
304 call wrgrads(igrads,1,ztau,'taut ','taut ' )
305 call wrgrads(igrads,llm,ucov,'u ','u ' )
306 call wrgrads(igrads,llm,ucovrea2,'ua ','ua ' )
307 call wrgrads(igrads,llm,teta,'T ','T ' )
308 call wrgrads(igrads,llm,tetarea2,'Ta ','Ta ' )
309 call wrgrads(igrads,llm,qrea2,'Qa ','Qa ' )
310 call wrgrads(igrads,llm,q,'Q ','Q ' )
311
312 call wrgrads(igrads,llm,qsat,'QSAT ','QSAT ' )
313
314 endif
315 else
316 count_no_rea=count_no_rea+1
317 endif
318
319 !-----------------------------------------------------------------------
320 ! Guidage
321 ! x_gcm = a * x_gcm + (1-a) * x_reanalyses
322 !-----------------------------------------------------------------------
323
324 if(ini_anal) print*,'ATTENTION !!! ON PART DU GUIDAGE'
325
326 ditau=real(itau)
327 dday_step=real(day_step)
328
329
330 tau=4*ditau/dday_step
331 tau=tau-aint(tau)
332
333 ! ucov
334 if (guide_u) then
335 do l=1,llm
336 do ij=1,ip1jmp1
337 a=(1.-tau)*ucovrea1(ij,l)+tau*ucovrea2(ij,l)
338 ucov(ij,l)=(1.-alpha_u(ij))*ucov(ij,l)+alpha_u(ij)*a
339 if (first.and.ini_anal) ucov(ij,l)=a
340 enddo
341 enddo
342 endif
343
344 ! teta
345 if (guide_T) then
346 do l=1,llm
347 do ij=1,ip1jmp1
348 a=(1.-tau)*tetarea1(ij,l)+tau*tetarea2(ij,l)
349 teta(ij,l)=(1.-alpha_T(ij))*teta(ij,l)+alpha_T(ij)*a
350 if (first.and.ini_anal) teta(ij,l)=a
351 enddo
352 enddo
353 endif
354
355 ! P
356 if (guide_P) then
357 do ij=1,ip1jmp1
358 a=(1.-tau)*psrea1(ij)+tau*psrea2(ij)
359 ps(ij)=(1.-alpha_P(ij))*ps(ij)+alpha_P(ij)*a
360 if (first.and.ini_anal) ps(ij)=a
361 enddo
362 CALL pression(ip1jmp1,ap,bp,ps,p)
363 CALL massdair(p,masse)
364 endif
365
366
367 ! q
368 if (guide_Q) then
369 do l=1,llm
370 do ij=1,ip1jmp1
371 a=(1.-tau)*qrea1(ij,l)+tau*qrea2(ij,l)
372 ! hum relative en % -> hum specif
373 a=qsat(ij,l)*a*0.01
374 q(ij,l)=(1.-alpha_Q(ij))*q(ij,l)+alpha_Q(ij)*a
375 if (first.and.ini_anal) q(ij,l)=a
376 enddo
377 enddo
378 endif
379
380 ! vcov
381 if (guide_v) then
382 do l=1,llm
383 do ij=1,ip1jm
384 a=(1.-tau)*vcovrea1(ij,l)+tau*vcovrea2(ij,l)
385 vcov(ij,l)=(1.-alpha_v(ij))*vcov(ij,l)+alpha_v(ij)*a
386 if (first.and.ini_anal) vcov(ij,l)=a
387 enddo
388 if (first.and.ini_anal) vcov(ij,l)=a
389 enddo
390 endif
391
392 ! call dump2d(iip1,jjp1,tetarea1,'TETA REA 1 ')
393 ! call dump2d(iip1,jjp1,tetarea2,'TETA REA 2 ')
394 ! call dump2d(iip1,jjp1,teta,'TETA ')
395
396 first=.false.
397
398 return
399 end subroutine guide
400
401 !=======================================================================
402 subroutine tau2alpha(type,pim,pjm,factt,taumin,taumax,alpha)
403 !=======================================================================
404
405 use dimens_m
406 use paramet_m
407 use comconst, only: pi
408 use comgeom
409 use serre
410 implicit none
411
412 ! arguments :
413 integer type
414 integer pim,pjm
415 real factt,taumin,taumax
416 real dxdy_,alpha(pim,pjm)
417 real dxdy_min,dxdy_max
418
419 ! local :
420 real alphamin,alphamax,gamma,xi
421 save gamma
422 integer i,j,ilon,ilat
423
424 logical first
425 save first
426 data first/.true./
427
428 real zdx(iip1,jjp1),zdy(iip1,jjp1)
429
430 real zlat
431 real dxdys(iip1,jjp1),dxdyu(iip1,jjp1),dxdyv(iip1,jjm)
432 common/comdxdy/dxdys,dxdyu,dxdyv
433
434 if (first) then
435 do j=2,jjm
436 do i=2,iip1
437 zdx(i,j)=0.5*(cu_2d(i-1,j)+cu_2d(i,j))/cos(rlatu(j))
438 enddo
439 zdx(1,j)=zdx(iip1,j)
440 enddo
441 do j=2,jjm
442 do i=1,iip1
443 zdy(i,j)=0.5*(cv_2d(i,j-1)+cv_2d(i,j))
444 enddo
445 enddo
446 do i=1,iip1
447 zdx(i,1)=zdx(i,2)
448 zdx(i,jjp1)=zdx(i,jjm)
449 zdy(i,1)=zdy(i,2)
450 zdy(i,jjp1)=zdy(i,jjm)
451 enddo
452 do j=1,jjp1
453 do i=1,iip1
454 dxdys(i,j)=sqrt(zdx(i,j)*zdx(i,j)+zdy(i,j)*zdy(i,j))
455 enddo
456 enddo
457 do j=1,jjp1
458 do i=1,iim
459 dxdyu(i,j)=0.5*(dxdys(i,j)+dxdys(i+1,j))
460 enddo
461 dxdyu(iip1,j)=dxdyu(1,j)
462 enddo
463 do j=1,jjm
464 do i=1,iip1
465 dxdyv(i,j)=0.5*(dxdys(i,j)+dxdys(i+1,j))
466 enddo
467 enddo
468
469 call dump2d(iip1,jjp1,dxdys,'DX2DY2 SCAL ')
470 call dump2d(iip1,jjp1,dxdyu,'DX2DY2 U ')
471 call dump2d(iip1,jjp1,dxdyv,'DX2DY2 v ')
472
473 ! coordonnees du centre du zoom
474 call coordij(clon,clat,ilon,ilat)
475 ! aire de la maille au centre du zoom
476 dxdy_min=dxdys(ilon,ilat)
477 ! dxdy maximale de la maille
478 dxdy_max=0.
479 do j=1,jjp1
480 do i=1,iip1
481 dxdy_max=max(dxdy_max,dxdys(i,j))
482 enddo
483 enddo
484
485 if (abs(grossismx-1.).lt.0.1.or.abs(grossismy-1.).lt.0.1) then
486 print*,'ATTENTION modele peu zoome'
487 print*,'ATTENTION on prend une constante de guidage cste'
488 gamma=0.
489 else
490 gamma=(dxdy_max-2.*dxdy_min)/(dxdy_max-dxdy_min)
491 print*,'gamma=',gamma
492 if (gamma.lt.1.e-5) then
493 print*,'gamma =',gamma,'<1e-5'
494 stop
495 endif
496 print*,'gamma=',gamma
497 gamma=log(0.5)/log(gamma)
498 endif
499 endif
500
501 alphamin=factt/taumax
502 alphamax=factt/taumin
503
504 do j=1,pjm
505 do i=1,pim
506 if (type.eq.1) then
507 dxdy_=dxdys(i,j)
508 zlat=rlatu(j)*180./pi
509 elseif (type.eq.2) then
510 dxdy_=dxdyu(i,j)
511 zlat=rlatu(j)*180./pi
512 elseif (type.eq.3) then
513 dxdy_=dxdyv(i,j)
514 zlat=rlatv(j)*180./pi
515 endif
516 if (abs(grossismx-1.).lt.0.1.or.abs(grossismy-1.).lt.0.1) then
517 ! pour une grille reguliere, xi=xxx**0=1 -> alpha=alphamin
518 alpha(i,j)=alphamin
519 else
520 xi=((dxdy_max-dxdy_)/(dxdy_max-dxdy_min))**gamma
521 xi=min(xi,1.)
522 if(lat_min_guide.le.zlat .and. zlat.le.lat_max_guide) then
523 alpha(i,j)=xi*alphamin+(1.-xi)*alphamax
524 else
525 alpha(i,j)=0.
526 endif
527 endif
528 enddo
529 enddo
530
531
532 return
533 end subroutine tau2alpha
534
535 end module guide_m

  ViewVC Help
Powered by ViewVC 1.1.21