source: codes/icosagcm/branches/SATURN_DYNAMICO/LMDZ.COMMON/libf/phystd/mufract.F @ 263

Last change on this file since 263 was 222, checked in by ymipsl, 10 years ago

Creating temporary dynamico/lmdz/saturn branche

YM
File size: 2.1 KB
Line 
1      SUBROUTINE mufract(jjm,pdecli, plat, pmu, pfract)
2      IMPLICIT NONE
3c
4c=======================================================================
5c
6c   Calcul of equivalent solar angle and and fraction of day without 
7c   diurnal cycle.
8c
9c   parmeters :
10c   -----------
11c
12c      Input :
13c      -------
14c         jjm              number of points
15c         pdecli           solar declinaison
16c         plat(jjm)        latitude
17c
18c      Output :
19c      --------
20c         pmu(jjm)          equivalent cosinus of the solar angle
21c         pfract(jjm)       fractionnal day
22c
23c
24c=======================================================================
25c
26c-----------------------------------------------------------------------
27c
28c    0. Declarations :
29c    -----------------
30c
31c     Arguments :
32c     -----------
33
34      INTEGER jjm
35      REAL plat(jjm),pmu(jjm), pfract(jjm)
36      REAL pdecli
37c
38c     Local variables :
39c     -----------------
40
41      INTEGER j
42      REAL pi
43      REAL z,cz,sz,tz,phi,cphi,sphi,tphi
44      REAL ap,a,t,b
45c
46c=======================================================================
47c
48      pi = 4. * atan(1.0)
49      z = pdecli
50      cz = cos (z*pi/180.)
51      sz = sin (z*pi/180.)
52c
53      DO 20 j = 1, jjm
54c
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         IF(t.eq.0.) THEN
65            t=0.5*pi
66         ELSE
67            IF (a.lt.0.) a = 0.
68            t = sqrt(a) / t
69            IF (t.lt.0.) THEN
70               t = -atan (-t) + pi
71            ELSE
72               t = atan(t)
73            ENDIF
74         ENDIF
75         pmu(j) = (sphi*sz*t) / pi + b*sin(t)/pi
76         pfract(j) = t / pi
77         IF (ap .lt.0.) THEN
78            pmu(j) = sphi * sz
79            pfract(j) = 1.0
80         ENDIF
81         IF (pmu(j).le.0.0) pmu(j) = 0.0
82         pmu(j) = pmu(j) / pfract(j)
83         IF (pmu(j).eq.0.) pfract(j) = 0.
84c
85   20 CONTINUE
86c
87      RETURN
88      END
89c
90c* end of mufract
91c=======================================================================
92c
Note: See TracBrowser for help on using the repository browser.