source: LMDZ6/trunk/libf/phylmd/ocean_albedo_mod.f90 @ 6087

Last change on this file since 6087 was 6048, checked in by fhourdin, 3 months ago

Renommage des nom de fichiers incluant des modules
File size: 10.3 KB
Line 
1!
2! $Id$
3!
4MODULE ocean_albedo_mod
5
6CONTAINS
7
8SUBROUTINE ocean_albedo(knon,zrmu0,knindex,pwind,SFRWL,alb_dir_new,alb_dif_new)
9!$gpum horizontal knon klon
10
11!!
12!!****  *ALBEDO_RS14* 
13!!
14!!    PURPOSE
15!!    -------
16!!     computes the direct & diffuse albedo over open water
17!!     
18!!**  METHOD
19!!    ------
20!!
21!!    EXTERNAL
22!!    --------
23!!
24!!    IMPLICIT ARGUMENTS
25!!    ------------------
26!!
27!!    REFERENCE
28!!    ---------
29!!
30!!    AUTHOR
31!!    ------
32!!      R. Séférian           * Meteo-France *
33!!
34!!    MODIFICATIONS
35!!    -------------
36!!      Original    03/2014
37!!                  05/2014 R. Séférian & B. Decharme :: Adaptation to spectral
38!!                  computation for diffuse and direct albedo
39!!                  08/2014 S. Baek :: for wider wavelength range 200-4000nm and
40!!                  adaptation to LMDZ + whitecap effect by Koepke + chrolophyll
41!!                  map from climatology file
42!!                  10/2016 O. Boucher :: some optimisation following R.
43!!                  Seferian's work in the CNRM Model
44!!
45!-------------------------------------------------------------------------------
46!
47!*           DECLARATIONS
48!            ------------
49!
50USE ocean_albedo_para, ONLY : nnwl, xakwl, xakrefm, xakachl, xakaw3, xakbw, xaw440, xfrwl, xrwc   
51USE phys_state_var_mod, ONLY : chl_con
52USE clesphys_mod_h, ONLY: nsw, ok_chlorophyll
53!
54!
55IMPLICIT NONE
56!
57!*      0.1    declarations of arguments
58!              -------------------------
59!
60!
61INTEGER, INTENT(IN) :: knon
62INTEGER, DIMENSION(knon), INTENT(IN) :: knindex
63REAL, DIMENSION(knon), INTENT(IN) :: zrmu0         !--cos(SZA) on full vector
64REAL, DIMENSION(knon), INTENT(IN) :: pwind         !--wind speed on compressed vector
65REAL, DIMENSION(6),INTENT(IN) :: SFRWL
66REAL, DIMENSION(knon,nsw), INTENT(OUT) :: alb_dir_new, alb_dif_new
67!
68!*      0.2    declarations of local variables
69!              -------------------------
70!
71REAL, DIMENSION(knon)           :: ZCHL        ! surface chlorophyll
72REAL, DIMENSION(knon)           :: ZCOSZEN     ! Cosine of the zenith solar angle
73!
74INTEGER                         :: JWL, INU    ! indexes
75INTEGER                         :: JI
76REAL                            :: ZWL         ! input parameter: wavelength and diffuse/direct fraction of light
77REAL:: ZCHLABS, ZAW, ZBW, ZREFM, ZYLMD, ZUE, ZUE2 ! scalar computation variables
78!
79REAL, DIMENSION(knon) :: ZAP, ZXX2, ZR00, ZRR0, ZRRR               ! computation variables
80REAL, DIMENSION(knon) :: ZR22, ZR11DF                              ! computation variables
81REAL, DIMENSION(knon) :: ZBBP, ZNU, ZHB                            ! computation variables
82REAL, DIMENSION(knon) :: ZR11, ZRW, ZRWDF, ZRDF                    ! 4 components of the OSA
83REAL, DIMENSION(knon) :: ZSIG, ZFWC, ZWORK1, ZWORK2, ZWORK3
84!
85!--initialisations-------------
86!
87
88IF (knon==0) RETURN ! A verifier pourquoi on en a besoin...
89
90alb_dir_new(:,:) = 0.
91alb_dif_new(:,:) = 0.
92!
93! Initialisation of chlorophyll content
94! ZCHL(:) = CHL_CON!0.05 ! averaged global values for surface chlorophyll
95IF (ok_chlorophyll) THEN
96  ZCHL(1:knon)=CHL_CON(knindex(1:knon))
97ELSE
98  ZCHL(1:knon) = 0.05
99ENDIF
100
101! variables that do not depend on wavelengths
102! loop over the grid points
103! functions of chlorophyll content
104ZWORK1(1:knon)= EXP(LOG(ZCHL(1:knon))*0.65)
105ZWORK2(1:knon)= 0.416 * EXP(LOG(ZCHL(1:knon))*0.766)
106ZWORK3(1:knon)= LOG10(ZCHL(1:knon))
107! store the cosine of the solar zenith angle
108ZCOSZEN(1:knon) = zrmu0(1:knon)
109! Compute sigma derived from wind speed (Cox & Munk reflectance model)
110ZSIG(1:knon)=SQRT(0.003+0.00512*PWIND(1:knon))
111! original : correction for foam (Eq 16-17)
112! has to be update once we have information from wave model (discussion with G. Madec)
113ZFWC(1:knon)=3.97e-4*PWIND(1:knon)**1.59 ! Salisbury 2014 eq(2) at 37GHz, value in fraction
114!
115DO JWL=1,NNWL           ! loop over the wavelengths
116  !
117  !---------------------------------------------------------------------------------
118  ! 0- Compute baseline values
119  !---------------------------------------------------------------------------------
120   
121  ! Get refractive index for the correspoding wavelength
122  ZWL=XAKWL(JWL)      !!!--------- wavelength value
123  ZREFM= XAKREFM(JWL) !!!--------- refraction index value
124 
125  !---------------------------------------------------------------------------------
126  ! 1- Compute direct surface albedo (ZR11)
127  !---------------------------------------------------------------------------------
128  !
129  ZXX2(1:knon)=SQRT(1.0-(1.0-ZCOSZEN(1:knon)**2)/ZREFM**2)
130  ZRR0(1:knon)=0.50*(((ZXX2(1:knon)-ZREFM*ZCOSZEN(1:knon))/(ZXX2(1:knon)+ZREFM*ZCOSZEN(1:knon)))**2 +  &
131               ((ZCOSZEN(1:knon)-ZREFM*ZXX2(1:knon))/(ZCOSZEN(1:knon)+ZREFM*ZXX2(1:knon)))**2)
132  ZRRR(1:knon)=0.50*(((ZXX2(1:knon)-1.34*ZCOSZEN(1:knon))/(ZXX2(1:knon)+1.34*ZCOSZEN(1:knon)))**2 + &
133               ((ZCOSZEN(1:knon)-1.34*ZXX2(1:knon))/(ZCOSZEN(1:knon)+1.34*ZXX2(1:knon)))**2)
134  ZR11(1:knon)=ZRR0(1:knon)-(0.0152-1.7873*ZCOSZEN(1:knon)+6.8972*ZCOSZEN(1:knon)**2-8.5778*ZCOSZEN(1:knon)**3+ &
135               4.071*ZSIG(1:knon)-7.6446*ZCOSZEN(1:knon)*ZSIG(1:knon)) *  &
136               EXP(0.1643-7.8409*ZCOSZEN(1:knon)-3.5639*ZCOSZEN(1:knon)**2-2.3588*ZSIG(1:knon)+ &
137               10.0538*ZCOSZEN(1:knon)*ZSIG(1:knon))*ZRR0(1:knon)/ZRRR(1:knon)
138  !
139  !---------------------------------------------------------------------------------
140  ! 2- Compute surface diffuse albedo (ZRDF)
141  !---------------------------------------------------------------------------------
142  ! Diffuse albedo from Jin et al., 2006 + estimation from diffuse fraction of
143  ! light (relying later on AOD). CNRM model has opted for Eq 5b
144  ZRDF(1:knon)=-0.1482-0.012*ZSIG(1:knon)+0.1609*ZREFM-0.0244*ZSIG(1:knon)*ZREFM ! surface diffuse (Eq 5a)
145  !!ZRDF(1:knon)=-0.1479+0.1502*ZREFM-0.0176*ZSIG(1:knon)*ZREFM   ! surface diffuse (Eq 5b)
146 
147  !---------------------------------------------------------------------------------
148  ! *- Determine absorption and backscattering
149  ! coefficients to determine reflectance below the surface (Ro) once for all
150  !
151  ! *.1- Absorption by chlorophyll
152  ZCHLABS= XAKACHL(JWL)
153  ! *.2- Absorption by seawater
154  ZAW= XAKAW3(JWL)
155  ! *.3- Backscattering by seawater
156  ZBW= XAKBW(JWL)
157  ! *.4- Backscattering by chlorophyll
158  ZYLMD = EXP(0.014*(440.0-ZWL))
159  ZAP(1:knon) = 0.06*ZCHLABS*ZWORK1(1:knon) +0.2*(XAW440+0.06*ZWORK1(1:knon))*ZYLMD
160   
161!!  WHERE ( ZCHL(1:knon) > 0.02 )
162!!    ZNU(:)=MIN(0.0,0.5*(ZWORK3(:)-0.3))
163!!    ZBBP(:)=(0.002+0.01*(0.5-0.25*ZWORK3(:))*(ZWL/550.)**ZNU(:))*ZWORK2(:)
164!!  ELSEWHERE
165!!    ZBBP(:)=0.019*(550./ZWL)*ZWORK2(:)       !ZBBPf=0.0113 at chl<=0.02
166!!  ENDWHERE
167
168    do JI = 1, knon
169      IF (ZCHL(JI) > 0.02) THEN
170        ZNU(JI)=MIN(0.0,0.5*(ZWORK3(JI)-0.3))
171        ZBBP(JI)=(0.002+0.01*(0.5-0.25*ZWORK3(JI))*(ZWL/550.)**ZNU(JI)) &
172                  *ZWORK2(JI)
173      ELSE
174        ZBBP(JI)=0.019*(550./ZWL)*ZWORK2(JI)       !ZBBPf=0.0113 at chl<=0.02
175      ENDIF
176    ENDDO
177
178  ! Morel-Gentili(1991), Eq (12)
179  ! ZHB=h/(h+2*ZBBPf*(1.-h))       
180  ZHB(1:knon)=0.5*ZBW/(0.5*ZBW+ZBBP(1:knon))
181   
182  !---------------------------------------------------------------------------------
183  ! 3- Compute direct water-leaving albedo (ZRW)
184  !---------------------------------------------------------------------------------
185  ! Based on Morel & Gentilli 1991 parametrization
186  ZR22(1:knon)=0.48168549-0.014894708*ZSIG(1:knon)-0.20703885*ZSIG(1:knon)**2
187
188  ! Use Morel 91 formula to compute the direct reflectance
189  ! below the surface
190  ZR00(1:knon)=(0.5*ZBW+ZBBP(1:knon))/(ZAW+ZAP(1:knon)) *  &
191               (0.6279-0.2227*ZHB(1:knon)-0.0513*ZHB(1:knon)**2 + &
192               (-0.3119+0.2465*ZHB(1:knon))*ZCOSZEN(1:knon))
193  ZRW(1:knon)=ZR00(1:knon)*(1.-ZR22(1:knon))/(1.-ZR00(1:knon)*ZR22(1:knon))
194
195  !---------------------------------------------------------------------------------
196  ! 4- Compute diffuse water-leaving albedo (ZRWDF)
197  !---------------------------------------------------------------------------------
198  ! as previous water-leaving computation but assumes a uniform incidence of
199  ! shortwave at surface (ue)
200  ZUE=0.676               ! equivalent u_unif for diffuse incidence
201  ZUE2=SQRT(1.0-(1.0-ZUE**2)/ZREFM**2)
202  ZRR0(1:knon)=0.50*(((ZUE2-ZREFM*ZUE)/(ZUE2+ZREFM*ZUE))**2 +((ZUE-ZREFM*ZUE2)/(ZUE+ZREFM*ZUE2))**2)
203  ZRRR(1:knon)=0.50*(((ZUE2-1.34*ZUE)/(ZUE2+1.34*ZUE))**2 +((ZUE-1.34*ZUE2)/(ZUE+1.34*ZUE2))**2)
204  ZR11DF(1:knon)=ZRR0(1:knon)-(0.0152-1.7873*ZUE+6.8972*ZUE**2-8.5778*ZUE**3+4.071*ZSIG(1:knon)-7.6446*ZUE*ZSIG(1:knon)) * &
205                 EXP(0.1643-7.8409*ZUE-3.5639*ZUE**2-2.3588*ZSIG(1:knon)+10.0538*ZUE*ZSIG(1:knon))*ZRR0(1:knon)/ZRRR(1:knon)
206
207  ! Use Morel 91 formula to compute the diffuse
208  ! reflectance below the surface
209  ZR00(1:knon) = (0.5*ZBW+ZBBP(1:knon)) / (ZAW+ZAP(1:knon)) &
210       * (0.6279-0.2227*ZHB(1:knon)-0.0513*ZHB(1:knon)**2 &
211       + (-0.3119+0.2465*ZHB(1:knon))*ZUE)
212  ZRWDF(1:knon)=ZR00(1:knon)*(1.-ZR22(1:knon))*(1.-ZR11DF(1:knon))/(1.-ZR00(1:knon)*ZR22(1:knon))
213   
214  ! get waveband index inu for each nsw band
215  SELECT CASE(nsw)
216  CASE(2)
217    IF (JWL.LE.49) THEN       ! from 200  to 680 nm
218     inu=1
219    ELSE                      ! from 690  to 4000 nm
220     inu=2
221    ENDIF
222  CASE(4)
223    IF (JWL.LE.49) THEN       ! from 200  to 680 nm
224     inu=1
225    ELSE IF (JWL.LE.99) THEN  ! from 690  to 1180 nm
226     inu=2
227    ELSE IF (JWL.LE.218) THEN ! from 1190 to 2370 nm
228     inu=3
229    ELSE                      ! from 2380 to 4000 nm
230     inu=4
231    ENDIF
232  CASE(6)
233    IF (JWL.LE.5) THEN        ! from 200  to 240 nm
234     inu=1
235    ELSE IF (JWL.LE.24) THEN  ! from 250  to 430 nm
236     inu=2
237    ELSE IF (JWL.LE.49) THEN  ! from 440  to 680 nm
238     inu=3
239    ELSE IF (JWL.LE.99) THEN  ! from 690  to 1180 nm
240     inu=4
241    ELSE IF (JWL.LE.218) THEN ! from 1190 to 2370 nm
242     inu=5
243    ELSE                      ! from 2380 to 4000 nm
244     inu=6
245    ENDIF
246  END SELECT
247
248  ! partitionning direct and diffuse albedo
249  ! excluding diffuse albedo ZRW on ZDIR_ALB
250
251  !--direct
252  alb_dir_new(1:knon,inu)=alb_dir_new(1:knon,inu) + &
253                          ( XFRWL(JWL) * ((1.-ZFWC(1:knon)) * (ZR11(1:knon)+ZRW(1:knon))   + ZFWC(1:knon)*XRWC(JWL)) )/SFRWL(inu)
254  !--diffuse
255  alb_dif_new(1:knon,inu)=alb_dif_new(1:knon,inu) + &
256                          ( XFRWL(JWL) * ((1.-ZFWC(1:knon)) * (ZRDF(1:knon)+ZRWDF(1:knon)) + ZFWC(1:knon)*XRWC(JWL)) )/SFRWL(inu)
257
258ENDDO ! ending loop over wavelengths
259
260END SUBROUTINE ocean_albedo
261
262END MODULE ocean_albedo_mod
Note: See TracBrowser for help on using the repository browser.