anchor.diff:          real,dimension(:,:),intent(inout) :: exchange_C
   
           ! Richardson number
  -        where(wind==0) wind=1e-10
  +        where(wind==0) wind=1e-5
  +        exchange_C = 0
   
           Ri = gravity/airt(:,1,:) * (airt(:,1,:)-tskin)*z_atm/(wind**2)
   
  -!         print*,"--------------------------------------------------"
  -!         print*, "Surface Richardson number"
  +        ! print*,"--------------------------------------------------"
  +        ! print*, "Surface Richardson number"
           where(Ri<0)  exchange_C = lnz_atm_term * (1.0-(15.0*Ri)/(1.0+(base_exchange_term * sqrt((-1.0)*Ri))))
           where(Ri>=0) exchange_C = lnz_atm_term * 1.0/((1.0+15.0*Ri)*sqrt(1.0+5.0*Ri))
   
  
positive.code: 
              endif
  
              ! accumulate soil moisture over the entire column
              domain%soil_totalmoisture%data_2d = domain%soil_water_content%data_3d(:,1,:) * DZS(1) * 1000
              do i = 2,num_soil_layers
                  domain%soil_totalmoisture%data_2d = domain%soil_totalmoisture%data_2d + domain%soil_water_content%data_3d(:,i,:) * DZS(i)
              enddo
  
  
              ! 2m Air T and Q are not well defined if Tskin is not coupled with the surface fluxes
              if (options%physics%landsurface > kLSM_BASIC) then
                  call surface_diagnostics(domain%sensible_heat%data_2d,    &
                                           QFX,                             &
                                           domain%skin_temperature%data_2d, &
  
negative.code: 
  
          ! initial guesses (not needed?)
          domain%temperature_2m%data_2d = 280
          domain%humidity_2m%data_2d = 0.001  ! domain%water_vapor%data_3d(:,1,:)
  
          if (options%physics%landsurface==kLSM_SIMPLE) then
              write(*,*) "    Simple LSM (may not work?)"
  