source: LMDZ6/branches/DYNAMICO-conv/libf/phymar/PHY_Atm_S0_RUN.f90 @ 3106

Last change on this file since 3106 was 2089, checked in by Laurent Fairhead, 10 years ago

Inclusion de la physique de MAR


Integration of MAR physics

File size: 17.4 KB
Line 
1      subroutine PHY_Atm_S0_RUN
2
3!------------------------------------------------------------------------------+
4!                                                         Mon 17-Jun-2013  MAR |
5!   MAR          PHY_Atm_S0_RUN                                                |
6!     subroutine PHY_Atm_S0_RUN performs    Insolation Computation             |
7!                                                                              |
8!     version 3.p.4.1 created by H. Gallee,               Thu 25-Apr-2013      |
9!           Last Modification by H. Gallee,               Mon 17-Jun-2013      |
10!                                                                              |
11!------------------------------------------------------------------------------+
12!                                                                              |
13!     REFER.:   Ch.  Tricot, personal communication                            |
14!     ^^^^^^^   M.F. Loutre, personal communication and thesis (1993)          |
15!                                                                              |
16!     INPUT :   Mon_TU, Day_TU        : Month and Day of the Year              |
17!     ^^^^^^^   HourTU, MinuTU, Sec_TU: Hour, Minute, and Second               |
18!               lon__r(kcolp) : Latitude                             [radians] |
19!               lat__h(kcolp) : Longitude                              [hours] |
20!                                                                              |
21!     OUTPUT:   rsunS0        : Insolation normal to Atmosphere Top (W/m2)     |
22!     ^^^^^^^   csz0S0        : Cosinus of the Zenithal Distance               |
23!                                                                              |
24!------------------------------------------------------------------------------+
25
26
27! Global Variables
28! ================
29
30      use Mod_Real
31      use Mod_PHY____dat
32      use Mod_PHY____grd
33      use Mod_PHY____kkl
34      use Mod_PHY_S0_ctr
35      use Mod_PHY_S0_dat
36      use Mod_PHY_S0_grd
37      use Mod_PHY_S0_kkl
38
39
40
41! LOCAL VARIABLES
42! ===============
43
44      use Mod_Atm_S0_RUN
45
46
47      IMPLICIT NONE
48
49
50      integer                                       ::   i, j  !
51      integer                                       ::  ikl    !
52      integer                                       ::     nj  !
53      integer                                       ::  lhc    !
54
55      real(kind=real8)                              ::  dlamm  ! mean long. sun for ma-ja
56      real(kind=real8)                              ::    anm  !
57      real(kind=real8)                              ::   ranm  !
58      real(kind=real8)                              ::   ranv  ! anomalie  vraie                          [radian]
59      real(kind=real8)                              ::    anv  ! anomalie  vraie                          [degree]
60      real(kind=real8)                              ::    tls  ! longitude vraie                          [degree]
61      real(kind=real8)                              ::  rlam   ! longitude vraie                          [radian]
62      real(kind=real8)                              ::    sd   ! sinus of solar declination angle (delta)      [-]
63      real(kind=real8)                              ::    cd   ! cosin of solar declination angle (delta)      [-]
64      real(kind=real8)                              ::  deltar ! Solar Declination (angle sun at equator) [radian]
65      real(kind=real8)                              ::  delta  ! Solar Declination (angle sun at equator) [degree]
66
67      real(kind=real8)                              ::  ddt    !
68      real(kind=real8)                              ::  arg    !
69      real(kind=real8)                              ::  et     !
70      real(kind=real8)                              ::  c1     !
71      real(kind=real8)                              ::  c2     !
72      real(kind=real8)                              ::  c3     !
73      real(kind=real8)                              ::  s1     !
74      real(kind=real8)                              ::  s2     !
75      real(kind=real8)                              ::  s3     !
76
77      real(kind=real8)                              ::  argc   !
78      real(kind=real8)                              ::  ahc    !
79
80      real(kind=real8)                              ::  LenDay !  Day     Length                               [h]
81      real(kind=real8)                              ::  SunRis !  Time of Sun Rise                             [h]
82      real(kind=real8)                              ::  SunSet !  Time of Sun Set                              [h]
83      real(kind=real8)                              ::  timh   !
84      real(kind=real8)                              ::  ahor   !  Time Angle                                [hour]
85      real(kind=real8)                              ::  ahorr  !  Time Angle                              [radian]
86      real(kind=real8)                              ::  chor   !
87      real(kind=real8)                              ::  zenitr !                                                   RUN
88
89      real(kind=real8)                              ::  omesun !  Sun       Azimuth                       [radian] RUN
90
91      real(kind=real8)                              ::  comes  !
92      real(kind=real8)                              ::  somes  !
93      real(kind=real8)                              ::  omeswr !
94
95      real(kind=real8)                              ::  r_azim !  azimuth         , fine resolution                ALL
96
97      integer                                       ::  k_azim !                                                   ALL
98      integer                                       ::  knazim !                                                   RUN
99
100
101
102
103!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
104!                                                                  !
105! ALLOCATION
106! ==========
107
108      IF (it_RUN.EQ.1 .OR. FlagDALLOC)                        THEN !
109          allocate ( cszMIN(kcolp) )                               !  MIN    cos(zenith.Dist), if Mountain MASK        RUN
110          allocate ( sin_sz(kcolp) )                               !     sin(Sun Zenith.Dist)                          RUN
111          allocate ( sin_ho(kcolp) )                               !     sin(Sun Hour  Angle)                          RUN
112      ENDIF
113!                                                                  !
114!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
115
116
117
118
119! Insolation at the Top of the Atmosphere (TIME       PARAMETERS)
120! ===============================================================
121
122! Solar declination : delta
123! -------------------------
124
125      nj      = Day_TU+njYear(Mon_TU)                                   ! Nb of Days since Begin of Year           [day]
126
127      dlamm   = xlam  +      (nj-80) * step                             ! mean long. sun for ma-ja
128        anm   = dlamm -       xl
129       ranm   =   anm *       Dg2Rad
130       ranv   =  ranm +(2.0*ecc-xe3/4.0)*sin(ranm)                     &! anomalie  vraie                       [radian]
131     &             +5.0/4.0*ecc*ecc     *sin(2.0*ranm)                 &!
132     &            +13.0/12.0   *xe3     *sin(3.0*ranm)                  !
133      anv     =  ranv /       Dg2Rad                                    ! anomalie  vraie                       [degree]
134      tls     =   anv +       xl                                        ! longitude vraie                       [degree]
135      rlam    =   tls *       Dg2Rad                                    ! longitude vraie                       [radian]
136
137      sd      =    so *   sin(rlam)                                     ! sinus of solar declination angle (delta)   [-]
138      cd      =          sqrt(un_1 - sd * sd)                           ! cosin of solar declination angle (delta)   [-]
139                                                                        ! sinus delta = sin(obl)*sin(lambda)
140                                                                        !                       with lambda = real longitude
141                                                                        !(Phd.thesis Marie-France Loutre, ASTR-UCL, Belgium, 1993)
142      deltar  =          atan(       sd / cd)                           !
143      delta   = deltar/Dg2Rad                                           ! Solar Declination (degrees, angle sun at equator)
144
145
146! Eccentricity Effect
147! -------------------
148
149      dST_UA  =(1.0-ecc*ecc)/(1.0+ecc*cos(ranv))
150      ddt     = 1.0/dST_UA                                              ! 1 / normalized earth's sun distance
151
152
153! Insolation normal to the atmosphere (W/m2)
154! ------------------------------------------
155
156      rsunS0  = ddt *ddt *1368.
157
158
159! Time Equation (Should maybe be modified in case other than present
160! -------------  conditions are used, minor impact)
161
162      arg = om*nj
163      c1  = cos(     arg)
164      c2  = cos(2.d0*arg)
165      c3  = cos(3.d0*arg)
166      s1  = sin(     arg)
167      s2  = sin(2.d0*arg)
168      s3  = sin(3.d0*arg)
169
170      et=0.0072d0*c1 -0.0528d0*c2 -0.0012d0*c3                         &! difference between true solar and mean solar hour angles
171     &  -0.1229d0*s1 -0.1565d0*s2 -0.0041d0*s3                          ! (connected to the earth orbital rotation speed)
172
173
174
175
176! Insolation at the Top of the Troposphere (Auxiliary Variables)
177! ==============================================================
178
179
180! Day Length, Time Sunrise and Sunset at Sounding Grid Point (i_x0,j_y0)
181! ----------------------------------------------------------------------
182
183      IF (LDayS0)                                                   THEN
184          i = i_x0
185          j = j_y0
186
187               argc   = -tan(lat__r(ikl0))*tan(deltar)
188       if (abs(argc).gt. 1.d0)                                      then
189           ahc        =  0.d0
190        if    (argc .gt. 1.d0)                                      then
191           lhc       =  -1
192           LenDay    =   0.d0                                           ! Polar  Night
193        else
194           lhc       =   1
195           LenDay    =  24.d0                                           ! Midnight Sun
196        end if
197           SunRis    =   0.d0
198           SunSet    =   0.d0
199       else
200           ahc       =  acos(argc)
201           lhc       =   0
202
203        if(ahc.lt.0.d0) ahc  =  -ahc
204           ahc       =           ahc * 12.d0 / piNmbr
205           LenDay    =           ahc *  2.d0
206           SunRis    =  12.d0  - ahc - et
207           SunSet    =  SunRis + LenDay
208       end if
209
210      END IF
211
212! Time Angle
213! ----------
214
215         timh  = HourTU + MinuTU      / 60.d0                           ! time         [hour]
216
217      DO ikl=1,kcolp
218         ahor  = timh   + lon__h(ikl) - 12.d0 - et                      ! time angle   [hour]
219         ahorr = ahor   * piNmbr      / 12.d0                           ! time angle [radian]
220         chor  = cos(ahorr)
221
222
223
224
225! Solar Zenithal Distance zenitr (radians) and
226! Insolation (W/m2) at the Atmosphere Top  ===
227! =======================================
228
229         csz0S0(ikl) =      sinLat(ikl) *sd                             &
230     &                    + cosLat(ikl) *cd *chor
231! Martin modif to avoid cos(sza)=0 for LMDZ:
232!          csz0SV(ikl) = max(1E-6,csz0S0(ikl))
233         csz0S0(ikl) =  max(csz0S0(ikl) ,zer0)
234         csz_S0(ikl) =      csz0S0(ikl)
235
236! Martin control
237!        PRINT*,'csz0S0(',ikl,')=',csz0S0(ikl)
238! Martin control
239
240      ENDDO
241
242
243
244! Sun   Azimuth
245! =============
246
247      IF   (FaceS0)                                                 THEN
248        DO ikl = 1,kcolp
249
250
251! Slope Impact
252! ------------
253
254                  zenitr      = acos(csz0S0(ikl))                       !  Solar Zenithal Distance [radians]
255                  sin_sz(ikl) =  sin(zenitr)
256                  sin_ho(ikl) =  sin(ahorr)
257                  sin_sz(ikl) =  max(eps6      ,sin_sz(ikl))
258
259                  comes       =(sd-sinLat(ikl) *csz0S0(ikl))           &!  Cosine of Sun Azimuth
260     &                        /(   cosLat(ikl) *sin_sz(ikl))            !
261                  somes       =(cd*sin_ho(ikl)) /sin_sz(ikl)            !    Sine of Sun Azimuth
262
263          IF (abs(comes).gt.zer0)                                   THEN
264                  omesun = atan(somes/comes)
265            IF   (comes .lt.zer0)                                      &
266     &            omesun =  omesun + piNmbr
267
268            IF   (omesun.gt.         piNmbr)                           &
269     &            omesun =  -2.0d0 * piNmbr + omesun
270            IF   (omesun.lt.        -piNmbr)                           &
271     &            omesun =   2.0d0 * piNmbr + omesun
272
273          ELSE
274            IF   (somes .gt.zer0)                                   THEN
275                  omesun =   0.5d0 * piNmbr
276            ELSE
277                  omesun =   1.5d0 * piNmbr
278            END IF
279          END IF
280
281          IF (ikl.eq.ikl0)                                             &
282      &           omeswr =  omesun / Dg2Rad
283                  omesun =  -2.0d0 * piNmbr + omesun + DirAxX * Dg2Rad  !  Sun Azimuth (in MAR Reference Frame)
284                                                                        !           (positive counterclockwise)
285
286! Minimum Zenithal Distance
287! ~~~~~~~~~~~~~~~~~~~~~~~~~
288                  cszMIN(ikl) = 0.0d00
289          IF       (MMskS0)                                         THEN
290                                r_azim = omesun  / d_azim
291                                k_azim =       int(r_azim)
292
293! Mountain S0 MASKING
294! ~~~~~~~~~~~~~~~~~~~
295           IF(k_azim.le.     0)                                     THEN
296                                r_azim = r_azim  + n_azim
297                                k_azim = k_azim  + n_azim
298           END IF
299                                knazim = k_azim  + 1
300           IF(knazim.gt.n_azim) knazim = knazim  - n_azim
301
302                  cszMIN(ikl) = cszkS0(ikl,k_azim)+(r_azim  - k_azim)  &
303     &                        *(cszkS0(ikl,knazim)-cszkS0(ikl,k_azim))
304          END IF ! (MMskS0)
305
306! Cosine of Solar Normal Angle
307! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
308                  csz_S0(ikl) = slopS0(ikl)            *csz0S0(ikl)    &
309     &          + sin_sz(ikl) * slopS0(ikl)*sqrt(un_1  -slopS0(ikl))   &
310     &                        *              cos(omesun-omenS0(ikl))
311                  csz_S0(ikl) =              max(zer0  ,csz_S0(ikl))
312          IF     (csz0S0(ikl).le.cszMIN(ikl))                          &
313     &            csz_S0(ikl) = 0.0d00
314        END DO
315      END IF     ! (FaceS0)
316
317
318
319
320! Output Insolation Characteristics
321! =================================
322
323      IF (MinuTU .eq. 0 .and. Sec_TU .eq. 0)                           THEN
324
325         ahor  = timh   + lon__h(ikl0)  - 12.d0 - et
326         zenitr=     acos(csz0S0(ikl0)) / Dg2Rad
327! #AZ    anormr=     acos(csz_S0(ikl0)) / Dg2Rad
328! #AZ    omenwr= DirAxX - omenS0(ikl0)  / Dg2Rad
329! #AZ    if (omenwr.lt.  0.) omenwr = omenwr + 360.d0
330! #AZ    if (omenwr.gt.360.) omenwr = omenwr - 360.d0
331! #AZ                        omeswr = 360.d0 - omeswr
332! #AZ    if (omeswr.lt.  0.) omeswr = omeswr + 360.d0
333! #AZ    if (omeswr.gt.360.) omeswr = omeswr - 360.d0
334
335        write(4,1)           lat__r(ikl0)/Dg2Rad,lon__h(ikl0)           &
336     &            ,Day_TU,Mon_TU,HourTU,MinuTU,Sec_TU
337 1      format(/,' lat.=',f6.1,3x,'long.=',f7.1,4x,'date :',i3,'-',i2,  &
338     &           ' / ',i2,' h.UT',i3,' min.',i3,' sec.')
339        write(4,2) i_x0,j_y0,lat__r(ikl0)/Dg2Rad,lon__h(ikl0)
340 2      format(' Sounding at (',i3,i3,') / (',f6.2,'dg,',f6.2,'ho)')
341        write(4,3) rsunS0*csz_S0(ikl0)     ,ahor,zenitr                 &
342! #AZ&            ,omeswr,omenwr           ,     anormr                 &
343     &            ,delta
344 3      format(' Insolation [W/m2]  = ',f7.2,'   Hor.Angle = ',f7.2,    &
345     &         '   Zenith.Angle = ',f7.2                                &
346! #AZ&      ,/,'                      ', 7x ,'   Sol.Azim. = ',f7.2     &
347! #AZ&      ,/,'                      ', 7x ,'   Nrm.Azim. = ',f7.2     &
348! #AZ&        ,'   Normal Angle = ',f7.2                                &
349     &      ,/,' Solar Declination  = ',f7.2)
350
351      IF       (LDayS0)                                             THEN
352       if (lhc.eq.-1)                                                   &
353     &  write(4,4) SunRis,LenDay,SunSet
354 4      format(' Sun Rise Time [h]  = ',f7.2,'   Day Leng. = ',f7.2,    &
355     &         '   Sun Set Time = ',f7.2,'  -- POLAR  NIGHT --')
356       if (lhc.eq. 0)                                                   &
357     &  write(4,5) SunRis,LenDay,SunSet
358 5      format(' Sun Rise Time [h]  = ',f7.2,'   Day Leng. = ',f7.2,    &
359     &         '   Sun Set Time = ',f7.2,'  -- SOLAR  TIME  --')
360       if (lhc.eq. 1)                                                   &
361     &  write(4,6) SunRis,LenDay,SunSet
362 6      format(' Sun Rise Time [h]  = ',f7.2,'   Day Leng. = ',f7.2,    &
363     &         '   Sun Set Time = ',f7.2,'  -- MIDNIGHT SUN --')
364      END IF ! (LDayS0)
365
366      END IF !
367
368
369
370
371!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
372!                                                                  !
373! DE-ALLOCATION
374! =============
375
376      IF (FlagDALLOC)                                         THEN !
377          deallocate ( cszMIN )                                    !  Sum of cos(zenith.Dist) on Azimut Intervals      RUN
378          deallocate ( sin_sz )                                    !     sin(Sun zenith.Dist)                          RUN
379          deallocate ( sin_ho )                                    !     sin(Sun Hour  Angle)                          RUN
380      ENDIF
381!                                                                  !
382!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
383
384
385
386
387      return
388      end subroutine PHY_Atm_S0_RUN
Note: See TracBrowser for help on using the repository browser.