!       ifort -132 -assume byterecl get_Rw0p_lambda_Cnorm_theta.F90
	program get_Rwp
	
	implicit none
	real bup0, b_b0, b_bp, b_b, bb0bb, gamma, b_0, K0l, bb
	real Kl, Xl, C, el, al, u, Tw, Rw0pl, pi, pi2, u1, R, Rp, Fu, dFudu, R0
	real sunz, theta, rad 
	integer nlmorel, nlsmith, nlt, lnc, l, n, i, ii, k, it
	character ct*2

	parameter(nlmorel=91, nlsmith=nlmorel, nlt=nlmorel, lnc=34)
	real Ci(15)
	real Cr, du, u3, uf0, Rf0	

	real Rwp(nlmorel), Rmo(nlmorel)				
	real Lambda, Kwl, Kdl, mu0
	real Lambda0(nlmorel), Kwl0(nlmorel), Xl0(nlmorel), el0(nlmorel) 
	real Lambdas(nlsmith), Kwls(nlsmith), aws(nlsmith), bmsw(nlsmith), bmfw(nlsmith)

        data Ci/0.00,0.02,0.03,0.05,0.10,0.20,0.30,0.50,1.0,2.0,3.0,5.0,10.0,20.0,30.0/

	pi = dacos(-1.0D0)
	pi2 = pi*2.0
	rad = 180.0D0/pi

	u1     = 0.75
	Tw     = 0.526
	bup0   = 0.30	! (1/m)

	open(10, file='Table_Morel_5nm_300-750.data', form='formatted', status='old')	!91 5nm
	open(12, file='Table_Smith_5nm_300-750.data', form='formatted', status='old')	!91 5nm 
!!
	do l=1,nlmorel
	read(10,*,ERR=999)Lambda0(l), Kwl0(l), Xl0(l), el0(l)			! Morel table
	write(8,'(i3,4F10.4)')l, Lambda0(l), Kwl0(l), Xl0(l), el0(l)
!	el = el0(l)	!Table2_morel.txt
!	cp2=(C**el)
	enddo
!!
      do l=1,nlsmith
	read(12,*,ERR=999)Lambdas(l), Kwls(l),aws(l), bmsw(l)			! Smith table 
	write(9,'(i3,5F10.4)')l, Lambdas(l), Kwls(l),aws(l), bmsw(l)
	enddo
	
	theta = 0.		! zenith angle

	do k=1,3
	theta = 30.*(k-1)
	it = theta 
	write(ct,'(i2.2)')it
	sunz  = theta
	mu0 = cos(sunz/rad)

	open(34,file='Rwp_chl_'//ct//'_nolog_e_5nm.datt', form='unformatted', &
	     access='direct', status='unknown', recl=4*91)

	do i = 1,15
	C   = Ci(i)

	do l = 1,nlmorel
	b_0 = bmsw(l)
	K0l = Kwl0(l)
	Xl  = Xl0(l)
	el  = el0(l)
	Lambda = Lambda0(l)

	b_b0 = 0.5*bmsw(l)		!! b_0: table Smith 1981; Morel (1974)
	if (C .eq. 0.)then
	b_bp = 0.
	else
!*wrong	b_bp = bup0*(C**0.62)*(0.002+0.02*(0.5-0.25*log10(C))/(550./lambda))
	b_bp = bup0*(C**0.62)*(0.002+0.02*(0.5-0.25*log10(C))*(550./lambda))
	endif 

	b_b  = b_b0+b_bp

	bb0bb = b_b0/b_b 
	gamma = 0.6279-0.2227*(bb0bb)-0.0513*(bb0bb*bb0bb)-0.3119*mu0 + 0.2465*(bb0bb)*mu0

	if (C .eq. 0.)then
	Kl = K0l
	else
	Kl = K0l+Xl*(C**el)	!! K0l, Xl, el: table Morel (1974)
	endif
!* 
!* Iterative solution: solve for Morel's u
!*
	bb  = b_b
	Kdl = Kl
	Cr = 0.33*(bb/Kdl)
	u = u1

	do n = 1,10			! Iterative sol
	R = Cr/u			! (0.33/u)*(bb/Kdl)
	Rp = -Cr/(u*u)
	Fu = u - 0.90*(1.-R)/(1.+2.25*R)
	dFudu = 1.0 + 0.90*3.25*Rp/((1.+2.25*R)*(1.+2.25*R))
	u3 = u - (Fu/dFudu)
	du = u3-u
	write(6,'(a,3f10.5,i3)')'Itera Sol: u3 du R n = ', u3, du, R, n
	if(abs(du) .le. 1e-8) goto 10
	u = u3
	enddo

   10	continue

	R0 = Cr/u3
	uf0 = 1./u3	! Frouin def
	al = Kdl/uf0
!xx
!!	u3 = 0.75
!!	al = Kdl*u3 	! Morel def
!xx
	Rwp(l) = Tw*gamma*b_b/al
	write(7,'(a,8f12.5)')'Tw ...  =', Tw,gamma,Kdl,b_b,al,R0,uf0,u3

	write(6,'(a,i3,f10.5,i5,f10.0,f15.5)')'i C l lm Rw0pl =', i,C,l,lambda,Rw0pl
	enddo		!! l = lambda

	write(34,rec=i)(Rwp(ii),ii=1,91)	!! No log10 values; 300-750nm; 91 lambdas	

	write(6,'(x)')
!!
	enddo 		!! i = Chl
!!
	close (34)
	enddo		!! k=1,3 theta
  999	continue
	end
