Changeset 1161 for trunk/LMDZ.GENERIC/libf/phystd/rings.F90

Timestamp:

Jan 21, 2014, 4:07:42 PM (11 years ago)

Author:

msylvestre

Message:

LMDZ.GENERIC/UNIVERSAL. rings shadow added in the case diurnal=.false. (daily average of the shadow)

File:

• trunk/LMDZ.GENERIC/libf/phystd/rings.F90 (modified) (6 diffs)

trunk/LMDZ.GENERIC/libf/phystd/rings.F90

@@ -40 +40 @@
     REAL, DIMENSION(:), ALLOCATABLE:: x, y, z ! cartesian coordinates of the points on the planet
     REAL :: xs, ys, zs                        ! cartesian coordinates of the points of the subsolar point
-    REAL, DIMENSION(:), ALLOCATABLE :: k      ! parameter (?)
+    REAL, DIMENSION(:), ALLOCATABLE :: k
     REAL, DIMENSION(:), ALLOCATABLE :: N      ! parameter to compute cartesian coordinates on a ellipsoidal planet
     REAL, DIMENSION(:), ALLOCATABLE :: r      ! distance at which the incident ray of sun crosses the equatorial plane
@@ -48 +48 @@
     ! equinox --> no shadow (AS: why is this needed?)
     if(declin .eq. 0.) then
-        eclipse = 0.
+        eclipse(:) = 0.
         return 
     endif 
@@ -63 +63 @@
     ALLOCATE(N(ngrid))
     ALLOCATE(r(ngrid))
-    eclipse = 2000.
+    eclipse(:) = 2000.
 
     ! Size of the rings
@@ -85 +85 @@
 
     ! Convert to cartesian coordinates
-    N = rad/sqrt(1-(e**2)*sinlat**2)
-    x = N*coslat*coslon
-    y = N*coslat*sinlon
-    z = N*(1-e**2)*sinlat
+    N(:) = rad / sqrt(1-(e**2)*sinlat(:)**2)
+    x(:) = N(:)*coslat(:)*coslon(:)
+    y(:) = N(:)*coslat(:)*sinlon(:)
+    z(:) = N(:)*(1-e**2)*sinlat(:)
 
 ! 2) LOCATION OF THE SUBSOLAR POINT 
@@ -95 +95 @@
     ! SG: the minus sign is important! ... otherwise subsolar point adopts a reverse rotation
     phi_S = -(ptime - 0.5)*2.*pi 
-    write(*,*) 'subsol point coords : ', declin*180./pi, phi_S*180./pi
+!    write(*,*) 'subsol point coords : ', declin*180./pi, phi_S*180./pi
 
     ! subsolar latitude is declin (declination of the sun)
@@ -106 +106 @@
 ! 3) WHERE DOES THE INCIDENT RAY OF SUN CROSS THE EQUATORIAL PLAN ?
 
-    k = -z/zs
-    r = (k*xs + x)**2 + (k*ys + y)**2 
-    r = sqrt(r)
+    k(:) = -z(:)/zs
+    r(:) = (k(:)*xs + x(:))**2 + (k(:)*ys + y(:))**2 
+    r(:) = sqrt(r(:))
 
 ! 4) SO WHERE ARE THE SHADOW OF THESE RINGS ?
 
     ! Summer hemisphere is not under the shadow of the rings
-    where(lati*declin .gt. 0.)
-       eclipse = 1000.
+    where(lati(:)*declin .gt. 0.)
+       eclipse(:) = 1000.
     end where
 
     ! No shadow of the rings by night
-    where(x*xs + y*ys + z*zs .lt. 0.)
-       eclipse = 1000.
+    where(x(:)*xs + y(:)*ys + z(:)*zs .lt. 0.)
+       eclipse(:) = 1000.
     end where
 
     ! if the incident rays of sun cross a ring, there is a shadow
     do i=1, nb_A 
-        where(r .ge. A_Rint(i) .and. r .le. A_Rext(i) .and. eclipse .ne. 1000.)
-            eclipse = 1. - exp(-tau_A(i)/cos(declin))
+        where(r(:) .ge. A_Rint(i) .and. r(:) .le. A_Rext(i) .and. eclipse(:) .ne. 1000.)
+            eclipse(:) = 1. - exp(-tau_A(i)/cos(declin))
         end where
     end do 
 
     do i=1, nb_B 
-        where(r .ge. B_Rint(i) .and. r .le. B_Rext(i) .and. eclipse .ne. 1000.)
-            eclipse = 1. - exp(-tau_B(i)/cos(declin))
+        where(r(:) .ge. B_Rint(i) .and. r(:) .le. B_Rext(i) .and. eclipse(:) .ne. 1000.)
+            eclipse(:) = 1. - exp(-tau_B(i)/cos(declin))
         end where
     enddo
     
     do i=1, nb_C 
-        where(r .ge. C_Rint(i) .and. r .le. C_Rext(i) .and. eclipse .ne. 1000.)
-            eclipse = 1. - exp(-tau_C(i)/cos(declin))
+        where(r(:) .ge. C_Rint(i) .and. r(:) .le. C_Rext(i) .and. eclipse(:) .ne. 1000.)
+            eclipse(:) = 1. - exp(-tau_C(i)/cos(declin))
         end where
     enddo
 
     ! At the other places and the excluded ones, eclipse is 0. 
-    where(eclipse .eq. 2000. .or. eclipse .eq. 1000.)
-        eclipse = 0. 
+    where(eclipse(:) .eq. 2000. .or. eclipse(:) .eq. 1000.)
+        eclipse(:) = 0. 
     end where