source: trunk/LMDZ.TITAN/libf/phytitan/call_rings.F90

Last change on this file was 1429, checked in by msylvestre, 10 years ago

Some computations relative to the rings' shadow have been moved from physiq.f90 to call_rings.f90

File size: 2.7 KB
Line 
1subroutine call_rings(ngrid, ptime, pday, diurnal)
2    ! A subroutine to compute the day fraction in case of rings shadowing.
3
4    use radcommon_h, only: eclipse
5    use comsaison_h, only: fract, declin
6    use comcstfi_mod, only: rad, pi
7    use comdiurn_h, only: coslat, sinlat, coslon, sinlon
8    use callkeys_mod, only: flatten
9
10    INTEGER, INTENT(IN) :: ngrid
11    REAL, INTENT(IN) :: ptime ! "universal time", given as fraction of sol (e.g.: 0.5 for noon)
12    REAL, INTENT(IN) :: pday  ! Number of days counted from the North. Spring equinoxe.
13    LOGICAL, INTENT(IN) :: diurnal
14!    REAL, DIMENSION(:), INTENT(INOUT) :: fract ! day fraction for each point of the planet
15
16!   to compute the daily average of rings shadowing
17    INTEGER, PARAMETER :: nb_hours = 1536 ! set how many times per day are used
18    REAL :: pas
19    INTEGER :: m
20    REAL :: ptime_day ! Universal time in sol fraction
21    REAL:: tmp_zls,tmp_dist_star, tmp_declin, tmp_right_ascen   ! tmp solar longitude, stellar dist, declin and RA
22    REAL :: ztim1, ztim2, ztim3
23    REAL, DIMENSION(:), ALLOCATABLE :: tmp_fract ! day fraction of the time interval
24    REAL, DIMENSION(:), ALLOCATABLE :: tmp_mu0 ! equivalent solar angle
25
26!! Eclipse incoming sunlight (e.g. Saturn ring shadowing)
27    ALLOCATE(eclipse(ngrid))
28
29    write(*,*) 'Rings shadow activated'
30       
31    if(diurnal .eqv. .false.) then ! we need to compute the daily average insolation (day fraction)
32        pas = 1./nb_hours
33        ptime_day = 0.
34        fract(:) = 0.
35        ALLOCATE(tmp_fract(ngrid))
36        ALLOCATE(tmp_mu0(ngrid))
37        tmp_fract(:) = 0.
38        eclipse(:) = 0.
39        tmp_mu0(:) = 0.
40                   
41        do m=1, nb_hours
42            ptime_day = m*pas
43            call stellarlong(pday+ptime_day,tmp_zls)
44            call orbite(tmp_zls,tmp_dist_star,tmp_declin,tmp_right_ascen)
45           
46            ztim1=SIN(tmp_declin)
47            ztim2=COS(tmp_declin)*COS(2.*pi*(pday+ptime_day-.5))
48            ztim3=-COS(tmp_declin)*SIN(2.*pi*(pday+ptime_day-.5))
49
50            call stelang(ngrid,sinlon,coslon,sinlat,coslat,    &
51                        ztim1,ztim2,ztim3,tmp_mu0,tmp_fract, flatten)       
52            call rings(ngrid, tmp_declin, ptime_day, rad, flatten, eclipse)
53            fract(:) = fract(:) + (1.-eclipse(:))*tmp_fract(:) !! fract takes into account the rings shadow and the day/night alternation
54
55        enddo       
56     
57        fract(:) = fract(:)/nb_hours
58
59        DEALLOCATE(tmp_fract)
60        DEALLOCATE(tmp_mu0)
61                 
62     else   ! instant insolation is weighted by the rings shadow
63            call rings(ngrid, declin, ptime, rad, 0., eclipse)
64            fract(:) = fract(:) * (1.-eclipse)
65    endif
66
67    IF (ALLOCATED(eclipse)) DEALLOCATE(eclipse)
68
69end subroutine call_rings
Note: See TracBrowser for help on using the repository browser.