ufo_coolskin_sim_mod.F90

Go to the documentation of this file.

module ufo_coolskin_sim_mod
  use kinds
  use ufo_constants_mod 
  implicit none
  public :: ufo_coolskin_sim, ufo_coolskin_jac
  private
 
contains
  subroutine ufo_coolskin_sim(Ts,dTc,S_ns,H_I,H_s,R_nl,Tdc,u0)
  use kinds 
  implicit none
 
  real(kind=kind_real),        intent(inout) :: ts,dtc      !dTc to return
  real(kind=kind_real),        intent(in)    :: s_ns,h_i,h_s,r_nl,u0,tdc
  
  ! local variables
  integer:: n_i,i
  real(kind=kind_real) :: delta,fc, u, lamda,q0,qb,td
  
 
  u     = max(0.0002, u0) !friction velocity over water
  dtc   = 0.0
  n_i   = 3
  td = tdc + 273.15 ! convert from C to K
 
  do i = 1,n_i
      q0     = h_i + h_s + (eps * sig * (td - dtc)**4 - r_nl)
      qb     = q0 + (s_b*cw/(alpha*l_e))*h_i
      lamda = 6.0*(1.0+((alpha*gr*qb/(rou*cw))*(16.0*rou**2.0 * cw**2.0 * v_w**3.0 / &
      k_t**2.0)* (1/u**4.0))**(3.0/4.0))**(-1.0/3.0)
      delta  =lamda * v_w / u
      fc  = 0.0685 + 11.0 * delta - 3.3e-5 /delta * (1.0-exp(-delta/(8.0e-4)))
      dtc  = (h_i + h_s + (eps * sig * (td - dtc)**4 - r_nl) - s_ns * fc) * delta/k_t
      dtc = max(dtc,0.0)
  enddo
 
  ts = td - dtc
  ts = ts - 273.15 ! convert from K to C
  
  end subroutine ufo_coolskin_sim
  
  !-----------------------------------------------------------------------
  
  subroutine ufo_coolskin_jac(jac,S_ns,H_I,H_s,R_nl,Tdc,u0)
  use kinds
  implicit none
  
  real(kind=kind_real),        intent(in)    :: s_ns,h_i,h_s,r_nl,u0,tdc
  real(kind=kind_real),        intent(out)   :: jac(6) ! jac calculated for all inputs
 
  real(kind=kind_real) :: delta ,fc ,u, td
  real(kind=kind_real) :: lamda ,q0 ,qb ,ts ,dtc ,c0 ,y ,q
  real(kind=kind_real) :: const ,d_lamda_dqb, dfc_d_delta
 
  u     = max(0.0002, u0) !friction velocity over water
 
  call ufo_coolskin_sim(ts,dtc,s_ns,h_i,h_s,r_nl,tdc,u)
  
  td = tdc + 273.15  ! convert from C to K
  ts = ts + 273.15
  q0      = h_i + h_s + (eps * sig * ts**4 - r_nl)
   
  !constant apears in net heat equation
  c0      = s_b*cw/(alpha*l_e)   
  qb      = q0 + c0*h_i
 
  !Saunder’s constant
  lamda   = 6.0*(1.0+((alpha*gr*qb/(rou*cw))*(16.0*rou**2.0 * cw**2.0 * v_w**3.0 &
            / k_t**2.0)* (1/u**4.0))**(0.75))**(-1.0/3.0)
  
  ! cool layer thickness
  delta   = lamda * v_w / u
  
  !solar absorbption profile in cool layer
  fc      = 0.0685 + 11.0 * delta - 3.3e-5 /delta * (1.0-exp(-delta/(8.0e-4)))
  
  ! net heat in cool layer
  q       = h_i + h_s + (eps * sig * (td - dtc)**4 - r_nl) - s_ns * fc
  const   = ((alpha*gr/(rou*cw))*(16.0*rou**2.0 * cw**2.0 * v_w**3.0 / k_t**2.0))**(0.75)
  
  ! calculate d(fc)/d(delta)
  dfc_d_delta  = 11 + 3.3e-5 / delta**2 *(1.0-exp(-delta/(8.0e-4))) - 3.3e-5 &
  / delta * (1.0/(8.0e-4) * exp(-delta/(8.0e-4)))
  
  ! calculate d(lamda)/d(Qb)
  d_lamda_dqb  = -2.0 *(1.0+const *(qb/u**4) **(0.75))**(-4.0/3.0) * &
                 (0.75) * const * (qb/u**4) ** (-0.25)
  
  ! this apears in several of jacobians, I decided to calculate it once (4 eps * sig * Ts**3)
  y            = 4 * eps * sig * ts**3
 
  ! d(Ts)/d(S_ns)
  jac(1) = fc * delta /(k_t+y*(delta+(q - s_ns*delta *dfc_d_delta)*v_w/u*d_lamda_dqb))
  
  ! d(Ts)/d(H_I)
  jac(2) = -((1+c0)*(delta+q*v_w/u*d_lamda_dqb-s_ns*dfc_d_delta*v_w/u*d_lamda_dqb))/ &
    (y*(delta+q*v_w/u*d_lamda_dqb-s_ns*dfc_d_delta*v_w/u*d_lamda_dqb)+k_t)
  
  ! d(Ts)/d(H_s)
  jac(3) = -((delta+q*v_w/u*d_lamda_dqb-s_ns*dfc_d_delta*v_w/u*d_lamda_dqb))/ &
    (y*(delta+q*v_w/u*d_lamda_dqb -s_ns*dfc_d_delta*v_w/u*d_lamda_dqb)+k_t)
  
  ! d(Ts)/d(R_nl)
  jac(4) = (delta +(q-s_ns*dfc_d_delta*delta)*(v_w/u*d_lamda_dqb))/ &
    (k_t+(q-s_ns*dfc_d_delta*delta)*(y*v_w/u*d_lamda_dqb))
  
  ! d(Ts)/d(Td)
  jac(5) = k_t/(k_t + y *(delta - delta*s_ns*dfc_d_delta +q) * v_w/u * d_lamda_dqb)
    
  ! d(Ts)/d(u)
  jac(6) = - (q-s_ns* delta*dfc_d_delta)*(v_w/u**2.0 *lamda + (4.0*u**(-6.0))*d_lamda_dqb*qb)/ &
    (k_t+y *(delta +(q-delta*s_ns*dfc_d_delta)*(v_w*d_lamda_dqb*u**(-5.0))))
  
  end subroutine ufo_coolskin_jac
  
  
end module ufo_coolskin_sim_mod