source: branches/publications/ORCHIDEE_IPSLCM5A2.1.r5307/ORCHIDEE/src_global/solar.f90 @ 7683

Last change on this file since 7683 was 3851, checked in by josefine.ghattas, 8 years ago

For assimilation :

  • intersurf.f90, teststomate.f90 : renamed l_first_intersurf into lstep_init_intersurf
  • weather.f90, solar.f90, hydrolc : renamed firstcall_weater_xxx into lstep_init_xxx
  • removed first_call from constantes_soil.f90, constantes.f90, pft_parameters.f90 because these subroutines are only called onces, intented full module
  • Property svn:keywords set to Revision Date HeadURL Date Author Revision
File size: 20.8 KB
Line 
1! =================================================================================================================================
2! MODULE       : solar
3!
4! CONTACT      : orchidee-help _at_ ipsl.jussieu.fr
5!
6! LICENCE      : IPSL (2011)
7! This software is governed by the CeCILL licence see ORCHIDEE/ORCHIDEE_CeCILL.LIC
8!
9!>\BRIEF        This module define solar
10!!
11!!\n DESCRIPTION:
12!!
13!! RECENT CHANGE(S): None
14!!
15!! REFERENCE(S) : None
16!!
17!! SVN          :
18!! $HeadURL$
19!! $Date$
20!! $Revision$
21!! \n
22!_ ================================================================================================================================
23
24MODULE solar
25
26  USE defprec
27  USE constantes
28  USE ioipsl_para 
29
30  IMPLICIT NONE
31
32
33
34CONTAINS
35
36
37!! ================================================================================================================================
38!! SUBROUTINE   : solarang
39!!
40!>\BRIEF         This subroutine computes the solar angle according to the method used by GSWP and developed by J.C. Morill.
41!!
42!! DESCRIPTION  : None
43!!
44!! RECENT CHANGE(S): None
45!!
46!! MAIN OUTPUT VARIABLE(S): ::csang
47!!
48!! REFERENCE(S) : See for further details :
49!! http://www.atmo.arizona.edu/~morrill/swrad.html
50!!
51!! FLOWCHART    : None
52!! \n
53!_ ================================================================================================================================
54
55  SUBROUTINE solarang (julian, julian0, iim, jjm, lon, lat, csang)
56
57    USE calendar
58
59    IMPLICIT NONE
60
61    !! 0. Parameters and variables declaration
62   
63    !! 0.1 Input variables
64
65    INTEGER, INTENT(in)                      :: iim, jjm        !!
66    REAL, INTENT(in)                         :: julian          !!
67    REAL, INTENT(in)                         :: julian0         !!
68    REAL, DIMENSION(iim,jjm), INTENT(in)     :: lon, lat        !!
69
70    !! 0.2 Output variables
71
72    REAL, DIMENSION(iim,jjm), INTENT(out)    :: csang           !!
73
74    !! 0.4 Local variables
75
76    REAL                                     :: gamma           !!
77    REAL                                     :: dec             !!
78    REAL                                     :: decd            !!
79    REAL                                     :: et              !!
80    REAL                                     :: gmt             !!
81    REAL                                     :: le              !!
82    REAL                                     :: ls              !!
83    REAL                                     :: lcorr           !!
84    REAL                                     :: latime          !!
85    REAL                                     :: omega           !!
86    REAL                                     :: omegad          !!
87    REAL                                     :: llatd, llat     !!
88    INTEGER                                  :: igmt            !!
89    INTEGER                                  :: ilon, ilat      !!
90    INTEGER, SAVE, ALLOCATABLE, DIMENSION(:) :: zone            !!
91!$OMP THREADPRIVATE(zone)
92    REAL, SAVE, ALLOCATABLE, DIMENSION(:)    :: lhour           !!
93!$OMP THREADPRIVATE(lhour)
94    !
95    LOGICAL                                  :: check = .FALSE. !!
96!$OMP THREADPRIVATE(check)
97
98!_ ================================================================================================================================
99
100    IF (check) WRITE(numout,*) 'Calendar information', julian, julian0
101    !-
102    !- 1) Day angle gamma
103    !-
104    !   gamma = 2.*pi*MOD(julian,one_year)/one_year
105    gamma = 2.*pi*(julian-julian0)/one_year
106    !-
107    !- 2) Solar declination (assumed constant for a 24 hour period)  page 7
108    !-    in radians
109    !-
110    IF (check) WRITE(numout,*) 'Solar declination', gamma, one_year
111    !
112    dec = ( 0.006918-0.399912*COS(gamma)+0.070257*SIN(gamma) &
113         &       -0.006758*COS(2*gamma)+0.000907*SIN(2*gamma)      &
114         &       -0.002697*COS(3*gamma)+0.00148*SIN(3*gamma))
115    decd = dec*(180/pi)
116    !-
117    !- maximum error 0.0006 rad (<3'), leads to error
118    !- of less than 1/2 degree in zenith angle
119    !-
120    IF (check) WRITE(numout,*) 'Equation of time', decd
121    !- 3)  Equation of time  page 11
122    !-
123    et = ( 0.000075+0.001868*COS(gamma)-0.032077*SIN(gamma)&
124         &      -0.014615*COS(2*gamma)-0.04089*SIN(2*gamma))*229.18
125    !-
126    !- Get the time zones for the current time
127    !-
128    gmt = 24.*(julian-INT(julian))
129    IF (.NOT.ALLOCATED(zone))  ALLOCATE(zone(iim))
130    IF (.NOT.ALLOCATED(lhour)) ALLOCATE(lhour(iim))
131    !-
132    !igmt = NINT(gmt)
133    IF (check) WRITE(numout,*) 'Get time zone', et, gmt
134    CALL time_zone(gmt, iim, jjm, lon, zone, lhour)
135    !-
136    !- Start a loop through the grid
137    !-
138    IF (check) WRITE(numout,*) 'Start a loop through the grid'
139    DO ilon=1,iim
140       !---
141       !--- 4) Local apparent time  page 13
142       !---
143       !--- ls     standard longitude (nearest 15 degree meridian)
144       !--- le     local longtitude
145       !--- lhour  local standard time
146       !--- latime local apparent time
147       !--- lcorr  longitudunal correction (minutes)
148       !---
149       le = lon(ilon,1)
150       ls = ((zone(ilon)-1)*15)-180.
151       lcorr = 4.*(ls-le)*(-1)
152       latime = lhour(ilon)+lcorr/60.+et/60.
153       IF (latime <  0.) latime = latime+24
154       IF (latime > 24.) latime = latime-24
155       !---
156       !--- 5) Hour angle omega  page 15
157       !---
158       !--- hour angle is zero at noon, positive in the morning
159       !--- It ranges from 180 to -180
160       !---
161       !--- omegad is omega in degrees, omega is in radians
162       !---
163       omegad = (latime-12.)*(-15.)
164       omega  = omegad*pi/180.
165       !---
166       DO ilat=1,jjm
167          !-----
168          !----- 6)  Zenith angle  page 15
169          !-----
170          !----- csang cosine of zenith angle (radians)
171          !----- llatd =  local latitude in degrees
172          !----- llat  =  local latitude in radians
173          !-----
174          llatd = lat(1,ilat)
175          llat  = llatd*pi/180.
176          csang(ilon,ilat) = &
177               &  MAX(zero,SIN(dec)*SIN(llat)+COS(dec)*COS(llat)*COS(omega))
178       ENDDO
179    ENDDO
180    !----------------------
181  END SUBROUTINE solarang
182
183
184!! ================================================================================================================================
185!! SUBROUTINE   : time_zone
186!!
187!>\BRIEF         
188!!
189!! DESCRIPTION  : None
190!!
191!! RECENT CHANGE(S): None
192!!
193!! MAIN OUTPUT VARIABLE(S): ::zone, ::lhour
194!!
195!! REFERENCE(S) : None
196!!
197!! FLOWCHART    : None
198!! \n
199!_ ================================================================================================================================
200
201  SUBROUTINE time_zone (gmt, iim, jjm, lon, zone, lhour)
202
203    IMPLICIT NONE
204
205    !! 0. Parameters and variables declaration
206   
207    !! 0.1 Input variables
208
209    INTEGER, INTENT(in)                   :: iim, jjm        !!
210    REAL, DIMENSION(iim,jjm), INTENT(in)  :: lon             !!
211    REAL, INTENT(in)                      :: gmt             !!
212
213    !! 0.2 Output variables
214
215    INTEGER, DIMENSION(iim), INTENT(out)  :: zone            !!
216    REAL, DIMENSION(iim), INTENT(out)     :: lhour           !!
217
218    !! 0.4 Local variables
219
220    INTEGER                               :: deg             !!
221    !!??   REAL :: deg
222    INTEGER                               :: i, ilon, change !! Indices
223
224!_ ================================================================================================================================
225
226    DO ilon=1,iim
227       !---
228       !--- Convert longitude index to longtitude in degrees
229       !---
230       deg = lon(ilon,1)
231       !---
232       !--- Determine into which time zone (15 degree interval) the
233       !--- longitude falls.
234       !---
235       DO i=1,25
236          IF (deg < (-187.5+(15*i))) THEN
237             zone(ilon) = i
238             IF (zone(ilon) == 25)   zone(ilon) = 1
239             EXIT
240          ENDIF
241       ENDDO
242       !---
243       !--- Calculate change (in number of hours) from GMT time to
244       !--- local hour.  Change will be negative for zones < 13 and
245       !--- positive for zones > 13.
246       !---
247       !--- There is also a correction for lhour < 0 and lhour > 23
248       !--- to lhour between 0 and 23.
249       !---
250       IF (zone(ilon) < 13) THEN
251          change = zone(ilon)-13
252          lhour(ilon) = gmt+change
253       ENDIF
254       !---
255       IF (zone(ilon) == 13) THEN
256          lhour(ilon) = gmt
257       ENDIF
258       !---
259       IF (zone(ilon) > 13) THEN
260          change = zone(ilon)-13
261          lhour(ilon) = gmt+change
262       ENDIF
263       IF (lhour(ilon) <  0) lhour(ilon) = lhour(ilon)+24
264       IF (lhour(ilon) >= 24) lhour(ilon) = lhour(ilon)-24
265       !---
266    ENDDO
267    !-----------------------
268  END SUBROUTINE time_zone
269
270
271
272!! ================================================================================================================================
273!! SUBROUTINE   : downward_solar_flux
274!!
275!>\BRIEF         
276!!
277!! DESCRIPTION  : None
278!!
279!! RECENT CHANGE(S): None
280!!
281!! MAIN OUTPUT VARIABLE(S): ::solad, ::solai
282!!
283!! REFERENCE(S) :See for further details :
284!! http://www.atmo.arizona.edu/~morrill/swrad.html
285!!
286!! FLOWCHART    : None
287!! \n
288!_ ================================================================================================================================
289
290  SUBROUTINE downward_solar_flux (npoi,latitude,calendar_str,jday,rtime,cloud,nband,solad,solai)
291
292    IMPLICIT NONE
293
294    !! 0. Parameter and variables declaration
295   
296    !! 0.1 Input variables
297   
298    INTEGER, INTENT(in)                      :: npoi               !! number of grid points (unitless)
299    REAL, DIMENSION(npoi), INTENT(in)        :: latitude           !! latitude in degrees
300    CHARACTER(LEN=20),INTENT(in)             :: calendar_str       !! Calendar type
301    REAL, INTENT(in)                         :: jday
302    REAL,INTENT(in)                          :: rtime
303    REAL, DIMENSION(npoi), INTENT(in)        :: cloud              !! cloud fraction (0-1, unitless)
304    INTEGER,INTENT(in)                       :: nband              !! number of visible bands (unitless)
305
306    !! 0.2 Output variables
307
308    REAL, DIMENSION(npoi,nband), INTENT(out) :: solad              !! direct downward solar flux
309                                                                   !! @tex $(W.m^{-2})$ @endtex
310    REAL, DIMENSION(npoi,nband), INTENT(out) :: solai              !! diffuse downward solar flux
311                                                                   !! @tex $(W.m^{-2})$ @endtex
312
313    !! 0.4 Local variables
314
315    ! Parametres orbitaux:
316
317    REAL, SAVE                               :: ecc                !! Eccentricity
318!$OMP THREADPRIVATE(ecc)
319    REAL, SAVE                               :: perh               !! Longitude of perihelie
320!$OMP THREADPRIVATE(perh)
321    REAL, SAVE                               :: xob                !! obliquity
322!$OMP THREADPRIVATE(xob)
323    REAL, PARAMETER                          :: pir = pi/180.      !!
324    REAL, SAVE                               :: step               !!
325!$OMP THREADPRIVATE(step)
326    REAL                                     :: xl                 !!
327    REAL                                     :: so                 !!
328    REAL                                     :: xllp               !!
329    REAL                                     :: xee                !!
330    REAL                                     :: xse                !! 
331    REAL                                     :: xlam               !!
332    REAL                                     :: dlamm              !!
333    REAL                                     :: anm                !!
334    REAL                                     :: ranm               !!
335    REAL                                     :: ranv               !!
336    REAL                                     :: anv                !!
337    REAL                                     :: tls                !!
338    REAL                                     :: rlam               !!
339    REAL                                     :: sd                 !!
340    REAL                                     :: cd                 !!
341    REAL                                     :: deltar             !!
342    REAL                                     :: delta              !!
343    REAL                                     :: Dis_ST             !! 
344    REAL                                     :: ddt                !!
345    REAL                                     :: orbit              !!
346    REAL                                     :: angle              !!
347    REAL                                     :: xdecl              !!
348    REAL                                     :: xlat               !!
349    REAL, DIMENSION(npoi)                    :: trans              !!
350    REAL, DIMENSION(npoi)                    :: fdiffuse           !!
351    REAL, DIMENSION(npoi)                    :: coszen             !! cosine of solar zenith angle
352    REAL                                     :: sw                 !!
353    REAL                                     :: frac               !!
354    INTEGER                                  :: i,ib               !!
355    LOGICAL, SAVE                            :: lstep_init = .TRUE. !!
356!$OMP THREADPRIVATE(lstep_init)
357
358!_ ================================================================================================================================
359   
360    IF (lstep_init) THEN
361       IF ( TRIM(calendar_str) .EQ. 'gregorian' ) THEN 
362          step = 1.0
363       ELSE
364          step = one_year/365.2425
365       ENDIF
366       !-
367       ! Read Orbital Parameters
368       !-
369
370       !Config Key   = ECCENTRICITY
371       !Config Desc  = Use prescribed values
372       !Config If    = ALLOW_WEATHERGEN
373       !Config Def   = 0.016724
374       !Config Help  =
375       !Config Units = [-]
376       ecc = 0.016724
377       CALL getin_p ('ECCENTRICITY',ecc)
378       WRITE(numout,*) 'ECCENTRICITY: ',ecc
379       !
380       !Config Key  = PERIHELIE
381       !Config Desc = Use prescribed values
382       !Config If   = ALLOW_WEATHERGEN
383       !Config Def  = 102.04
384       !Config Help  =
385       !Config Units = [-]
386       perh = 102.04
387       CALL getin_p ('PERIHELIE',perh)
388       WRITE(numout,*) 'PERIHELIE: ',perh
389       !
390       !Config Key  = OBLIQUITY
391       !Config Desc = Use prescribed values
392       !Config If   = ALLOW_WEATHERGEN
393       !Config Def  = 23.446
394       !Config Help  =
395       !Config Units = [Degrees]
396       xob = 23.446
397       CALL getin_p ('OBLIQUITY',xob)
398       WRITE(numout,*) 'OBLIQUITY: ',xob
399
400       lstep_init = .FALSE.
401    ENDIF
402
403    !-
404    ! calendar and orbital calculations
405    !-
406
407    !-
408    ! calculate the earth's orbital angle (around the sun) in radians
409    orbit = 2.0*pi*jday/365.2425
410    !-
411    ! calculate the solar hour angle in radians
412    angle = 2.0*pi*(rtime-12.0)/24.0
413    !-
414    ! calculate the current solar declination angle
415    ! ref: global physical climatology, hartmann, appendix a
416    !
417    ! xdecl = 0.006918 &
418    !        -0.399912*cos(orbit) &
419    !        +0.070257*sin(orbit) &
420    !        -0.006758*cos(2.0*orbit) &
421    !        +0.000907*sin(2.0*orbit) &
422    !        -0.002697*cos(3.0*orbit) &
423    !        +0.001480*sin(3.0*orbit)
424    !
425    ! calculate the effective solar constant,
426    ! including effects of eccentricity
427    ! ref: global physical climatology, hartmann, appendix a
428    !
429    ! sw = 1370.*( 1.000110 &
430    !             +0.034221*cos(orbit) &
431    !             +0.001280*sin(orbit) &
432    !             +0.000719*cos(2.0*orbit) &
433    !             +0.000077*sin(2.0*orbit))
434    !
435    ! correction Marie-France Loutre
436    !
437    !    orbital parameters
438    !
439    !    ecc = 0.016724
440    !    perh = 102.04
441    !    xob = 23.446
442    !-
443    xl = perh+180.0
444    ! so : sinus of obliquity
445    so = sin(xob*pir)
446    !-
447    xllp = xl*pir
448    xee  = ecc*ecc
449    xse  = sqrt(1.0d0-xee)
450    ! xlam : true long. sun for mean long. = 0
451    xlam = (ecc/2.0+ecc*xee/8.0d0)*(1.0+xse)*sin(xllp)-xee/4.0 &
452         &      *(0.5+xse)*sin(2.0*xllp)+ecc*xee/8.0*(1.0/3.0+xse) &
453         &      *sin(3.0*xllp)
454    xlam  =2.0d0*xlam/pir
455    ! dlamm : mean long. sun for ma-ja
456    dlamm =xlam+(INT(jday)-79)*step
457    anm   = dlamm-xl
458    ranm  = anm*pir
459    xee   = xee*ecc
460    ! ranv : anomalie vraie    (radian)
461    ranv = ranm+(2.0*ecc-xee/4.0)*sin(ranm)+5.0/4.0*ecc*ecc &
462         &      *sin(2.0*ranm)+13.0/12.0*xee*sin(3.0*ranm)
463    ! anv  : anomalie vraie    (degrees)
464    anv  = ranv/pir
465    ! tls  : longitude vraie   (degrees)
466    tls  = anv+xl
467    ! rlam : longitude vraie   (radian)
468    rlam = tls*pir
469    ! sd and cd: cosinus and sinus of solar declination angle (delta)
470    ! sinus delta = sin (obl)*sin(lambda) with lambda = real longitude
471    ! (Phd. thesis of Marie-France Loutre, ASTR-UCL, Belgium, 1993)
472    sd    = so*sin(rlam)
473    cd    = sqrt(1.0d0-sd*sd)
474    ! delta : Solar Declination (degrees, angle sun at equator)
475    deltar = atan(sd/cd)
476    delta  = deltar/pir
477    !-
478    ! Eccentricity Effect
479    !-
480    Dis_ST  =(1.0-ecc*ecc)/(1.0+ecc*cos(ranv))
481    ! ddt :    1 / normalized earth's sun distance
482    ddt = 1.0/Dis_ST
483    !-
484    ! Insolation normal to the atmosphere (W/m2)
485    !-
486    sw  = ddt *ddt *1370.d0
487    !-
488    xdecl = deltar
489    !-
490    ! solar calculations
491    !-
492    do i=1,npoi
493       !---
494       !-- calculate the latitude in radians
495       !---
496       xlat = latitude(i)*pi/180.0
497       !---
498       !-- calculate the cosine of the solar zenith angle
499       !---
500       coszen(i) = MAX(0.0, (sin(xlat)*sin(xdecl) &
501            &                      + cos(xlat)*cos(xdecl)*cos(angle)))
502       !---
503       !-- calculate the solar transmission through the atmosphere
504       !-- using simple linear function of tranmission and cloud cover
505       !---
506       !-- note that the 'cloud cover' data is typically obtained from
507       !-- sunshine hours -- not direct cloud observations
508       !---
509       !-- where, cloud cover = 1 - sunshine fraction
510       !---
511       !-- different authors present different values for the slope and
512       !-- intercept terms of this equation
513       !---
514       !-- Friend, A: Parameterization of a global daily weather generator
515       !-- for terrestrial ecosystem and biogeochemical modelling,
516       !-- Ecological Modelling
517       !---
518       !-- Spitters et al., 1986: Separating the diffuse and direct component
519       !-- of global radiation and its implications for modeling canopy
520       !-- photosynthesis, Part I: Components of incoming radiation,
521       !-- Agricultural and Forest Meteorology, 38, 217-229.
522       !---
523       !-- A. Friend       : trans = 0.251+0.509*(1.0-cloud(i))
524       !-- Spitters et al. : trans = 0.200+0.560*(1.0-cloud(i))
525       !---
526       !-- we are using the values from A. Friend
527       !---
528       trans(i) = 0.251+0.509*(1.0-cloud(i))
529       !---
530       !-- calculate the fraction of indirect (diffuse) solar radiation
531       !-- based upon the cloud cover
532       !---
533       !-- note that these relationships typically are measured for either
534       !-- monthly or daily timescales, and may not be exactly appropriate
535       !-- for hourly calculations -- however, in ibis, cloud cover is fixed
536       !-- through the entire day so it may not make much difference
537       !---
538       !-- method i --
539       !---
540       !-- we use a simple empirical relationships
541       !-- from Nikolov and Zeller (1992)
542       !---
543       !-- Nikolov, N. and K.F. Zeller, 1992:  A solar radiation algorithm
544       !-- for ecosystem dynamics models, Ecological Modelling, 61, 149-168.
545       !---
546       fdiffuse(i) = 1.0045+trans(i)*( 0.0435+trans(i) &
547            &                                 *(-3.5227+trans(i)*2.6313))
548       !---
549       IF (trans(i) > 0.75) fdiffuse(i) = 0.166
550       !---
551       !-- method ii --
552       !---
553       !-- another method was suggested by Spitters et al. (1986) based on
554       !-- long-term data from the Netherlands
555       !--
556       !-- Spitters et al., 1986: Separating the diffuse and direct component
557       !-- of global radiation and its implications for modeling canopy
558       !-- photosynthesis, Part I: Components of incoming radiation,
559       !-- Agricultural and Forest Meteorology, 38, 217-229.
560       !---
561       !--       if ((trans == 0.00).and.(trans < 0.07)) then
562       !--         fdiffuse = 1.0
563       !--       else if ((trans >= 0.07).and.(trans < 0.35)) then
564       !--         fdiffuse = 1.0-2.3*(trans-0.07)**2
565       !--       else if ((trans >= 0.35).and.(trans < 0.75)) then
566       !--         fdiffuse = 1.33-1.46*trans
567       !--       ELSE
568       !--         fdiffuse = 0.23
569       !--       endif
570       !---
571    ENDDO
572    !-----------------------
573    !-
574    ! do for each waveband
575    !-
576    DO ib=1,nband
577       !---
578       !-- calculate the fraction in each waveband
579       !---
580       !-- GK010200
581       IF (nband == 2) then
582          frac = 0.46+0.08*REAL(ib-1)
583       ELSE
584          frac = 1./REAL(nband)
585       ENDIF
586       !---
587       DO i=1,npoi
588          !-----
589          !---- calculate the direct and indirect solar radiation
590          !-----
591          solad(i,ib) = sw*coszen(i)*frac*trans(i)*(1.-fdiffuse(i))
592          solai(i,ib) = sw*coszen(i)*frac*trans(i)*fdiffuse(i)
593       ENDDO
594    ENDDO
595
596
597  END SUBROUTINE downward_solar_flux
598
599 
600END MODULE solar
Note: See TracBrowser for help on using the repository browser.