source: trunk/LMDZ.GENERIC/libf/phystd/rings.F90 @ 1579

Last change on this file since 1579 was 1543, checked in by emillour, 9 years ago

All models: Further adaptations to keep up with changes in LMDZ5 concerning
physics/dynamics separation:

  • dyn3d:
  • adapted gcm.F so that all physics initializations are now done in iniphysiq.
  • dyn3dpar:
  • adapted gcm.F so that all physics initializations are now done in iniphysiq.
  • updated calfis_p.F to follow up with changes.
  • copied over updated "bands.F90" from LMDZ5.
  • dynphy_lonlat:
  • calfis_p.F90, mod_interface_dyn_phys.F90, follow up of changes in phy_common/mod_* routines
  • phy_common:
  • added "geometry_mod.F90" to store information about the grid (replaces phy*/comgeomphy.F90) and give variables friendlier names: rlond => longitude , rlatd => latitude, airephy => cell_area, cuphy => dx , cvphy => dy
  • added "physics_distribution_mod.F90"
  • updated "mod_grid_phy_lmdz.F90", "mod_phys_lmdz_mpi_data.F90", "mod_phys_lmdz_para.F90", "mod_phys_lmdz_mpi_transfert.F90", "mod_grid_phy_lmdz.F90", "mod_phys_lmdz_omp_data.F90", "mod_phys_lmdz_omp_transfert.F90", "write_field_phy.F90" and "ioipsl_getin_p_mod.F90" to LMDZ5 versions.
  • phy[venus/titan/mars/std]:
  • removed "init_phys_lmdz.F90", "comgeomphy.F90"; adapted routines to use geometry_mod (longitude, latitude, cell_area, etc.)

EM

File size: 6.8 KB
Line 
1SUBROUTINE rings(ngrid, declin, ptime, rad, flat, eclipse)
2! Calculates Saturn's rings shadowing
3! Includes rings opacities measured by Cassini/UVIS
4! Authors: M. Sylvestre, M. Capderou, S. Guerlet, A. Spiga
5
6    use comdiurn_h, only: sinlat, sinlon, coslat, coslon
7    use geometry_mod, only: latitude ! (rad)
8 
9    implicit none   
10
11    INTEGER, INTENT(IN) :: ngrid  ! horizontal grid dimension
12    REAL, INTENT(IN) :: declin    ! latitude of the subsolar point
13    REAL, INTENT(IN) :: ptime     ! UTC time in sol fraction : ptime=0.5 at noon
14    REAL, INTENT(IN) :: rad       ! equatorial radius of the planet
15    REAL, INTENT(IN) :: flat      ! flattening of the planet
16    REAL, DIMENSION(ngrid), INTENT(OUT) :: eclipse ! absorption of the light by the rings   
17   
18    REAL :: rpol   ! polar radius of the planet
19    REAL :: e      ! shape excentricity of the planet : (1-e*e) = (1-f)*(1-f)   
20    INTEGER, PARAMETER :: nb_a = 4 ! number of subdivisions of the A ring
21    INTEGER, PARAMETER :: nb_b = 3 ! number of subdivisions of the B ring
22    INTEGER, PARAMETER :: nb_c = 3 ! number of subdivisions of the C ring
23    INTEGER, PARAMETER :: nb_ca = 2 ! number of subdivisions in the Cassini division
24    INTEGER :: i
25
26    ! arrays for the rings. TBD: dynamical?
27    REAL, DIMENSION(nb_a) :: A_Rint ! internal radii of the subdivisions of the A ring
28    REAL, DIMENSION(nb_a) :: A_Rext ! external radii of the subdivisions of the A ring
29    REAL, DIMENSION(nb_b) :: B_Rint ! internal radii of the subdivisions of the B ring
30    REAL, DIMENSION(nb_b) :: B_Rext ! external radii of the subdivisions of the B ring
31    REAL, DIMENSION(nb_c) :: C_Rint ! internal radii of the subdivisions of the C ring
32    REAL, DIMENSION(nb_c) :: C_Rext ! external radii of the subdivisions of the C ring
33    REAL, DIMENSION(nb_ca) :: Ca_Rint ! internal radii of the subdivisions of the Cassini Division
34    REAL, DIMENSION(nb_ca) :: Ca_Rext ! external radii of the subdivisions of the Cassini Division
35
36    ! Opacities of the rings : for each one we can give different opacities for each part
37    REAL, DIMENSION(nb_a) :: tau_A ! opacity of the A ring
38    REAL, DIMENSION(nb_b) :: tau_B ! opacity of the B ring
39    REAL, DIMENSION(nb_c) :: tau_C ! opacity of the C ring
40    REAL, DIMENSION(nb_ca) :: tau_Ca ! opacity of the Cassini Division
41
42    ! Parameters used to calculate if a point is under a ring subdivision's shadow
43    REAL :: phi_S                             ! subsolar point longitude
44    REAL, PARAMETER :: pi=acos(-1.0)   
45    REAL, DIMENSION(:), ALLOCATABLE:: x, y, z ! cartesian coordinates of the points on the planet
46    REAL :: xs, ys, zs                        ! cartesian coordinates of the points of the subsolar point
47    REAL, DIMENSION(:), ALLOCATABLE :: k
48    REAL, DIMENSION(:), ALLOCATABLE :: N      ! parameter to compute cartesian coordinates on a ellipsoidal planet
49    REAL, DIMENSION(:), ALLOCATABLE :: r      ! distance at which the incident ray of sun crosses the equatorial plane
50                                              ! measured from the center of the planet   
51    REAL :: Ns                                ! (same for the subsolar point)
52   
53    ! equinox --> no shadow (AS: why is this needed?)
54    if(declin .eq. 0.) then
55        eclipse(:) = 0.
56        return
57    endif
58
59! 1) INITIALIZATION
60
61    ! Generic
62    rpol = (1.- flat)*rad
63    e = sqrt(2*flat - flat**2)
64    ALLOCATE(x(ngrid))
65    ALLOCATE(y(ngrid))
66    ALLOCATE(z(ngrid))
67    ALLOCATE(k(ngrid))
68    ALLOCATE(N(ngrid))
69    ALLOCATE(r(ngrid))
70    eclipse(:) = 2000.
71
72! Model of the rings with Cassini/UVIS opacities
73
74    ! Size of the rings
75    A_Rint(1) = 2.03*rad
76    A_Rext(1) = 2.06*rad
77    A_Rint(2) = 2.06*rad
78    A_Rext(2) = 2.09*rad
79    A_Rint(3) = 2.09*rad
80    A_Rext(3) = 2.12*rad
81    A_Rint(4) = 2.12*rad
82    A_Rext(4) = 2.27*rad
83
84    B_Rint(1) = 1.53*rad
85    B_Rext(1) = 1.64*rad
86    B_Rint(2) = 1.64*rad
87    B_Rext(2) = 1.83*rad
88    B_Rint(3) = 1.83*rad
89    B_Rext(3) = 1.95*rad
90   
91    C_Rint(1) = 1.24*rad
92    C_Rext(1) = 1.29*rad
93    C_Rint(2) = 1.29*rad
94    C_Rext(2) = 1.43*rad
95    C_Rint(3) = 1.43*rad
96    C_Rext(3) = 1.53*rad
97
98    Ca_Rint(1) = 1.95*rad
99    Ca_Rext(1) = 1.99*rad
100    Ca_Rint(2) = 1.99*rad
101    Ca_Rext(2) = 2.03*rad
102
103
104    ! Opacities of the rings
105    tau_A(1) = 1.24
106    tau_A(2) = 0.81
107    tau_A(3) = 0.67
108    tau_A(4) = 0.58
109               
110    tau_B(1) = 1.29
111    tau_B(2) = 5.13
112    tau_B(3) = 2.84
113   
114    tau_C(1) = 0.06
115    tau_C(2) = 0.10
116    tau_C(3) = 0.14
117
118    tau_Ca(1) = 0.06
119    tau_Ca(2) = 0.24
120
121    ! Convert to cartesian coordinates
122    N(:) = rad / sqrt(1-(e**2)*sinlat(:)**2)
123    x(:) = N(:)*coslat(:)*coslon(:)
124    y(:) = N(:)*coslat(:)*sinlon(:)
125    z(:) = N(:)*(1-e**2)*sinlat(:)
126
127! 2) LOCATION OF THE SUBSOLAR POINT
128 
129    ! subsolar longitude is deduced from time fraction ptime
130    ! SG: the minus sign is important! ... otherwise subsolar point adopts a reverse rotation
131    phi_S = -(ptime - 0.5)*2.*pi
132!    write(*,*) 'subsol point coords : ', declin*180./pi, phi_S*180./pi
133
134    ! subsolar latitude is declin (declination of the sun)
135    ! now convert in cartesian coordinates :
136    Ns = rad/sqrt(1-(e**2)*sin(declin)**2)
137    xs = Ns*cos(declin)*cos(phi_S)
138    ys = Ns*cos(declin)*sin(phi_S)
139    zs = Ns*(1-e**2)*sin(declin)
140
141! 3) WHERE DOES THE INCIDENT RAY OF SUN CROSS THE EQUATORIAL PLAN ?
142
143    k(:) = -z(:)/zs
144    r(:) = (k(:)*xs + x(:))**2 + (k(:)*ys + y(:))**2
145    r(:) = sqrt(r(:))
146
147! 4) SO WHERE ARE THE SHADOW OF THESE RINGS ?
148
149    ! Summer hemisphere is not under the shadow of the rings
150    where(latitude(:)*declin .gt. 0.)
151       eclipse(:) = 1000.
152    end where
153
154    ! No shadow of the rings by night
155    where(x(:)*xs + y(:)*ys + z(:)*zs .lt. 0.)
156       eclipse(:) = 1000.
157    end where
158
159    ! if the incident rays of sun cross a ring, there is a shadow
160    do i=1, nb_A
161        where(r(:) .ge. A_Rint(i) .and. r(:) .le. A_Rext(i) .and. eclipse(:) .ne. 1000.)
162            eclipse(:) = 1. - exp(-tau_A(i)/abs(sin(declin)))
163        end where
164    end do
165
166    do i=1, nb_B
167        where(r(:) .ge. B_Rint(i) .and. r(:) .le. B_Rext(i) .and. eclipse(:) .ne. 1000.)
168            eclipse(:) = 1. - exp(-tau_B(i)/abs(sin(declin)))
169        end where
170    enddo
171   
172    do i=1, nb_C
173        where(r(:) .ge. C_Rint(i) .and. r(:) .le. C_Rext(i) .and. eclipse(:) .ne. 1000.)
174            eclipse(:) = 1. - exp(-tau_C(i)/abs(sin(declin)))
175        end where
176    enddo
177
178    do i=1, nb_ca
179        where(r(:) .ge. Ca_Rint(i) .and. r(:) .le. Ca_Rext(i) .and. eclipse(:) .ne. 1000.)
180            eclipse(:) = 1. - exp(-tau_Ca(i)/abs(sin(declin)))
181        end where
182    enddo
183
184    ! At the other places and the excluded ones, eclipse is 0.
185    where(eclipse(:) .eq. 2000. .or. eclipse(:) .eq. 1000.)
186        eclipse(:) = 0.
187    end where
188
189! 5) CLEAN THE PLACE
190    DEALLOCATE(x)
191    DEALLOCATE(y)
192    DEALLOCATE(z)
193    DEALLOCATE(k)
194    DEALLOCATE(N)
195    DEALLOCATE(r)
196
197END SUBROUTINE rings
Note: See TracBrowser for help on using the repository browser.