/[lmdze]/trunk/libf/dyn3d/guide.f90
ViewVC logotype

Annotation of /trunk/libf/dyn3d/guide.f90

Parent Directory Parent Directory | Revision Log Revision Log


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

  ViewVC Help
Powered by ViewVC 1.1.21