source: LMDZ5/trunk/libf/phylmd/ocean_albedo.F90 @ 2680

Last change on this file since 2680 was 2680, checked in by fhourdin, 8 years ago

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