anchor.diff:              elseif (exchange_term==2) then
                   call calc_mahrt_holtslag_exchange_coefficient(windspd,domain%skin_temperature%data_2d,domain%temperature%data_3d,domain%roughness_z0%data_2d,CHS)
               endif
  -!             print*, CHS(128,103)
   
               ! --------------------------------------------------
               ! First handle the open water surface options
  
positive.code:         real,dimension(:,:),intent(inout) :: exchange_C
  
          ! Richardson number
          where(wind==0) wind=1e-10
  
          Ri = gravity/airt(:,1,:) * (airt(:,1,:)-tskin)*z_atm/(wind**2)
  
  !         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))
  
  
negative.code:             lsm_dt = domain%model_time%seconds() - last_model_time
              last_model_time = domain%model_time%seconds()
  
              domain%temperature%data_3d = domain%exner%data_3d * domain%potential_temperature%data_3d
  
              ! exchange coefficients
              windspd = sqrt(domain%u_10m%data_2d**2 + domain%v_10m%data_2d**2)
              if (exchange_term==1) then
  