source: trunk/LMDZ.GENERIC/libf/phystd/mucorr.F @ 603

Last change on this file since 603 was 135, checked in by aslmd, 14 years ago

CHANGEMENT ARBORESCENCE ETAPE 2 -- NON COMPLET

File size: 2.5 KB
Line 
1      SUBROUTINE mucorr(npts,pdeclin, plat, pmu, pfract,phaut,prad)
2      IMPLICIT NONE
3
4c=======================================================================
5c
6c   Calcul of equivalent solar angle and and fraction of day whithout
7c   diurnal cycle.
8c
9c   parmeters :
10c   -----------
11c
12c      Input :
13c      -------
14c         npts             number of points
15c         pdeclin          solar declinaison
16c         plat(npts)        latitude
17c         phaut            hauteur typique de l'atmosphere
18c         prad             rayon planetaire
19c
20c      Output :
21c      --------
22c         pmu(npts)          equivalent cosinus of the solar angle
23c         pfract(npts)       fractionnal day
24c
25c=======================================================================
26
27c-----------------------------------------------------------------------
28c
29c    0. Declarations :
30c    -----------------
31
32c     Arguments :
33c     -----------
34      INTEGER npts
35      REAL plat(npts),pmu(npts), pfract(npts)
36      REAL phaut,prad,pdeclin
37c
38c     Local variables :
39c     -----------------
40      INTEGER j
41      REAL pi
42      REAL z,cz,sz,tz,phi,cphi,sphi,tphi
43      REAL ap,a,t,b, tp
44      REAL alph
45
46c-----------------------------------------------------------------------
47
48      pi = 4. * atan(1.0)
49      z = pdeclin
50      cz = cos (z)
51      sz = sin (z)
52
53      DO 20 j = 1, npts
54
55         phi = plat(j)
56         cphi = cos(phi)
57         if (cphi.le.1.e-9) cphi=1.e-9
58         sphi = sin(phi)
59         tphi = sphi / cphi
60         b = cphi * cz
61         t = -tphi * sz / cz
62         a = 1.0 - t*t
63         ap = a
64   
65         IF(t.eq.0.) then
66            tp=0.5*pi
67         ELSE
68            IF (a.lt.0.) a = 0.
69            t = sqrt(a) / t
70            IF (t.lt.0.) then
71               tp = -atan (-t) + pi
72            ELSE
73               tp = atan(t)
74            ENDIF
75         ENDIF
76         t = tp
77   
78         pmu(j) = (sphi*sz*t) / pi + b*sin(t)/pi
79         pfract(j) = t / pi
80         IF (ap .lt.0.) then
81            pmu(j) = sphi * sz
82            pfract(j) = 1.0
83         ENDIF
84         IF (pmu(j).le.0.0) pmu(j) = 0.0
85         pmu(j) = pmu(j) / pfract(j)
86         IF (pmu(j).eq.0.) pfract(j) = 0.
87
88   20 CONTINUE
89
90c-----------------------------------------------------------------------
91c   correction de rotondite:
92c   ------------------------
93
94      alph=phaut/prad
95      DO 30 j=1,npts
96c !!!!!!
97         pmu(j)=sqrt(1224.*pmu(j)*pmu(j)+1.)/35.
98c    $          (sqrt(alph*alph*pmu(j)*pmu(j)+2.*alph+1.)-alph*pmu(j))
9930    CONTINUE
100
101      RETURN
102      END
Note: See TracBrowser for help on using the repository browser.