N        H    T   RH    U    V   T0   T0 is temperature in previous time step
       50.0  6.8 97.1  1.8  2.3  7.1  H is level's geographic height (above the see level)
      100.0  7.0 97.3  3.0  3.8  7.5
      200.0  6.5 83.5  3.5  3.4  6.8
      300.0  6.0 60.2  4.2  4.2  6.5
      450.0  5.4 51.3  5.2  4.4  5.4
      600.0  4.7 44.0  6.1  5.5  4.7
      800.0  3.3 33.0  7.2  5.0  3.3
     1000.0 -0.9 33.5  8.3  6.6 -0.9
     1300.0 -1.2 32.0  9.4  6.9 -1.2
     1600.0 -2.5 32.1  9.4  5.4 -2.5
     2000.0 -3.4 31.5  9.6  4.6 -3.4
     2300.0 -4.9 30.7  9.6  6.7 -4.9
     2600.0 -5.2 24.1  9.8  2.9 -5.2
     2900.0 -6.4 20.6  9.9  1.5 -6.4
       HSFC  T2m RH2M  U10  V10 T2m0
       30.0  5.5 97.5  1.0  1.2  6.7
    20000.0   20000.0                   clout_top and base heights (above the ground), very large value 20000 menas no cloud
   -0.5000/10000                         adv (in [g/kg]/s) = -[u(dq/dx)+v(dq/dy)], adv>0:wet advection, < 0 dry advection






=============================================================================================================
                                        CODE DESCRIPTION
==============================================================================================================
**Bcckground**
This code is based on the ipaper "Asmptotic Analysis of Equilibrium  of Radiation Fog" by 
Zhou and Ferrier (JAMC, 2008). To extend this theory from radiation fog to advection fog, 
moisture advection term and low cloud to deal with advection fog and low stratus are added.
The LWC computation is based on saturated levels (RH can be < 100%) and cooling rate in whatever
case where LWC exists or not. This means, even if there is no LWC but if the layers
are saturated and there is cooling there, the fog LWC can be computed out. If there
already exist, the existing LWC is only used to compute advection.

The major subroutine is "get_new_fog" which can be used in a mesoscale fog directly but you
should be careful to make the levels match well and get advection computed before calling
this subroutine.

Note: Above data are from model grid output, you can use model output or observed profile data.
  H,T,RH,U,V,T0 are vertical profiles for level height(above sea level,m), temperature (C), RH (%), 
x and y wind components (m/s), and previous hours time step's temperature (C). HSFC, T2m,RH2M, U10, V10,T2m0  
are Surface height, 2m temperature (C), 2m RH(%), 10m wind speed in x and y components, and 2m previous hours 
time step's temperature (C).

20000.0 and 20000.0 are cloud top and base heights (above the surface level, m).
20000.0 means there is no cloud in this case. 


-0.5000/10000 is LWC/moisture advection term in g/kg/sec (>0 wet advection, <0 dry advection),
where the "10000" will not be read into the code. If no advection, you set it be 0.0000/10000

There are 14 levels set for this code. If your data levels is not this case, you can modify
'lv= your levels'  in the code and data as well as levels  in this file. 

Another note: time step is in hour (3600 will be divided in side the code).    
This code assumes that the time step is 1 hour (interval=1). If your data is no this case,
modify the variable "interval= your case". 
 
The testing data here came from a profile at a grid point in an real model (NAM) execution.
But only 14 levels are cut-off to input here. The number of levels can be changed freely.
You also can change profile data (cloud top, base, adv, saturated layers, etc) to test the code.
95% is assumed as saturated in side the code, you can re-set based on your own situation.
If you won't consider cloud stratus, just set clout top to be 1000.0 in this file
If you won't consider advection, just set adv to be 0.0 in this file
Please note the data are formatted here, be sure make them in correct location in this file
(check the read format inside the code).

How to run the code:
compiler f77, f90, f95 on Linux
example: f77 test_new_fog.f
Execution: a.out<test_new_fog.data     output will be on screen 
or a.out<test_new_fog.data>output      output will be in the output file 
 
If you don't change anything for this data, running the code will get following output:

------------------------------------------------------------------------------ 
[wd20bz@lnx255]> test_new_fog.x < fog_test.office_note
 rhthr=  95.
 Input data:
     2900.0 -6.4 20.6  9.9  1.5 -6.4
     2600.0 -5.2 24.1  9.8  2.9 -5.2
     2300.0 -4.9 30.7  9.6  6.7 -4.9
     2000.0 -3.4 31.5  9.6  4.6 -3.4
     1600.0 -2.5 32.1  9.4  5.4 -2.5
     1300.0 -1.2 32.0  9.4  6.9 -1.2
     1000.0 -0.9 33.5  8.3  6.6 -0.9
      800.0  3.3 33.0  7.2  5.0  3.3
      600.0  4.7 44.0  6.1  5.5  4.7
      450.0  5.4 51.3  5.2  4.4  5.4
      300.0  6.0 60.2  4.2  4.2  6.5
      200.0  6.5 83.5  3.5  3.4  6.8
      100.0  7.0 97.3  3.0  3.8  7.5
       50.0  6.8 97.1  1.8  2.3  7.1
       30.0  5.5 97.5  1.0  1.2  6.7
   20000.0   20000.0
-.5000E-04
 Output:
 z10m=  40. z2m=  32.
 kp= 0
 z_RH-hsfc=  70.
 touch_down=  0.
 z_RH=  100. kfog= 2 rh2=  97.5 z2m=  32.
 depth_fog=  70. cooling_rate=  0.000236111082 tfog=  6.25
 beta=  0.45508799
  cooling_term=  0.00010745132 adv_term= -4.99999987E-05
 es=  9.5054636 Tk=  279.420013 beta=  0.45508799 lwc_max=  0.254684895
 wp=  4.84148741 w10=  1.56204998 dwdz=  0.0546572916 dTdz=  0.0220588241
 Ri=  0.376374274 FM=  0.0176481903
 Km=  0.26176554 Kc=  1.5253588 lamda=  16.4733562 FBL=  8.28872299
 LWC profile:
  0.m: LWC=  0.
  7.m: LWC=  0.088507168
  14.m: LWC=  0.148385376
  21.m: LWC=  0.175626352
  28.m: LWC=  0.18047525
  35.m: LWC=  0.172729746
  42.m: LWC=  0.157887682
  49.m: LWC=  0.138121158
  56.m: LWC=  0.113306493
  63.m: LWC=  0.0802838206
  70.m: LWC=  0.
 Surface 10m LWC(g/kg)=  0.088507168
-----------------------------------------------------------------------------------
How to compute the advection term:
Advection is layer averaged moist advection within fog. If no sounding data, use 2m
specific humidity and 10m winds. 

Mositure advection term, adv,  should be computed outside this code as following: 
This code requires total moisture (vapor + liquid) specific (g/kg) to compute this term

            do j=2,jm-1                                     !general upwind scheme
             do i=2,im-1                                    !
              qw(i,j)=q(i,j)+LWC(i,j)                       !              [i,j+1]
              qw(i+1,j)=q(i+1,j)+LWC(i+1,j)                 !                 |
              qw(i-1,j)=q(i-1,j)+LWC(i-1,j)                 !       [i-1,j]-[i,j]-[i+1,j]
              qw(i,j+1)=q(i,j)+LWC(i,j+1)                   !                 |
              qw(i,j-1)=q(i,j)+LWC(i,j-1)                   !              [i,j-1]
              if(u10(i,j).ge.0.0.and.v10(i,j).ge.0.0) then
                 qadv(i,j) =
     +            -u10(i,j)*(qw(i,j)-qw(i-1,j))/dx   !where dx,dx are grid space (in km)
     +            -v10(i,j)*(qw(i,j)-qw(i,j-1))/dy
         else if(u10(i,j).gt.0.0.and.v10(i,j).lt.0.0) then
                 qadv(i,j) =
     +            -u10(i,j)*(qw(i,j)-qw(i-1,j))/dx
     +            -v10(i,j)*(qw(i,j+1)-qw(i,j))/dy
         else if(u10(i,j).lt.0.0.and.v10(i,j).lt.0.0) then
                 qadv(i,j) =
     +            -u10(i,j)*(qw(i+1,j)-qw(i,j))/dx
     +            -v10(i,j)*(qw(i,j+1)-qw(i,j))/dy
                 end if
         else if(u10(i,j).lt.0.0.and.v10(i,j).gt.0.0) then
                  qadv(i,j) =
     +             -u10(i,j)*(qw(i+1,j)-qw(i,j))/dx
     +             -v10(i,j)*(qw(i,j)-qw(i,j-1))/dy
              end if
             end do
            end do
 

-   
Inside the subroutine, several parameters can be tuned based on your own model:
(1) Saturated threshold, current is RH=95%
(2) Turbulence paramters, current statbility function (Fm) is dependent on Ri number and
  takes form of ECMWF for Ri>0 and Hostlag for Ri<0. You can use your own formula if possible.



cc
cc     Testing get_new_fog subroutine
cc
 
        parameter (lv=14,time_step=3600.0)                     !14 levels, 1 hour time step
        REAL up(lv),vp(lv),hp(lv),rhp(lv),tp(lv,2),t2(2)
        REAL u10,v10,hsfc,lwc,rh2,rhthr,
     +    cooling_rate,Km,lamda,lwc_max

c        open(10,file='fog_test.data',status='old')

        read(*,*)
        do k=1,lv
         read(*,20) hp(k),tp(k,2),rhp(k),up(k),
     +      vp(k),tp(k,1)
20       format(5x,f6.1, 5f5.1) 
        end do
        read(*,*)
        read(*,20) hsfc,t2(2),rh2,u10,v10,t2(1)
        read(*,21) cldt, cldb
        read(*,'(f10.4)') adv
21      format(2f10.1) 

         rhr=95.0                          !assume RH>=95% as to be saturated
         dx=10000.0
         adv=adv/dx                        !in g/kg/sec

        call get_new_fog(up,vp,hp,u10,v10,hsfc,rhp,
     +       rh2,t2,tp,time_step,lv,cldt,cldb,adv,rhr,
     +       lwc)

        write(*,*) 'Surface 10m LWC(g/kg)=',lwc

        stop
        end

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
cc
cc  SUBROUTINE get_new_fog: implemented based on Zhou and Ferrier 's aysmptotic formulation
cc      2008 on JAMC but add moisture advection term                     
cc        
cc    Input: grid-wide u,v component profile and surface values, temperature and RH profiles
cc            and surface values, surface height, and previous time step temperature profile 
cc            and surface values (are used to compute cooling rate).
cc    Output: fog LWC
cc
cc    Application: added to NCEP SREF/VSREF product generator or WRF-NMM/ARW unified post
cc
cc    This code is for grid or point-wide computation as following steps:
cc 
cc         (1) Search for which (pressure)levels the surface height is located      
cc         (2) Search for which (pressure)levels the suturated top above the surface is located 
cc                and get the fog layer (saturated layer) depth, if its depth > 800m, not 
cc                a fog but reture, otherwise              
cc         (3) Search from the fog layer top downward to see if its bottom reach the ground
cc             if its bottom not touch the ground and > 800m, it is not a fog but return, 
cc             otherwise 
cc         (4) Compute averaged cooling rate within the fog layer using current and previous 
cc             time step's temperature profiles and surface (in C/hr) 
cc         (5) Based on temperaure at nearest level from the fog top and 2m to compute
cc             averag temperature with the fog layer, based on which parameter beta is computed
cc       (5.1) From cooling rate, fog depth and beta, the critical turb exch. ceoff Kc can 
cc               be computed
cc         (6) Based on the u,v at nearest level from the fog top and 10m u,v, to compute   
cc             wind speeds at the fog top and the surface
cc         (7) Based on T, wind speeds at nearest level from the fog top and 2m, Ri is computed          
cc             from which stability function Fm(Ri) is computed, from which turbulent excahnge
cc             coefficient K is computed
cc       (7.1) If won't further compute LWC, using K and Kc, fog exists or not can be determined
cc              at this step (K>Kc, K<Kc?), otherwise,     
cc         (8) With K, cooling rate, beta, fog boundary layer FBL can be computed
cc         (9) From the FBL, cooling rate, fog depth, LWC profile with the fog layer are computed 
cc        
cc    Plan to add:
cc        Add LWC advection term from pervious fog LWC field into the Kc and LWC computations
cc
cc    Language: Fortran 77/90
cc
cc    Origninal author: Binbin Zhou, EMC/NCEP/NOAA
cc                      August 5, 2010
cc
cc    Modifcation history: 
cc        Sept 9, 2010 by B. Zhou
cc            Add startus cloud layer checking (< 400m), if yes, use modeled cloud layer as fog layer
cc            instead of using RH to check saturated layer as fog layer. This imply the modeled cloud
cc            depth is assumed correct but just adjust its LWC distrbution.   
cc            If clout top < 400m and cloud bottm < 50 m above the ground,
cc            find whcih level the cloud top is located, then repeat the above from step(4).
cc
cc        Sept. 19, 2010 by B. Zhou
cc            Add advection term. Around a grid, if there are LWC in these grid, use LWC to computer
cc            advection, otherwise, use humidity specific to compute the advection. If no humidity specific,
cc            use RH to convert to humidity specific   
cc            
cc
cc   Description diagram:  
cc
cc   (1) Clear sky case but lower levels are saturated
cc
cc             kp and kp+1: the level just below and abobe the surface
cc                    kfog: fog top level 
cc
cc
cc  ---------------------------------------------------------------------------- k=14, hp(14)
cc
cc   .....................................................                       .....
cc                                        
cc  ------------ fog top---100%-------------- kfog (=4 in this case)------------ k=4, hp(4)
cc         ^                             ^   
cc         |                             |
cc   ------|-------------- 100%----------|--- kp+1 (=3 in this case)------------ k=3, hp(3)
cc      depth_fog=z_RH-hsfc              |   
cc         | ............ o-o (10m)      |
cc         |   ^  ....[2m] |            z_RH=hp(kfog)
cc      ___v___|____^___Y__|_____________|___________________________ surface level (hsfc)
cc     /       |    |                    |                           \
cc-------------|----|--------------------|---kp (=2 in this case) --------------- k=2, hp(2)    
cc   /         |   z2m=2+hsfc            |                             \
cc  /          |    |                    |                              \          
cc /        z10m=10+hsfc                 |                               \
cc ------------|----|--------------------|-----------------------------------------k=1, hp(1)
cc/            |    |                    |                                 \
cc ____________v____v____________________v__________________________________\______sea level
cc|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
cc
cc
cc  (2) Stratus cloud case (cloud top: cldt<400m)
cc
cc  ~~~~~~~cloud top (cldt < 400m)~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 
cc         ^                             ^
cc         |                             |
cc   -kfog-|- (level just below cldt) ---|---------------------------------------k=4. hp(4)
cc         |                             |
cc         |                             |
cc      depth_fog=cldt (i.e. z_RH-hsfc)  |
cc         |                             |
cc  -------|-----------------------------|----------------------------------------k=3, hp(3)
cc         |                             |
cc         | ............ o-o (10m)      |
cc         |   ^  ....[2m] |            z_RH=cldt+hsfc
cc      ___v___|____^___Y__|_____________|___________________________ surface level (hsfc)
cc     /       |    |                    |                           \
cc-------------|----|--------------------|---kp (=2 in this case) --------------- k=2, hp(2)        
cc   /         |   z2m=2+hsfc            |                             \
cc  /          |    |                    |                              \
cc /        z10m=10+hsfc                 |                               \
cc ------------|----|--------------------|-----------------------------------------k=1, hp(1)
cc/            |    |                    |                                 \
cc ____________v____v____________________v__________________________________\______sea level
cc||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c

        subroutine get_new_fog(up,vp,hp,u10,v10,hsfc,rhp,
     +         rh2,t2,tp,time_step,lv,cldt,cldb,adv,rhr,
     +             lwc)

         parameter (Rd = 287.0, Rv=416.0)

c                    Input
c         up,vp : u,v on pressure levels
c            hp : height of pressure levels
c       u10,v10 : 10m u and v
c          hsfc : sfc height
c            tp : previous and current temperature on  pressure levels
c                 tp(1:) presious, tp(2:) current temerature
c           rhp : RH on pressure levels
c            t2 : previous and current temperature at 2m
c                 t2(1) previous, t2(2) current temerature
c           rh2 : previous and current RH at 2m
c          cldt : cloud top height above the ground
c          cldb : cloud base height above the ground
c     time_step : time step interval (in second) between current and previous temperature    
c            lv : number of verical levels 
c           adv : advection term of LWC
c
c                    Output
c           lwc : fog LWC
c
c    Note: lv=14 is set for NCEP regional models, you can set others for observational data
c          or other models
c
c



        REAL up(lv),vp(lv),hp(lv),rhp(lv),tp(lv,2)
        REAL u10,v10,hsfc,lwc,rh2,t2(2),cooling_rate,
     +       Km,lamda_1,lamda,lwc_max,Kc,cldt,cldb

        real rhthr

         if(t2(2).lt.-20.0) then                   !ice fog case
            ttw=17.62*t2(2)/(243.12 + t2(2))     !over water
            tti=22.46*t2(2)/(272.62 + t2(2))     !over ice
                                                 !esi = 6.112 *2.71828 ** tti
                                                 !esw = 6.112 *2.71828 ** ttw
            rhthr = (2.71828 ** (tti-ttw))*100.0   !rhthr = esi / esw, ice fog saturation threshold
         else
          rhthr=rhr
         end if

        write(*,*) 'rhthr=',rhthr
        write(*,*) 'Input data:'
        do k=lv,1,-1
         write(*,30) hp(k),tp(k,2),rhp(k),up(k),
     +      vp(k),tp(k,1)
        end do
        write(*,30) hsfc,t2(2),rh2,u10,v10,t2(1)
        write(*,'(2f10.1)') cldt,cldb
        write(*,'(e10.4)')  adv 

30     format(5x,f6.1, 5f5.1)

        write(*,*) 'Output:'

        lwc = 0.0
c   search sfc in which pressure layers     
c         hp(kp)   : the  level just below the surface level
c         hp(kp+1) : the  level just above the surface level

         z10m = 10. + hsfc
         z2m  =  2. + hsfc

         write(*,*)'z10m=',z10m,' z2m=',z2m
          if (z10m .lt. hp(1) ) then
            kp = 0
          else
            do k = 1,lv-1
             if(z10m.ge.hp(k).and.z10m.lt.hp(k+1)) then
               kp = k
             end if
            end do
          end if

          write(*,*)'kp=',kp  
      
cc
cc  Now see if there is low stratus cloud near the surface
cc         
          if(cldt.lt.400.0.and.cldb.lt.50.0) then    !yes there is stratus cloud case
            kfog=0                                   !400, 50 can be re-tuned for low stratus 
            z_RH=cldt+hsfc
            !seach which level the cloud top is below
            do k = 1, lv-1
             if(z_RH.ge.hp(k).and.z_RH.lt.hp(k+1)) then
               kfog = k+1
             end if
            end do
           goto 1000                           
          end if

c   search RH=100% (95%) depth above the sfc and averaged cooling rate
c         hp(kfog) : the  level just above the fog top

          z_RH=hsfc  
          kfog=0
          do k = 1, lv-1
            if(rhp(k).ge.rhthr.and.
     +         rhp(k+1).lt.rhthr ) then
              kfog=k           !find which pressure level is saturated
              z_RH = hp(kfog)  
            end if
          end do

           write(*,*) 'z_RH-hsfc=',z_RH-hsfc
          if((z_RH-hsfc).gt.800.0) return            !not a fog layer but a stratus cloud

          !now search from fog top downward to find fog bottom reach the ground?
          kbottom=0
          if(kfog.ge.2) then
            do k = kfog,1,-1
             if (rhp(k).lt.rhthr) then
              kbottom=k+1
              goto 100
             end if
            end do
           else
           kbottom=1
           goto 100
          end if

100       continue
          if(rh2.gt.rhthr) then
            touch_down=0.
          else
            touch_down=hp(kbottom)-hsfc
          end if

          write(*,*) 'touch_down=',touch_down
          if(touch_down.gt.100.0) return               !not a fog layer but a stratus cloud 
     
1000      continue
 
          write(*,*)'z_RH=',z_RH,' kfog=',kfog,
     +    ' rh2=',rh2,' z2m=',z2m

          if(z_RH.le.z2m) then
            if(rh2.ge.rhthr) then 
              depth_fog = 2.0
              cooling_rate=(t2(1)-t2(2))/time_step       !in C/sec 
              tfog=t2(2)
c              if (cooling_rate.le.0.0) return 
            else
              depth_fog = 0.0
              return
            end if 
          else
             depth_fog = z_RH - hsfc
             cooling_rate=( (tp(kfog,1)-tp(kfog,2)) +     !average cooling in fog layer
     +         (t2(1)-t2(2)) )/time_step/2.0
             tfog=(tp(kfog,2)+t2(2))/2.0                  !average temperature in fog layer
          end if



          write(*,*)'depth_fog=',depth_fog,
     +      ' cooling_rate=',cooling_rate,
     +      ' tfog=',tfog

C    es:    Use WMO recommended formulation (2008)
C   Guide to Metteorological Instruments and Methods of
C   Observation (CIMO Guide):

          if( tfog.ge.0.0) then
            tt=17.62*tfog/(243.12 + tfog)     !over liquid water
          else
            tt=22.46*tfog/(272.62 + tfog)     !over ice, assume no supercooled water
          end if
          es = 6.112 *2.71828 ** tt                     !es:saturated vapor pressure 

          Tk=tfog + 273.17
          beta = 622.0*2.5e6*es/(Rv*Tk*Tk*1000)
          write(*,*) 'beta=',beta
          ct_term =beta*cooling_rate                !in g/kg/sec 
          adv_term = adv                            !in g/kg/sec
          write(*,*) ' cooling_term=',ct_term,
     +               ' adv_term=',adv_term
          c_adv= ct_term + adv_term               !total LWC's time increment due to cooling and advection

          if(c_adv.lt.0.0) then
            lwc=0.0
            return
          end if 
          lwc_max = sqrt(c_adv*depth_fog/0.062)

          write(*,*)'es=',es,' Tk=',Tk,
     +      ' beta=',beta,' lwc_max=',lwc_max

c   compute FBL (Fog Boundary Layer, delta)
          if ( depth_fog .le. 2.0 ) then
            wp=sqrt(up(kp+1)*up(kp+1)+vp(kp+1)*vp(kp+1))
            w10=sqrt(u10*u10+v10*v10)
            dwdz=(wp-w10)/(hp(kp+1)-z10m)
            dtdz=(tp(kp+1,2)-t2(2))/(hp(kp+1)-z2m)
          else
            wp=sqrt(up(kfog)*up(kfog)+vp(kfog)*vp(kfog))
            w10=sqrt(u10*u10+v10*v10)
            dwdz=(wp-w10)/(hp(kfog)-z10m)
            dtdz=(tp(kfog,2)-t2(2))/(hp(kfog)-z2m)
          end if
          
          dwdz=abs(dwdz)
          if(dwdz.eq.0.0) dwdz=0.001
          write(*,*)'wp=',wp,' w10=',w10,' dwdz=',dwdz,
     +       ' dTdz=',dtdz  

          Ri = (9.8/Tk)*(dtdz+0.01)/dwdz/dwdz 
          if (Ri.le.0.0) then                        !Unstable case: Beljaars: "The parameterization 
           z=hp(kp+1)-hsfc                           !of the planetary boundary layer", May 1992, ECMWF report
           y=(1.0+z/0.02)                            !You can use other formulation
           Ch=12.0*sqrt(y)/(log(y)*log(y))
           Fm=1.0-10.0*Ri/(1.+Ch*sqrt(-Ri))
          else                                       !Stable case: Holtslag, et al. BLM 2006, vol. 118
c           Fm=1.0/(1.0+10.0*Ri)                     !long tail format, Ric may > 0.25 this leads to higher turbulence
                                                     !and fog LWC becomes smaller or dispersed
                                                     !You can use other formulation

            if(Ri.ge.0.0.and.Ri.lt.0.1) then         !Sharp tail (Ric ~ 0.25)
             Fm=(1.0-5.0*Ri)*(1.0-5.0*Ri)
            else
             Fm=(1.0/20.0/Ri)*(1.0/20.0/Ri)
            end if
          end if

           write(*,*)'Ri=',Ri,' FM=',FM

         
         lamda_1=1./(0.4*(depth_fog+0.02)) + 1./40.0
         lamda=1.0/lamda_1
          Km=lamda*lamda*dwdz*Fm
          FBL=Km/2.0/
     +     sqrt(0.062*c_adv*depth_fog)

          Kc=1.38*sqrt(0.062*c_adv)*
     +        depth_fog**1.5

          write(*,*)'Km=',Km,' Kc=',Kc,
     +     ' lamda=',lamda,' FBL=',FBL

c  at last compute LWC near the surface (compute LWC  at 1m for 2m fog, at 0.2 of depth of other fog)
         if(depth_fog.le.2.0) then
            z=1.0
         else
            z=0.1*depth_fog
         end if
         
         ax=sqrt(1.0-z/depth_fog)
         cx=2.0/(1.0+exp(z/FBL))

         lwc=lwc_max*(ax - cx)                            

 
         write(*,*)'LWC profile:'         
         dz=depth_fog/10.0
         do k=0,10
          z=dz*k
          prof_lwc=lwc_max*(sqrt(1.0-z/depth_fog) -       
     +       2.0/(1.0+exp(z/FBL)))
          if (prof_lwc.lt.1.0E-2) prof_lwc=0.0
          write(*,*) z, 'm: LWC=',prof_lwc
         end do

         if (lwc.lt.1.0E-2) lwc = 0.0

       return
       end