source: branches/publications/ORCHIDEE_GLUC_r6545/src_stomate/grassland_fonctions.f90 @ 6737

Last change on this file since 6737 was 3771, checked in by albert.jornet, 8 years ago

Merge: new Grassland Management Module (GRM) from revision [3759/perso/jinfeng.chang/ORCHIDEE-MICT/ORCHIDEE/ORCHIDEE]. Done by Jinfeng.
File size: 16.6 KB
Line 
1! =================================================================================================================================
2! MODULE       : grassland_fonctions
3!
4! CONTACT      : orchidee-help _at_ ipsl.jussieu.fr
5!
6! LICENCE      : IPSL (2006)
7! This software is governed by the CeCILL licence see
8! ORCHIDEE/ORCHIDEE_CeCILL.LIC
9!
10!>\BRIEF       This module defined basic aggregation functions used
11!! in grassland management module
12!!
13!!\n DESCRIPTION : None
14!!
15!! RECENT CHANGE(S) : None
16!!
17!! REFERENCE(S) : None
18!!
19!! \n
20!_
21!================================================================================================================================
22MODULE grassland_fonctions
23
24  USE grassland_constantes
25  USE constantes
26
27  ! Used functions
28  ! Euler_funct
29  ! linreg_pasim
30
31  IMPLICIT NONE
32
33  PUBLIC :: Euler_funct
34
35  INTEGER(i_std), SAVE :: emplacement 
36  ! permet de se souvenir du dernier emplacement dans l'interpolation
37
38CONTAINS
39
40  !!!!!!!!!!!!!!!!!!!!!!!!
41  !!!!!! EULER
42  !!!!!!!!!!!!!!!!!!!!!!!!
43
44  SUBROUTINE Euler_funct(npts, dt, dY, Y)
45    ! Y(n+1) = Y(n) + h*F(tn,Y(n))
46    ! d'ou :
47    ! Y = Yp + h*dY
48    ! ou h = dt
49
50    INTEGER(i_std), INTENT(in)                   :: npts   
51    REAL(R_std), INTENT(in)                      :: dt     
52    REAL(R_std), DIMENSION(npts), INTENT(in)     :: dY     
53    ! result
54    REAL(R_std), DIMENSION(npts), INTENT(inout)  :: Y
55
56    Y(:) = Y(:) + dt*dY(:)
57
58  END SUBROUTINE Euler_funct
59
60
61  SUBROUTINE Euler_X(npts, ncol, dt, dY, Y)
62    ! matrix
63    INTEGER(i_std), INTENT(in)   :: npts 
64    INTEGER(i_std), INTENT(in)   :: ncol ! number colonne of matrix
65    REAL(R_std), INTENT(in)      :: dt   
66    REAL(R_std), DIMENSION(npts, ncol) :: dY   
67    REAL(R_std), DIMENSION(npts, ncol) :: Y   
68
69    INTEGER(i_std) :: i
70
71    DO i=1, ncol
72      Y(:,i) = Y(:,i) + dt*dY(:,i)
73    END DO
74
75  END SUBROUTINE Euler_X
76
77  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
78  !!!!!!!!      FTSIGM
79  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
80
81  SUBROUTINE fTsigm(npts, Tform, &
82     T0form, T0pform, qftform, &
83     fTsigmform)
84
85    ! functiong of temperature equation 3.11
86    INTEGER(i_std), INTENT(in)                :: npts
87    REAL(R_std), DIMENSION(npts), INTENT(in)  :: Tform
88    REAL(R_std), INTENT(in)                   :: T0form
89    REAL(R_std), INTENT(in)                   :: T0pform
90    REAL(R_std), INTENT(in)                   :: qftform
91    REAL(r_std), DIMENSION(npts), INTENT(out) :: fTsigmform
92   
93    INTEGER(i_std) :: i
94
95    DO i=1, npts
96      IF ((Tform(i) .GT. T0form) .AND. (Tform(i) .LT. T0pform)) THEN
97          fTsigmform(i) = ((Tform(i) - T0form)**qftform)*(T0pform - Tform(i))/ &
98             (((293.0 - T0form)**qftform)*(T0pform - 293.0)) 
99      ELSE
100          fTsigmform(i) = 0.0
101      ENDIF
102    ENDDO
103
104  END SUBROUTINE fTsigm
105
106
107  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
108  !!!!!!!!!!!!!!!! INTERPOLATION EXTRAPOLATION
109  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
110
111  SUBROUTINE interpolation_extrapolation (npts, nbpoints, table, absc, nouv_annee, resultat)
112
113! Principe suivi : cette fonction n'est appelée que pour faire l'interpolation
114!sur les fichiers météo (ta, ea, pa, iatmtot, nh3a, co2, u).
115! Elle est appelée à chaque pas de temps. Il faut faire deux remarques :
116! 1)- pour tous les points de simulations l'interpolation se fait dans
117!le même intervalle de temps pour un fichier.
118! 2)- Tous les fichiers météo étant construit sur le même format,
119!ils font tous leur interpolation dans le même intervalle pour un pas de temps donné.
120! Nous allons donc conserver en mémoire le dernier emplacement
121!(dans quel intervalle) d'interpolation, en sachant que la suivante sera soit au même
122!endroit soit dans l'intervalle suivant.
123
124    INTEGER(i_std), INTENT(in) :: npts
125    INTEGER(i_std), INTENT(in) :: nbpoints 
126    REAL(r_std), DIMENSION(npts,2,nbpoints), INTENT(in) :: table
127    REAL(r_std), INTENT(in)  :: absc
128    LOGICAL, intent(in) :: nouv_annee
129    REAL(r_std), DIMENSION(npts), INTENT(out) :: resultat
130
131    INTEGER(i_std) :: i
132    INTEGER(i_std) :: j
133    LOGICAL, DIMENSION(npts) :: calcul
134
135    calcul(:) = .FALSE. 
136
137    ! remarquons que tous les tableaux appelant cette fonction sont au même niveau pour l'interpolation
138    ! remarquons que pour i=1,npts tous les points sont interpolés au même endroit
139
140    IF (nouv_annee) THEN   ! if it's the firts call during a new_year
141       
142      ! table has npts line, but all points are the make on the same model,
143      ! so we search on point 1 and it's the same on all other points
144      IF ((absc .GE. table(1,1,1)) .AND. (absc .LE. table(1,1,nbpoints))) THEN 
145
146        ! if table (1,2,i+1) = -999 it's the end of the array table
147        DO i=1,nbpoints-1       
148          IF ((absc .GE. table(1,1,i)) .AND. (absc .LE. table(1,1,i+1)) .AND.&
149            (table(1,2,i+1).NE. -999.0)) THEN
150               
151            calcul(:) = .TRUE. 
152            resultat(:) = table(:,2,i)*(table(:,1,i+1)-absc)/(table(:,1,i+1)-table(:,1,i)) + &
153                   table(:,2,i+1)*(absc - table(:,1,i))/(table(:,1,i+1)-table(:,1,i))
154            emplacement = i
155               
156          ELSE IF ((table(1,2,i+1).EQ. -999.0).AND. (calcul(1) .EQ. .FALSE.)) THEN
157            calcul(:) = .TRUE.
158            resultat(:) = table(:,2,i-1)*(table(:,1,i)-absc)/ &
159                   (table(:,1,i)-table(:,1,i-1)) + &
160                   table(:,2,i)*(absc - table(:,1,i-1))/&
161                   (table(:,1,i)-table(:,1,i-1))
162            emplacement = i
163
164          END IF
165        END DO
166      ELSE
167        IF (absc .LT. table(1,1,1)) THEN
168
169          calcul(:) = .TRUE.
170          resultat(:) = table(:,2,1)*(table(:,1,2)-absc)/(table(:,1,2)-table(:,1,1)) + &
171                 table(:,2,2)*(absc - table(:,1,1))/(table(:,1,2)-table(:,1,1))
172          emplacement = 0
173             
174        ELSE
175
176          calcul(:) = .TRUE.
177          resultat(:) = table(:,2,nbpoints-1)*(table(:,1,nbpoints)-absc)/ &
178                 (table(:,1,nbpoints)-table(:,1,nbpoints-1)) + &
179                 table(:,2,nbpoints)*(absc - table(:,1,nbpoints-1))/&
180                 (table(:,1,nbpoints)-table(:,1,nbpoints-1))
181          emplacement = nbpoints
182             
183        END IF
184      END IF
185    ELSE              ! si ce n'est pas le premier appel de l'année
186      ! dans ce cas nous connaissons le dernier emplacement utilisé. Nous repartons de là
187      IF ((emplacement .NE. 0) .AND. (emplacement .NE. nbpoints)) THEN
188        IF ((absc .GE. table(1,1,emplacement)) .AND. (absc .LE. table(1,1,emplacement+1)) .AND.&
189             (table(1,2,emplacement+1).NE. -999.0)) THEN
190              ! emplacement ne change pas
191
192          calcul(:) = .TRUE. 
193          resultat(:) = &
194                table(:,2,emplacement)*(table(:,1,emplacement+1)-absc)/&
195                (table(:,1,emplacement+1)-table(:,1,emplacement)) + &
196                 table(:,2,emplacement+1)*(absc - table(:,1,emplacement))/&
197                 (table(:,1,emplacement+1)-table(:,1,emplacement))
198             
199         
200        ELSE IF ((table(1,2,emplacement+1).EQ. -999.0).AND. (calcul(1) .EQ. .FALSE.)) THEN
201              !emplacement ne change pas
202          calcul(:) = .TRUE.
203          resultat(:) = table(:,2,emplacement-1)*(table(:,1,emplacement)-absc)/ &
204                 (table(:,1,emplacement)-table(:,1,emplacement-1)) + &
205                 table(:,2,emplacement)*(absc - table(:,1,emplacement-1))/&
206                 (table(:,1,emplacement)-table(:,1,emplacement-1))
207
208           resultat(:) = table(:,2,emplacement)
209             
210        ELSE IF (emplacement .NE. nbpoints-1) THEN
211          ! dans ce cas là l'interpolation est dans la case suivante
212          emplacement = emplacement + 1
213          IF ((absc .GE. table(1,1,emplacement)) .AND. (absc .LE. table(1,1,emplacement+1)) .AND.&
214               (table(1,2,emplacement+1).NE. -999.0)) THEN
215            ! emplacement ne change pas
216            calcul(:) = .TRUE. 
217            resultat(:) = &
218                     table(:,2,emplacement)*(table(:,1,emplacement+1)-absc)/&
219                     (table(:,1,emplacement+1)-table(:,1,emplacement)) + &
220                     table(:,2,emplacement+1)*(absc - table(:,1,emplacement))/&
221                     (table(:,1,emplacement+1)-table(:,1,emplacement))
222             
223                 
224          ELSE IF ((table(1,2,emplacement+1).EQ. -999.0).AND. (calcul(1) .EQ. .FALSE.)) THEN
225            !emplacement ne change pas
226            calcul(:) = .TRUE.
227            resultat(:) = table(:,2,emplacement-1)*(table(:,1,emplacement)-absc)/ &
228                     (table(:,1,emplacement)-table(:,1,emplacement-1)) + &
229                     table(:,2,emplacement)*(absc - table(:,1,emplacement-1))/&
230                     (table(:,1,emplacement)-table(:,1,emplacement-1))
231
232            resultat(:) = table(:,2,emplacement)
233                 
234          END IF
235             
236        ELSE IF (emplacement .EQ. nbpoints -1) then
237          emplacement = emplacement + 1
238          calcul(:) = .TRUE.
239          resultat(:) = table(:,2,nbpoints-1)*(table(:,1,nbpoints)-absc)/ &
240               (table(:,1,nbpoints)-table(:,1,nbpoints-1)) + &
241               table(:,2,nbpoints)*(absc - table(:,1,nbpoints-1))/&
242               (table(:,1,nbpoints)-table(:,1,nbpoints-1))
243          emplacement = nbpoints
244           
245
246        END IF
247      ELSE IF (emplacement .EQ. 0) THEN
248        IF (absc .LT. table(1,1,1) ) THEN
249          ! emplacement ne change pas
250          calcul(:) = .TRUE.
251          resultat(:) = table(:,2,1)*(table(:,1,2)-absc)/(table(:,1,2)-table(:,1,1)) + &
252               table(:,2,2)*(absc - table(:,1,1))/(table(:,1,2)-table(:,1,1))
253          emplacement = 0
254        ELSE
255          emplacement = 1
256          calcul(:) = .TRUE. 
257          resultat(:) = &
258               table(:,2,emplacement)*(table(:,1,emplacement+1)-absc)/&
259               (table(:,1,emplacement+1)-table(:,1,emplacement)) + &
260               table(:,2,emplacement+1)*(absc - table(:,1,emplacement))/&
261               (table(:,1,emplacement+1)-table(:,1,emplacement))
262        END IF
263
264      ELSE IF (emplacement .EQ. nbpoints)THEN
265
266        calcul(:) = .TRUE.
267        resultat(:) = table(:,2,nbpoints-1)*(table(:,1,nbpoints)-absc)/ &
268           (table(:,1,nbpoints)-table(:,1,nbpoints-1)) + &
269           table(:,2,nbpoints)*(absc - table(:,1,nbpoints-1))/&
270           (table(:,1,nbpoints)-table(:,1,nbpoints-1))
271        emplacement = nbpoints
272      ELSE
273        STOP 'erreur avec l''interpolation'
274
275      END IF
276
277    ENDIF
278
279  END SUBROUTINE interpolation_extrapolation
280
281  !!!!!!!!!!!!!!!!!!!!!!!!!!!!
282  !!!!!! PM
283  !!!!!!!!!!!!!!!!!!!!!!!!!!!!
284
285  SUBROUTINE Pm(npts, Pmform,alphaform,Tlform, Co2curr, pm20, pmn, pmc)
286
287    INTEGER(i_std), INTENT(in) :: npts
288    REAL(R_std), DIMENSION(npts), INTENT(in) :: pm20
289    REAL(R_std), DIMENSION(npts), INTENT(in) :: pmc
290    REAL(R_std), DIMENSION(npts), INTENT(in) :: pmn
291    REAL(R_std), DIMENSION(npts), INTENT(out) :: Pmform
292    REAL(R_std), DIMENSION(npts), INTENT(out) :: alphaform
293    REAL(R_std), DIMENSION(npts), INTENT(in)  :: Tlform
294    REAL(R_std), DIMENSION(npts), INTENT(in)  :: Co2curr
295   
296
297    REAL(R_std)  :: PmT0p
298    REAL(R_std), DIMENSION(npts) :: PmT
299    REAL(R_std), DIMENSION(npts) :: Vcmax
300    REAL(R_std), DIMENSION(npts) :: Vomax
301    REAL(R_std), DIMENSION(npts) :: Kc
302    REAL(R_std), DIMENSION(npts) :: Ko
303    REAL(R_std), DIMENSION(npts) :: Oi
304    REAL(R_std), DIMENSION(npts) :: Ci
305    REAL(R_std), DIMENSION(npts) :: Ci350
306    REAL(R_std), DIMENSION(npts) :: GAMMAstar
307    REAL(R_std), DIMENSION(npts) :: Amax
308    REAL(R_std), DIMENSION(npts) :: Amax350
309    REAL(R_std), DIMENSION(npts) :: PmCO2
310    REAL(R_std), DIMENSION(npts) :: Eps
311    REAL(R_std), DIMENSION(npts) :: Eps350
312    REAL(R_std), DIMENSION(npts) :: alphaCO2
313    REAL(R_std), DIMENSION(npts) :: Tlcel
314   
315    REAL(R_std), PARAMETER :: pmqft        = 1.5
316    REAL(R_std), PARAMETER :: pmtopt       = 303.0
317    REAL(R_std), PARAMETER :: pmt0         = 273.0
318
319
320
321
322    !Abhaengigkeit von Klimaaenderung (CO2 & T)
323
324
325    Vcmax = 98.0*EXP(68000.0*(Tlform-298.15)/ (298.15*Tlform*Rgas))*SQRT(Tlform/298.15)
326    Vomax = 0.21*Vcmax
327    Kc    = 460.0*EXP(65800*(Tlform-298.15)/(298.15*Tlform*Rgas))*SQRT(Tlform/298.15)
328    Ko    = 330.0*EXP(1400*(Tlform-298.15)/(298.15*Tlform*Rgas))*SQRT(Tlform/298.15)
329    Tlcel = Tlform - 273.15
330    Oi    = 210*(0.047 - 0.0013087*Tlcel + 0.000025603*Tlcel**2 - 0.00000021441*Tlcel**3)/0.026934
331    Ci    = 0.7*CO2curr*((1.6740 - 0.061294*Tlcel + 0.0011688*Tlcel**2 - 0.0000088741*Tlcel**3)/0.73547)
332    Ci350 = 0.7*350.0*((1.6740 - 0.061294*Tlcel + 0.0011688*Tlcel**2 - 0.0000088741*Tlcel**3)/0.73547)
333    GAMMAstar = 0.5*Vomax*Kc*Oi/(Vcmax*Ko)
334    Amax    = (Ci    - GAMMAstar)*Vcmax/(Ci    + Kc*(1.0 + Oi/Ko))
335    Amax350 = (Ci350 - GAMMAstar)*Vcmax/(Ci350 + Kc*(1.0 + Oi/Ko))
336    PmCO2 = Vcmaxadap*Amax/Amax350
337    Eps    = (ABSORvl/2.1)*(Ci    - GAMMAstar)/ (4.5*Ci    + 10.5*GAMMAstar)
338    Eps350 = (ABSORvl/2.1)*(Ci350 - GAMMAstar)/ (4.5*Ci350 + 10.5*GAMMAstar)
339    alphaCO2 = Eps/Eps350
340
341    alphaform = alpha350*alphaCO2
342
343    !Abhaengigkeit von der Temperatur
344
345    PmT0p = ((1.0 + Pmqft)*PmTopt - PmT0)/Pmqft
346    CALL fTsigm(npts, Tlform, PmT0, PmT0p, Pmqft, PmT)
347
348    Pmform = Pm20*PmN*PmC*PmCO2*PmT
349
350  END SUBROUTINE Pm
351
352
353  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
354  !  Fonctions calculant es et s : la pression de vapeur à saturation
355  ! et sa dérivée
356  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
357
358
359  SUBROUTINE es_fonct(npts,temp, resultat)
360
361    INTEGER(i_std), INTENT(in) :: npts
362    REAL(r_std), DIMENSION(npts), INTENT(in)  :: temp
363    REAL(r_std), DIMENSION(npts), INTENT(out) :: resultat
364
365
366    resultat = 0.6112*EXP(17.67*(temp/(temp+246.2)))
367
368  END SUBROUTINE es_fonct
369
370
371  SUBROUTINE s_fonct(npts,temp, resultat)
372   
373    INTEGER(i_std), INTENT(in) :: npts
374    REAL(r_std), DIMENSION(npts), INTENT(in)  :: temp
375    REAL(r_std), DIMENSION(npts), INTENT(out) :: resultat
376   
377    REAL(r_std), DIMENSION(npts) :: es_cal
378
379    CALL es_fonct(npts, temp, es_cal)
380    resultat = es_cal*17.67*246.2/(temp+246.2)**2
381  END SUBROUTINE s_fonct
382
383  ! Subroutine linreg_pasim for autogestion
384
385  SUBROUTINE linreg_pasim( npts, MaxElems, xx, yy, nDataPairs, misVal, &
386     muX, muY, sigX, sigY, a, b, cor )
387
388    INTEGER(i_std), INTENT(in)                                :: npts
389    INTEGER(i_std), INTENT(in)                                :: MaxElems
390    REAL(r_std), DIMENSION(npts,MaxElems), INTENT(in) :: xx     
391    ! IN :x-data
392    REAL(r_std), DIMENSION(npts,MaxElems), INTENT(in) :: yy       
393    ! IN :y-data
394    INTEGER(i_std), INTENT(in)                        :: nDataPairs 
395    ! IN : number o valid x-y pairs
396    REAL(r_std), INTENT(in)                           :: misVal 
397    ! IN :missing value
398    REAL(r_std), DIMENSION(npts), INTENT(out)         :: muX, muY 
399    ! OUT: meanof x and y
400    REAL(r_std), DIMENSION(npts), INTENT(out)         :: sigX, sigY
401    ! OUT:standard deviation of x and y
402    REAL(r_std), DIMENSION(npts), INTENT(out)         :: a, b   
403    ! OUT:parameters of least-square                                                               
404    !      fit y=a*x+b,
405    REAL(r_std), DIMENSION(npts), INTENT(out)         :: cor     
406    ! OUT:correlation between x and y
407
408    INTEGER(i_std)     , DIMENSION(npts)  ::  nn
409    INTEGER(i_std) ::  ii, i
410    REAL(r_std), DIMENSION(npts)  ::  sxy
411
412
413    ! initialize variables
414    nn   = 0
415    muX  = 0.0
416    muY  = 0.0
417    sigX = 0.0
418    sigY = 0.0
419    sxy  = 0.0
420
421    DO ii=1,nDataPairs
422
423      WHERE ((xx(:,ii).NE.misVal).AND.(yy(:,ii).NE.misVal))
424        nn(:)   = nn(:) + 1
425        muX(:)  = muX(:) + xx(:,ii)
426        muY(:)  = muY(:) + yy(:,ii)
427        sigX(:) = sigX(:) + xx(:,ii)*xx(:,ii)
428        sigY(:) = sigY(:) + yy(:,ii)*yy(:,ii)
429      END WHERE
430
431    END DO
432    DO i=1,npts
433      IF (nn(i) .GT. 1) THEN
434        muX(i)  = muX(i)/nn(i)
435        muY(i)  = muY(i)/nn(i)
436        sigX(i) = sigX(i)-nn(i)*muX(i)*muX(i)
437        sigY(i) = sigY(i)-nn(i)*muY(i)*muY(i)
438
439        IF (sigX(i) .LT. 0.0)  sigX(i) = 0.0
440        IF (sigY(i) .LT. 0.0)  sigY(i) = 0.0
441
442        DO ii=1,nDataPairs
443          IF ((xx(i,ii).NE.misVal) .AND. (yy(i,ii) .NE. misVal)) THEN
444              sxy(i) = sxy(i)+(xx(i,ii)-muX(i))*(yy(i,ii)-muY(i))
445          END IF
446        END DO
447
448        IF ((sigX(i).GT.0.0).AND.(sigY(i).GT.0.0)) THEN
449            a(i)   = sxy(i)/sigX(i)
450            b(i)   = muY(i) - a(i)*muX(i)
451            cor(i) = sxy(i)/SQRT(sigX(i)*sigY(i))
452        ELSE
453            a(i)   = misVal
454            b(i)   = misVal
455            cor(i) = misVal
456        END IF
457
458        sigX(i) = SQRT(sigX(i)/REAL(nn(i)-1))
459        sigY(i) = SQRT(sigY(i)/REAL(nn(i)-1))
460
461      ELSE
462
463        muX(i)  = misVal
464        muY(i)  = misVal
465        sigX(i) = misVal
466        sigY(i) = misVal
467        a(i)    = misVal
468        b(i)    = misVal
469        cor(i)  = misVal
470
471      END IF
472    END DO
473
474
475  END SUBROUTINE linreg_pasim  ! LinReg
476
477END MODULE grassland_fonctions
478
479
480
481
482
483
484
Note: See TracBrowser for help on using the repository browser.