!    ==========================================================
program THC_4box_with_seaice
!  ==========================================================
  implicit none
  include 'THC_4box_with_seaice.h'
  real*8 day, emp_flux(2)
  integer iocn_seaice,icoupling_period

  call open_files
  call get_parameters

  !  ---------------------------
  !  step coupled model in time:
  !  ---------------------------
  do icoupling_period = 1,ncoupling_periods
     call air_sea_fluxes
     ! step atmospheric model for one coupling_period
     if (run_mode.eq.'cpl') call atmosphere
     ! step ocean and seaice for one coupling_period
     do iocn_seaice=1,ocn_seaice_steps_per_coupling
        ! step ocean one time step:
        call sea_ice
        ! step seaice one time step:
        call ocean
     end do
  end do

  call write_restart

  stop
end Program THC_4box_with_seaice


subroutine ocean
  !  Do a ONE TIME STEP of 
  !  a 4-box thermohaline circulation model.
  !
  !       -------------     -----                 	}
  !       | 1         |---> | 2 |                 	} upper boxes
  !       |           |     |   |                 	}
  !       -------------     -----                 	}
  !        ^                              
  !        |                  |           
  !        |                  |           
  !        |                  V           
  !       -------------     -----                 	}
  !       | 4         | <-- | 3 |                 	}
  !       |           |     |   |                 	}  deep boxes
  !       -------------     -----                 	}
  !       EQ                 NP

  implicit none
  include 'THC_4box_with_seaice.h'
  integer nt,i,m
  real*8 f12(1:ntracers_ocn),f23(1:ntracers_ocn),f34(1:ntracers_ocn),f41(1:ntracers_ocn)
  real*8 sigma,P_sigma,T_sigma,S_sigma
  real*8 rhoUpperNorth, rhoUpperSouth, rhoLowerNorth, rhoLowerSouth
  real*8 hdiff12(1:ntracers_ocn), hdiff34(1:ntracers_ocn) &
       , vdiff14(1:ntracers_ocn), vdiff23(1:ntracers_ocn)
  real*8 kappa_H_srfc1
  
  !     ==================================
  !     Step ocean model ONE STEP in time:
  !     ==================================
  do nt=ntime_steps_ocn+1,ntime_steps_ocn+1

     !        Time since begining of model run:
     time=dt_ocn*(nt-1)

     !        ---------
     !        Calc rho:
     !        ---------
     if (linear_eqn_state) then
        rho(1:nboxes_ocn)= &
        rho_0*(1+alpha*(T(2,1:nboxes_ocn,1)-10)+beta*(T(2,1:nboxes_ocn,2)-35))
     else
        do i=1,nboxes_ocn
           P_sigma=0.0
           T_sigma=T(2,i,1)
           S_sigma=T(2,i,2)
           rho(i)=sigma(P_sigma,T_sigma,S_sigma)
        end do
     end if
 
     !        -----------------------------------
     !        Calc mass transports between boxes:
     !        -----------------------------------

     rhoUpperSouth=rho(1)
     rhoUpperNorth=rho(2)
     rhoLowerNorth=rho(3)
     rhoLowerSouth=rho(4)

     !        ---------------------------------------------
     !        Check static stability and convect if needed:
     !        ---------------------------------------------

     if (rhoUpperNorth.gt.rhoLowerNorth) then
        ! Convect:
        T(2,2,1:ntracers_ocn)=(T(2,2,1:ntracers_ocn)*V_ocn(2) + T(2,3,1:ntracers_ocn) &
             *V_ocn(3))/(V_ocn(2)+V_ocn(3))
        T(2,3,1:ntracers_ocn)=T(2,2,1:ntracers_ocn)
     end if

     if (rhoUpperSouth.gt.rhoLowerSouth) then
        ! Convect:
        T(2,1,1:ntracers_ocn)=(T(2,1,1:ntracers_ocn)*V_ocn(1) + T(2,4,1:ntracers_ocn) &
             *V_ocn(4))/(V_ocn(1)+V_ocn(4))
        T(2,4,1:ntracers_ocn)=T(2,1,1:ntracers_ocn)
     end if

     ! Recalculate density:
     if (rhoUpperNorth.gt.rhoLowerNorth  &
          .or. &
          rhoUpperSouth.gt.rhoLowerSouth) then
        if (linear_eqn_state) then
           rho(1:nboxes_ocn)= &
                rho_0*(1+alpha*(T(2,1:nboxes_ocn,1)-10)+beta*(T(2,1:nboxes_ocn,2)-35))
        else
           do i=1,nboxes_ocn
              P_sigma=0.0
              T_sigma=T(2,i,1)
              S_sigma=T(2,i,2)
              rho(i)=sigma(P_sigma,T_sigma,S_sigma)
           end do
        end if
        rhoUpperSouth=rho(1)
        rhoUpperNorth=rho(2)
        rhoLowerNorth=rho(3)
        rhoLowerSouth=rho(4)
     end if

     u_thermohaline=C_u*g*0.25*( &
          DepthUpper_ocn*(rhoUpperNorth-rhoUpperSouth)+ &
          DepthLower_ocn*(rhoLowerNorth-rhoLowerSouth) &
          )

     if (stable_stratification_stop_THC) then
        ! Thermally-dominant THC is active only if sinking
        ! area density is unstable:
        if (u_thermohaline.gt.0.0) then
!           u_thermohaline=u_thermohaline &
!                * 0.5*(1+tanh((rhoUpperNorth-rhoLowerNorth)/delta_rho_shutoff))
           if (rhoUpperNorth < rhoLowerNorth-delta_rho_shutoff) then
              u_thermohaline=0.0
           end if
        elseif (u_thermohaline.lt.0.0) then
           if (rhoUpperSouth < rhoLowerSouth-delta_rho_shutoff) then
              u_thermohaline=0.0
           end if
        end if
     end if
     
     ! don't allow for salinity dominant THC:
     if (no_reversed_THC.and.u_thermohaline.lt.0.0) then
        u_thermohaline=0.0
     end if

! calculate mass fluxes that go into the advection of T,S:
     u12=u_thermohaline
     u23=u_thermohaline
     u34=u_thermohaline
     u41=u_thermohaline


     !        -----------------------------------
     !        Calc heat/salt transports between boxes:
     !        -----------------------------------

     if (u12.ge.0.0) then
        f12(1:ntracers_ocn)=u12*T(2,1,1:ntracers_ocn)*[rho_0*Cp_ocn, 1.0d0]
     else
        f12(1:ntracers_ocn)=u12*T(2,2,1:ntracers_ocn)*[rho_0*Cp_ocn, 1.0d0]
     end if
     if (u23.ge.0.0) then
        f23(1:ntracers_ocn)=u23*T(2,2,1:ntracers_ocn)*[rho_0*Cp_ocn, 1.0d0]
     else
        f23(1:ntracers_ocn)=u23*T(2,3,1:ntracers_ocn)*[rho_0*Cp_ocn, 1.0d0]
     end if
     if (u34.ge.0.0) then
        f34(1:ntracers_ocn)=u34*T(2,3,1:ntracers_ocn)*[rho_0*Cp_ocn, 1.0d0]
     else
        f34(1:ntracers_ocn)=u34*T(2,4,1:ntracers_ocn)*[rho_0*Cp_ocn, 1.0d0]
     end if
     if (u41.ge.0.0) then
        f41(1:ntracers_ocn)=u41*T(2,4,1:ntracers_ocn)*[rho_0*Cp_ocn, 1.0d0]
     else
        f41(1:ntracers_ocn)=u41*T(2,1,1:ntracers_ocn)*[rho_0*Cp_ocn, 1.0d0]
     end if

     !XX: set to .true. for hysteresis run:
     if (.false.) then
        if (time<tend/2) then
           kappa_H_srfc1=3d3+97d3*time/(tend/2-tbegin)
        else
           kappa_H_srfc1=100d3-97d3*(time-tend/2)/(tend/2-tbegin)
        end if
     else
        kappa_H_srfc1=kappa_H_srfc
     end if
     hdiff12(1:ntracers_ocn)=kappa_H_srfc1 *(T(1,1,1:ntracers_ocn)-T(1,2,1:ntracers_ocn)) & 
          *(Width_ocn*DepthUpper_ocn) *[rho_0*Cp_ocn, 1.0d0]/(0.5*(L_ocn(1)+L_ocn(2)))
     hdiff34(1:ntracers_ocn)=kappa_H_deep *(T(1,3,1:ntracers_ocn)-T(1,4,1:ntracers_ocn)) & 
          *(Width_ocn*DepthLower_ocn) *[rho_0*Cp_ocn, 1.0d0]/(0.5*(L_ocn(3)+L_ocn(4)))
     vdiff14(1:ntracers_ocn)=kappa_V_south*(T(1,1,1:ntracers_ocn)-T(1,4,1:ntracers_ocn)) & 
          *(Width_ocn*LengthSouth_ocn)*[rho_0*Cp_ocn, 1.0d0]/(0.5*(D_ocn(1)+D_ocn(4)))
     vdiff23(1:ntracers_ocn)=kappa_V_north*(T(1,2,1:ntracers_ocn)-T(1,3,1:ntracers_ocn)) & 
          *(Width_ocn*LengthNorth_ocn)*[rho_0*Cp_ocn, 1.0d0]/(0.5*(D_ocn(2)+D_ocn(3)))

     !   ---------
     !   step T,S:
     !   ---------
     T(3,1,1:ntracers_ocn)=T(1,1,1:ntracers_ocn)+2.*dt_ocn*( &
          -f12(1:ntracers_ocn) &
          +f41(1:ntracers_ocn) &
          +air_to_sea_flux(1,1:ntracers_ocn)*[1.0d0,-S0]*A_ocn(1) &
          -ocean_to_seaice_flux(1,1:ntracers_ocn)*[1.0d0,-S0]*A_ocn(1) &
          -vdiff14(1:ntracers_ocn) &
          -hdiff12(1:ntracers_ocn) &
          )/(V_ocn(1)*[rho_0*Cp_ocn, 1.0d0])

     T(3,2,1:ntracers_ocn)=T(1,2,1:ntracers_ocn)+2.*dt_ocn*( &
          +f12(1:ntracers_ocn) &
          -f23(1:ntracers_ocn) &
          +air_to_sea_flux(2,1:ntracers_ocn)*[1.0d0,-S0]*A_ocn(2) &
          -ocean_to_seaice_flux(2,1:ntracers_ocn)*[1.0d0,-S0]*A_ocn(2) &
          -vdiff23(1:ntracers_ocn) &
          +hdiff12(1:ntracers_ocn) &
          )/(V_ocn(2)*[rho_0*Cp_ocn, 1.0d0])

     T(3,3,1:ntracers_ocn)=T(1,3,1:ntracers_ocn)+2.*dt_ocn*( &
          +f23(1:ntracers_ocn) &
          -f34(1:ntracers_ocn) &
          +vdiff23(1:ntracers_ocn) &
          -hdiff34(1:ntracers_ocn) &
          )/(V_ocn(3)*[rho_0*Cp_ocn, 1.0d0])

     T(3,4,1:ntracers_ocn)=T(1,4,1:ntracers_ocn)+2.*dt_ocn*( &
          +f34(1:ntracers_ocn) &
          -f41(1:ntracers_ocn) &
          +vdiff14(1:ntracers_ocn) &
          +hdiff34(1:ntracers_ocn) &
          )/(V_ocn(4)*[rho_0*Cp_ocn, 1.0d0])

     !   ---------------------------
     !   Prepare for next time step:
     !   ---------------------------
     ! Robert filter:
     T(2,1:nboxes_ocn,1:ntracers_ocn) = T(2,1:nboxes_ocn,1:ntracers_ocn) &
          + robert_filter_coeff * &
          (   T(3,1:nboxes_ocn,1:ntracers_ocn) &
          - 2*T(2,1:nboxes_ocn,1:ntracers_ocn) &
          +   T(1,1:nboxes_ocn,1:ntracers_ocn) &
          )

     !  normal time step:
     T(1,1:nboxes_ocn,1:ntracers_ocn)=T(2,1:nboxes_ocn,1:ntracers_ocn)
     T(2,1:nboxes_ocn,1:ntracers_ocn)=T(3,1:nboxes_ocn,1:ntracers_ocn)

     !        ------------------
     !        Print diagnostics:
     !        ------------------
     if ((nt.eq.1).or. &
          ( abs(time - int(time/tprint_short)*tprint_short) &
          .lt. dt_ocn/2.0)) then
        write(18,'(100g14.6)') time/(3600.*24.*365.) &
             ,u_thermohaline/Sv &
             ,(T(3,i,1),i=1,nboxes_ocn) &
             ,(T(3,i,2),i=1,nboxes_ocn) &
             ,(rho(i),i=1,nboxes_ocn) &
             ,(frac_seaice(i),i=1,nsurface_boxes_ocn) &
             ,(air_to_sea_flux(i,1),i=1,nsurface_boxes_ocn) &
             ,(air_to_sea_flux(i,2),i=1,nsurface_boxes_ocn)
     end if

     if ((nt.eq.1).or. &
          abs(time-int(time/tprint)*tprint) &
          .lt.dt_ocn/2.0) then
        write(16,'(" ",/," ==>> Time (ocean, years)=",g14.6)') time/(3600.*24.*365.)
        write(16,'('' nt='',i12,/ &
             &,''; u_thermohaline='',g14.6 &
             &,/,''  ; T_t='',4g14.6 &
             &,/,''  ; S_t='',4g14.6 &
             &,/,''  ; T  ='',4g14.6 &
             &,/,''  ; S  ='',4g14.6 &
             &,/,''  ; rho='',4g14.6 &
             &)') &
             nt &
             ,u_thermohaline/Sv&
             ,(T(3,i,1)-T(1,i,1),i=1,nboxes_ocn) &
             ,(T(3,i,2)-T(1,i,2),i=1,nboxes_ocn) &
             ,(T(3,i,1),i=1,nboxes_ocn) &
             ,(T(3,i,2),i=1,nboxes_ocn) &
             ,(rho(i),i=1,nboxes_ocn)
        write(16,'('' H(1,2)='',4g14.6,'' watts/m^2'')') &
             air_to_sea_flux(1,1) &
             ,air_to_sea_flux(2,1)
        write(16,'('' E-P(1,2)='',2g14.6,'' m/year'')') &
             air_to_sea_flux(1,2)*year &
             ,air_to_sea_flux(2,2)*year
        write(16,'(" avg total H=",g14.6,"; avg total E-P=",g14.6)') &
           (air_to_sea_flux(1,1)*A_ocn(1)+air_to_sea_flux(2,1)*A_ocn(2)) &
           /(A_ocn(1)+A_ocn(2)), &
           (air_to_sea_flux(1,2)*A_ocn(1)+air_to_sea_flux(2,2)*A_ocn(2)) &
           *year/(A_ocn(1)+A_ocn(2))
        write(16,'('' ocean_to_seaice_flux(:,1)='',2g14.6,'' watts'')') &
          ocean_to_seaice_flux(1:nsurface_boxes_ocn,1)
        write(16,'('' ocean_to_seaice_flux(:,2)='',2g14.6,'' m/year'')') &
          ocean_to_seaice_flux(1:nsurface_boxes_ocn,2)*year

     end if
  end do
  ntime_steps_ocn=ntime_steps_ocn+1
  return
end subroutine ocean


subroutine write_restart
  implicit none
  integer n,i,m
  include 'THC_4box_with_seaice.h'

  if (write_restart_data) then
     rewind(36)
     !        write restart data for ocean model:
     write(36,*) " T:"
     write(36,*) (((T(n,i,m),n=1,3),i=1,nboxes_ocn),m=1,2)
     write(36,*) " T_star_equiv:"
     write(36,*) (T_star_equiv(i),i=1,nsurface_boxes_ocn)
     write(36,*) " air_to_sea_flux (H, watts/m^2):"
     write(36,*) (air_to_sea_flux(i,1),i=1,nsurface_boxes_ocn)
     write(36,*) " air_to_sea_flux (E-P, m/year):"
     write(36,*) (air_to_sea_flux(i,2)*year,i=1,nsurface_boxes_ocn)
     write(36,*) " transports (Sv):"
     write(36,*) u_thermohaline/Sv
     write(36,*) " atmospheric temperature:"
     write(36,*) ((T_atm(n,i),n=1,3),i=1,nboxes_atm)
     write(36,*) " atmospheric meridional fluxes:"
     write(36,*) f12_heat, f12_water
     write(36,*) " sea ice volume:"
     write(36,*) V_seaice
     write(36,*) " sea ice fraction:"
     write(36,*) frac_seaice
     write(36,*) " ocean_to_seaice_flux:"
     write(36,*) ocean_to_seaice_flux
     write(36,*) " shortwave_to_seaice_heat_flux:"
     write(36,*) shortwave_to_seaice_heat_flux

  end if
  return
end subroutine write_restart


! SIGMA FCN ******* AUG 3 1981 *****
function sigma(P,T,S)
  ! ************************************
  ! DERIVED FROM EOS80
  ! REFERENCES
  ! MILLERO ET AL 1980, DEEP-SEA RES.,27A,255-264
  ! JPOTS NINTH REPORT 1978,TENTH REPORT 1980
  ! UNITS:
  !C      PRESSURE        P        DECIBARS
  !       TEMPERATURE     T        DEG CELSIUS (IPTS-68)
  !       SALINITY        S        NSU (IPSS-78)
  !C                      SIGMA    (10.**-3)G/CM**3
  !
  !C UNITS CHANGED TO DECIBARS AND CGS FOR USE AT MIT
  ! N FOFONOFF REVISED OCT 7 1980
  ! *****************************
  IMPLICIT REAL*8 (A-H,O-Z)
  !*Eli: July 21, 1993: changed sigma and P,T,S to REAL*8
  real*8 sigma
  REAL*8 P,T,S

  REAL*8 K,K0,KW
  ! EQUIV
  EQUIVALENCE (E,D,B1,R4),(BW,B,R3),(C,A1,R2)
  EQUIVALENCE (AW,A,R1,R0),(KW,K0,K)
  !C CHANGE PRESSURE FROM INPUT UNITS OF DECIBARS TO BARS
  DP=DBLE(P)
  DT=DBLE(T)
  DS=DBLE(S)
  DP=DP/10.D0
  ! SQUARE ROOT SALINITY.
  SR = DSQRT(DABS(DS))
  ! DENSITY PURE WATER AT ATM PRESS IN KG/M3 = (10**-3)GM/CM3
  R1 = ((((6.536332D-9*DT-1.120083D-6)*DT+1.001685D-4)*DT &
       -9.095290D-3)*DT+6.793952D-2)*DT- .157406D0
  ! SEAWATER DENSITY ATM PRESS.
  R2 = (((5.3875D-9*DT-8.2467D-7)*DT+7.6438D-5)*DT-4.0899D-3)*DT &
       +8.24493D-1
  R3 = (-1.6546D-6*DT+1.0227D-4)*DT-5.72466D-3
  R4 = 4.8314D-4
  !
  SIGMA = SNGL((R4*DS + R3*SR + R2)*DS + R1)
  !C DENSITY - 1000 AT ATMOSPHERIC PRESSURE
  IF(ABS(P).GE..00001)GO TO 1
  !C SIGMA REMAINS UNCHANGED, SINCE IS (10**-3)GM/CM3
  RETURN
  ! COMPUTE COMPRESSION TERMS
1 E = (9.1697D-10*DT+2.0816D-8)*DT-9.9348D-7
  BW = (5.2787D-8*DT-6.12293D-6)*DT+8.50935D-5
  B = BW + E*DS
  !
  D = 1.91075D-4
  C = (-1.6078D-6*DT-1.0981D-5)*DT+2.2838D-3
  AW = ((-5.77905D-7*DT+1.16092D-4)*DT+1.43713D-3)*DT &
       +3.239908D0
  A = (D*SR + C)*DS + AW
  !
  B1 = (-5.3009D-4*DT+1.6483D-2)*DT+7.944D-2
  A1 = ((-6.1670D-5*DT+1.09987D-2)*DT-0.603459D0)*DT+54.6746D0
  KW = (((-5.155288D-5*DT+1.360477D-2)*DT-2.327105D0)*DT &
       +148.4206D0)*DT+19652.21D0
  K0 = (B1*SR + A1)*DS + KW
  ! EVALUATE PRESSURE POLYNOMIAL AND RETURN
  K = (B*DP + A)*DP + K0
  SIGMA = SNGL((K*DBLE(SIGMA) + 1000.D0*DP)/(K-DP))
  !C SIGMA REMAINS UNCHANGED SINCE IS (10**-3)GM/CM**3
  RETURN
END FUNCTION SIGMA


subroutine open_files
  implicit none
  include 'THC_4box_with_seaice.h'
  open(15,file='THC_4box_with_seaice.data')
  open(16,file='Output/output.THC_4box_with_seaice')
  open(36,file='Output/restart')
  open(18,file='Output/transp.dat')
  return
end subroutine open_files


subroutine get_parameters
  implicit none
  include 'THC_4box_with_seaice.h'
  !     air_to_sea_flux is the one received by the ocean.
  !     flux is the one sent from the atmosphere:
  real*8 t_check,s_check,rho_check,P_check,sigma,my_random_number
  integer n,i,m,my_status

  !      write(*,*) " reading nl..."
  read(15,nl_coupled)
  read(15,nl_ocn)
  read(15,nl_atm)
  read(15,nl_seaice)
  ! 100  continue
  !      write(*,*) " ...done reading nl"
  !      close(15)

  write(16,'('' '')')
  write(16,'(''--------- coupled Input pars: ----------------'')')
  write(16,nl_coupled)
  write(16,nl_ocn)
  write(16,nl_atm)
  write(16,nl_seaice)
  write(16,'(''--------- End coupled input pars.-------------'')')


  ! set atmospheric model time step:
   dt_atm=dt_atm_box

  !     Parameters for the equation of state:
  !     -------------------------------------
  !     values for alpha and beta are from Gill, p. 603, for (T,S)=(10,35);
  !     alpha=\rho^{-1}\partial \rho/\partial T;
  !     beta=\rho^{-1}\partial \rho/\partial S;
  alpha=-1668.e-7
  beta=0.781/(1026.9)
  !     with these coeffs, the dimensional density, in cgs is:
  !     rho = rho_0*{1+alpha*(T-10.)+beta*(S-35.)}
  !     check value:
  T_check=25.
  S_check=37.
  P_check=0.0
  if (linear_eqn_state) then
     rho_check=rho_0*(1+alpha*(T_check-10.)+beta*(S_check-35.))
  else
     rho_check=sigma(P_check,T_check,S_check)
  end if
  write(16,'('' density check value at T,S='',2f4.0 &
       &,'' is '',g12.4)') &
       T_check,S_check,rho_check
  if (linear_eqn_state) then
     write(16,'('' alpha*T check value at T,S='',2f4.0 &
          &,'' is '',g12.4)') &
          T_check,S_check &
          ,rho_0*(  alpha*(T_check-10.)                   )
     write(16,'('' beta*S check value at T,S='',2f4.0 &
          &,'' is '',g12.4)') &
          T_check,S_check &
          ,rho_0*(                      beta*(S_check-35.))
  end if

  !     Initialize surface temp and sal for restoring conditions:
  !     ---------------------------------------------------------
  ! temperature:
  T_star(1,1)=22.0
  T_star(2,1)=0.5
  ! salinity:
  T_star(1,2)=36.
  T_star(2,2)=34.75

  D_ocn(1)=DepthUpper_ocn
  D_ocn(2)=DepthUpper_ocn
  D_ocn(3)=DepthLower_ocn
  D_ocn(4)=DepthLower_ocn
  L_ocn(1)=LengthSouth_ocn
  L_ocn(2)=LengthNorth_ocn
  L_ocn(3)=LengthNorth_ocn
  L_ocn(4)=LengthSouth_ocn
  W_ocn(1)=Width_ocn
  W_ocn(2)=Width_ocn
  W_ocn(3)=Width_ocn
  W_ocn(4)=Width_ocn

  L_atm(1)=LengthSouth_ocn
  L_atm(2)=LengthNorth_ocn
  W_atm(1) = W_ocn(1)/(1-frac_land(1))
  W_atm(2) = W_ocn(2)/(1-frac_land(2))

  !  Ocean box volumes:
  !  ------------------
  V_ocn(1:nboxes_ocn)=L_ocn(1:nboxes_ocn)*D_ocn(1:nboxes_ocn)*W_ocn(1:nboxes_ocn)

  !  area of boxes:
  !  --------------
  A_ocn(1:nboxes_ocn)=L_ocn(1:nboxes_ocn)*W_ocn(1:nboxes_ocn)
  A_atm(1:nboxes_atm)=L_atm(1:nboxes_atm)*W_atm(1:nboxes_atm)

  !  (T,S) initial conditions for the ocean:
  T(1:3,1,1)= 20.0
  T(1:3,2,1)= 1.0
  T(1:3,3,1)= 1.0
  T(1:3,4,1)= 1.0
  T(1:3,1,2)= 36.0
  T(1:3,2,2)= 34.5
  T(1:3,3,2)= 34.5
  T(1:3,4,2)= 34.5

  !  Temperature initial conditions for atm:
  T_atm(1:3,1)= 27.0+273.0
  T_atm(1:3,2)= -5.0+273.0

  !  initial fluxes for the ocean:
  air_to_sea_flux(1:nsurface_boxes_ocn,1:ntracers_ocn)=0.0

  ! initial conditions for seaice:
  V_seaice(:,:) = 0.0
  frac_seaice(:) = 0.0
  ocean_to_seaice_flux(:,:) = 0.0
  shortwave_to_seaice_heat_flux(:) = 0.0
  thickness_seaice(:)=seaice_growth_thickness

  !     initial equivalent atmos temperature:
  t_star_equiv(1:nsurface_boxes_ocn)= 0.0

  if (use_mixed_bc_and_equiv_t .and. .not.read_restart_data) then
     write(*,*) "*** ERROR: must have read_restart=true to use equivalent/mixed b.c."
     stop
  end if
  
  if (read_restart_data) then
     !        Read restart data:
     read(36,*,end=100,iostat=my_status) ! " T:"
100  continue
     if (my_status<0) then
        write(*,*) " *** ERROR: problem reading restart file."
        stop
     end if
     read(36,*) (((T(n,i,m),n=1,3),i=1,nboxes_ocn),m=1,2)
     read(36,*) ! " T_star_equiv:"
     read(36,*) (T_star_equiv(i),i=1,nsurface_boxes_ocn)
     read(36,*) ! " air_to_sea_flux (H, watts/m^2):"
     read(36,*) (air_to_sea_flux_from_restart(i,1),i=1,nsurface_boxes_ocn)
     read(36,*) ! " air_to_sea_flux (E-P, cm/sec):"
     read(36,*) (air_to_sea_flux_from_restart(i,2),i=1,nsurface_boxes_ocn)
     read(36,*) ! " transports (Sv):"
     read(36,*) u_thermohaline_from_restart
     read(36,*) ! " atmospheric temperature:"
     read(36,*) ((T_atm(n,i),n=1,3),i=1,nboxes_atm)
     read(36,*) ! " atmospheric meridional fluxes:"
     read(36,*) f12_heat, f12_water
     read(36,*) ! " sea ice volume:"
     read(36,*) V_seaice
     read(36,*) ! " sea ice fraction:"
     read(36,*) frac_seaice
     read(36,*) ! " ocean_to_seaice_flux:"
     read(36,*) ocean_to_seaice_flux
     read(36,*) ! " shortwave_to_seaice_heat_flux:"
     read(36,*) shortwave_to_seaice_heat_flux

     u_thermohaline_from_restart=u_thermohaline_from_restart*Sv

     !        Save restart values of temperature and salinity:
     T_restart(1:3,1:nboxes_ocn,1:ntracers_ocn)=T(1:3,1:nboxes_ocn,1:ntracers_ocn)

     !        add random noise to restart data, to check stability:
     !        (no noise is added to time level 3 which is used
     !        to calculate the equiv temperature).
     call random_seed()
     do n=1,2
        do i=1,nboxes_ocn
           ! m=1 is temperature, m=2 is salinity:
           call random_number(my_random_number)
           do m=1,ntracers_ocn
              T(n,i,m)=T(n,i,m) &
                   *(1.+amplitude_noise_IC(m)*(2*my_random_number-1))
           end do
        end do
     end do

     !  if not restoring conditions, restore E-P units and multiply
     !  by factor_EmP:
     if (run_mode.eq.'cpl' .or. .not.restoring_conditions) then
        air_to_sea_flux_from_restart(:,2)= &
             air_to_sea_flux_from_restart(:,2)*factor_EmP/year
     end if

     air_to_sea_flux(:,:)=air_to_sea_flux_from_restart(:,:)

  end if

  !        --------------------------
  !        Calculate equivalent temp:
  !        --------------------------
  !        (calculated from time level 3 to which noise is
  !        not added when restart data is read in.)

  if (use_mixed_bc_and_equiv_t) then
     do i=1,nsurface_boxes_ocn
        t_star_equiv(i)= &
             air_to_sea_flux(i,1)/(gamma_equiv*DepthUpper_ocn*Cp_ocn*rho_0) &
             +T(3,i,1)
     end do
     write(16,*) ' Calculated equivalent temperature:'
     write(16,'('' T_star      (1,2)='',2g14.6)') &
          (t_star(i,1),i=1,nsurface_boxes_ocn)
     write(16,'('' T_star_equiv(1,2)='',2g14.6)') &
          (t_star_equiv(i),i=1,nsurface_boxes_ocn)
 end if

  ! initialize time step counters:
  ntime_steps_ocn=0
  ntime_steps_atm=0
  ntime_steps_seaice=0
  ! Calculate number of time steps per coupling period:
  ncoupling_periods=(tend-tbegin)/coupling_period+1
  atm_steps_per_coupling_period=coupling_period/dt_atm
  ocn_seaice_steps_per_coupling=coupling_period/dt_ocn
  if (abs(coupling_period-atm_steps_per_coupling_period*dt_atm).gt.1.e-5) then
     write(*,*) "There should be an integer number of atm steps per coupling period"
     stop
  end if
  if (abs(coupling_period-ocn_seaice_steps_per_coupling*dt_ocn).gt.1.e-5) then
     write(*,*) "There should be an integer number of ocn_seaice steps per coupling period"
     stop
  end if
  if (dt_ocn.ne.dt_seaice) then
     write(*,*) "ocean and seaice time steps should be equal."
     stop
  end if


  return
end subroutine get_parameters


subroutine air_sea_fluxes
  implicit none
  include 'THC_4box_with_seaice.h'
  integer m,i
  real*8 R,sqrt1mR2,my_random_number
  logical first_call
  data first_call/.true./
  save first_call

  if (run_mode.eq.'cpl') then


     ! air sea heat flux: (watts/m^2)
     air_to_sea_flux(1:nsurface_boxes_ocn,1) = & 
          ! sensible:
          C_bulk_air_sea_heat_flux&
          * (T_atm(2,1:nboxes_atm)-273.0-T(2,1:nsurface_boxes_ocn,1)) &
           * ( &
          ! open ocean part:
          (1-frac_seaice(1:nsurface_boxes_ocn)) &
          ! sea ice covered part:
          + frac_seaice(1:nsurface_boxes_ocn)*0.05/(1.7+thickness_seaice(:)) &
          )

     ! air-sea fresh water flux: (m/sec)
     air_to_sea_flux(1:nsurface_boxes_ocn,2)=[-1.0d0, 1.0d0]*f12_water/A_ocn(1:nsurface_boxes_ocn)


  else if (run_mode.eq.'ocn') then


     if (restoring_conditions) then
        !           -------------------------------------------
        !           calc surface flux for restoring conditions:
        !           -------------------------------------------
        ! heat flux (watts/m^2):
        air_to_sea_flux(1:nsurface_boxes_ocn,1)=gamma(1)*DepthUpper_ocn*Cp_ocn*rho_0 &
             *(T_star(1:nsurface_boxes_ocn,1)-T(2,1:nsurface_boxes_ocn,1))
        ! water flux (m/sec):
        air_to_sea_flux(1:nsurface_boxes_ocn,2)=-(gamma(2)*DepthUpper_ocn/S0) &
             *(T_star(1:nsurface_boxes_ocn,2)-T(2,1:nsurface_boxes_ocn,2))
     elseif (use_mixed_bc_and_equiv_t) then
        !           ---------------------------------------------------------
        !           use restoring with equiv gamma and T for temperature, and
        !           flux for the salinity:
        !           ---------------------------------------------------------
        do i=1,nsurface_boxes_ocn
           air_to_sea_flux(i,1)=gamma_equiv*DepthUpper_ocn*Cp_ocn*rho_0 &
                *(T_star_equiv(i)-T(2,i,1))
        end do
     end if

     if (add_noise_to_surface_EmP) then
        !        add random noise
        R=exp(-dt_ocn/(14*86400)) ! correlation time specified here
        sqrt1mR2=sqrt(1-R**2)
        if (first_call) then
           first_call=.false.
           air_to_sea_flux_noise(1:nsurface_boxes_ocn,1:ntracers_ocn)=0
           call random_seed()
        end if
        call random_number(my_random_number)
        ! m=1 is temperature, m=2 is salinity:
        m=2
        air_to_sea_flux_noise(1,m)=R*air_to_sea_flux_noise(1,m) &
             +sqrt1mR2*amplitude_noise_EmP*(2*my_random_number-1)
        ! make sure total E-P added by noise at each time step is zero:
        air_to_sea_flux_noise(2,m)=-air_to_sea_flux_noise(1,m) &
             *(A_ocn(1)/A_ocn(2))
        !write(66,*) my_random_number, air_to_sea_flux_noise(1,m)
        
        do i=1,nsurface_boxes_ocn
           air_to_sea_flux(i,m)=air_to_sea_flux(i,m) &
                +air_to_sea_flux_noise(i,2)
        end do
     end if
  end if

end subroutine air_sea_fluxes


subroutine atmosphere
  !  A 2-box energy balance model.
  !
  !       -------------     -----  
  !       | 1         |<--->| 2 |  
  !       |           |     |   |  
  !       -------------     -----  
  !       EQ                 NP

  implicit none
  include 'THC_4box_with_seaice.h'
  integer nt,i
  real*8 ac1,R_gas_const,Ps_atm
  character*16 atm_printout_units

  !     =======================
  !     Step atm model in time:
  !     =======================
  do nt=ntime_steps_atm+1,ntime_steps_atm+atm_steps_per_coupling_period

     !        Time since begining of model run:
     time = dt_atm*(nt-1)


     !   ---------------------------------------------
     !   Calc heat/ moisture transports between boxes:
     !   ---------------------------------------------
     f12_heat(1) = (K_heat_atm*( T_atm(2,1)-T_atm(2,2) ) &
          /( 0.5*( L_atm(1)+L_atm(2) )))
     f12_heat(2) = - f12_heat(1)
 
     call calc_FW_flux

     albedo(:)= &
            albedo_clouds(:) &
          + (1-albedo_clouds(:)) &
          *     (  frac_land(:) * albedo_land(:) &
          +        (1-frac_land(:)) &
          *         (  frac_seaice(:)*albedo_seaice &
          +            (1-frac_seaice(:))*albedo_ocean &
          ))
     incoming_radiation(:)= &
          ! use total albedo to calculate radiation that remains in atmosphere:
          solar_incoming_radiation(:)*(1-albedo(:) ) &
          ! substract radiation that is absorbed by seaice:
          -frac_seaice(:)*shotrwave_penetration_to_seaice &
          *solar_incoming_radiation(:)*(1-albedo_clouds(:))*(1-frac_land(:))*frac_seaice(:) &
          *(1-albedo_seaice)

     longwave_outgoing(:) = sigma_SB*emissivity(:) &
          *T_atm(2,:)*T_atm(2,:)*T_atm(2,:)*T_atm(2,:)

     ! Gas constant for dry air (J/Kg*K)
     R_gas_const = 287.04
     ! surface pressure taken as 1000mb
     Ps_atm = 1.0d+5
     ! normalization factor from J/s into K/s:
     ac1=(2.**(R_gas_const/Cp_atm))*g/(Ps_atm*Cp_atm)     ! 2^(R/Cp)*g/(Ps*Cp)

     !     ---------
     !     step T:
     !     ---------
     T_atm(3,:)=T_atm(1,:)+2.*dt_atm*( &
          -f12_heat(:)/A_atm(:) &
          -air_to_sea_flux(:,1)*A_ocn(1:nsurface_boxes_ocn)/A_atm(:) &
          +incoming_radiation(:) &
          -longwave_outgoing(:) &
          )*ac1

     !    ---------------------------
     !    Prepare for next time step:
     !    ---------------------------
        ! Robert filter:
     T_atm(2,1:nboxes_atm) = T_atm(2,1:nboxes_atm) &
          + robert_filter_coeff * &
          (   T_atm(3,1:nboxes_atm) &
          - 2*T_atm(2,1:nboxes_atm) &
          +   T_atm(1,1:nboxes_atm) &
          )

     !  normal time step:
     T_atm(1,:)=T_atm(2,:)
     T_atm(2,:)=T_atm(3,:)

     !        ------------------
     !        Print diagnostics:
     !        ------------------
     if ((nt.eq.1).or. &
          (abs(time- &
          int(time/tprint_short)*tprint_short) &
          .lt. dt_atm/2.0)) then
     !   write(18,'(50g14.6)') time/(3600.*24.*365.) &
     !        ,u_thermohaline/Sv &
     !        ,(T_atm(3,i),i=1,nboxes_atm)
     end if
     if (print_atm_in_Pwatts) then
        convert_to_Pw(:)=A_atm(:)/1.e15
        atm_printout_units="Pwatts   "
     else
        convert_to_Pw(:)=1.0
        atm_printout_units="watts/m^2"
     endif
     if ((nt.eq.1).or. &
          abs(time-int(time/tprint)*tprint) &
          .lt.dt_atm/2.0) then
        write(16,'(  " ---------------------------------------------------------------")')
        write(16,'(" ",/," ==>> Time (atm, years)=",g14.6)') time/year
        write(16,'('' nt(atm)='',i12 &
             &,/,''  ; T(atm)_t='',2g14.6 &
             &,/,''  ; T(atm)  ='',2g14.6 &
             &,/,''  ; f12_heat(1) (atm, Pwatt) = '',g14.6,''; f12_water (Sv) ='',g14.6 &
             &,/,''  ; air_to_sea_flux(:,1)*A_ocn(:)/A_atm(:) (atm, '',a9,'') ='',2g14.6 &
             &,/,''  ; incoming_radiation(:) (atm, '',a9,'') ='',2g14.6 &
             &,/,''  ; longwave_outgoing(:)  (atm, '',a9,'') ='',2g14.6 &
             &,/,''  ; albedo(:)       ='',2g14.6 &
             &,/,''  ; solar_incoming_radiation(:)       ='',2g14.6 &
             &)') &
             nt &
             ,(T_atm(3,i)-T_atm(1,i),i=1,nboxes_atm) &
             ,(T_atm(3,i)-273,i=1,nboxes_atm) &
             ,f12_heat(1)/1.e15, f12_water/Sv &
             ,atm_printout_units,air_to_sea_flux(:,1)*A_ocn(1:nsurface_boxes_ocn)/A_atm(:)*convert_to_Pw(:) &
             ,atm_printout_units,incoming_radiation(:)*convert_to_Pw(:) &
             ,atm_printout_units,longwave_outgoing(:)*convert_to_Pw(:) &
             ,albedo(:) &
             ,solar_incoming_radiation(:)
     end if
  end do
  ntime_steps_atm=ntime_steps_atm+atm_steps_per_coupling_period

end subroutine atmosphere


subroutine sea_ice
  !  Do a ONE TIME STEP of 
  !  a thermodynamic sea ice box model.
  implicit none
  include 'THC_4box_with_seaice.h'
  integer nt, ibox

  !     ==========================
  !     Step seaice model in time:
  !     ==========================
  do nt=ntime_steps_seaice+1,ntime_steps_seaice+1

     !        Time since begining of model run:
     time = dt_seaice*(nt-1)

     do ibox=1,nsurface_boxes_ocn
        ! If ocean temperature is lower than freezing temperature, or if it 
        ! is higher but sea ice cover exists, then 
        if ( T(2,ibox,1) .lt. T_freezing &
             .or. &
             ( T(2,ibox,1).gt.T_freezing .and. frac_seaice(ibox) .gt. 0.0 ) &
             ) then
           ! calculate heat needed to restore ocean to freezing temp, which 
           ! will also be used for freezing/melting:
           ocean_to_seaice_flux(ibox,1) = &
                -(1/seaice_restoring_time)*Cp_ocn*rho_0*D_ocn(ibox) &
                *(T_freezing-T(2,ibox,1))
       else
          ocean_to_seaice_flux(ibox,1) = 0.0
       end if
     end do
     ! calculate contribution of short wave to melting of seaice:
     shortwave_to_seaice_heat_flux(:) = &
          frac_seaice(:)*shotrwave_penetration_to_seaice &
          *solar_incoming_radiation(:)*(1-albedo_clouds(:))*(1-frac_land(:))*frac_seaice(:) &
          *(1-albedo_seaice)

     !  ------------
     !  step seaice:
     !  ------------
     V_seaice(3,:)=V_seaice(1,:) &
          + 2.*dt_seaice*( &

          ! seaice growth due to heating/ cooling:
          (- ocean_to_seaice_flux(:,1) &
          - shortwave_to_seaice_heat_flux(:) &
          )/(rho_ice*ice_fusion_heat) &

          ! seaice growth due to percipitation falling on sea ice:
          ! dbgz: need to make sure that precipitation falling on land goes 
          ! to ocean, not sea ice:
          + air_to_sea_flux(:,2)*frac_seaice(:)*(rho_ice/rho_0) &

          )*A_ocn(1:nsurface_boxes_ocn)


     ! if seaice volume is negative, we need to form seaice to 
     ! bring its volume to zero.  This means cooling the ocean, so that 
     ! ocean_to_seaice_flux is positive:
     do ibox=1,nsurface_boxes_ocn
        if (V_seaice(3,ibox).lt.0.0) then
           ocean_to_seaice_flux(ibox,1)=ocean_to_seaice_flux(ibox,1) &
                - V_seaice(3,ibox)*(rho_ice*ice_fusion_heat)/(2.0*dt_seaice*A_ocn(ibox))
           V_seaice(3,ibox) = 0.0
        end if
     end do

     ! fresh water flux from ocean to seaice due to melting/ freezing (m/sec):
     ocean_to_seaice_flux(1:nsurface_boxes_ocn,2) = &
          (V_seaice(3,1:nsurface_boxes_ocn)-V_seaice(1,1:nsurface_boxes_ocn)) &
          *(rho_0/rho_ice)/(2.0*dt_seaice*A_ocn(1:nsurface_boxes_ocn))

     !    ---------------------------
     !    Prepare for next time step:
     !    ---------------------------
     ! Robert filter:
     V_seaice(2,1:nsurface_boxes_ocn) = V_seaice(2,1:nsurface_boxes_ocn) &
          + robert_filter_coeff * &
          (   V_seaice(3,1:nsurface_boxes_ocn) &
          - 2*V_seaice(2,1:nsurface_boxes_ocn) &
          +   V_seaice(1,1:nsurface_boxes_ocn) &
          )

     !  normal time step:
     V_seaice(1,:)=V_seaice(2,:)
     V_seaice(2,:)=V_seaice(3,:)

     ! calculate fraction of ocean box area covered by sea ice
     do ibox=1,nsurface_boxes_ocn
        if (V_seaice(3,ibox)/seaice_growth_thickness.lt.A_ocn(ibox)) then
           frac_seaice(ibox)=(V_seaice(3,ibox)/seaice_growth_thickness)/A_ocn(ibox)
           thickness_seaice(ibox)=seaice_growth_thickness
        else
           frac_seaice(ibox)=1.0
           thickness_seaice(ibox)=V_seaice(3,ibox)/A_ocn(ibox)
        endif
     end do

     !   ------------------
     !   Print diagnostics:
     !   ------------------
     if ((nt.eq.1).or. &
          (abs(time- &
          int(time/tprint_short)*tprint_short) &
          .lt. dt_seaice/2.0)) then
     !   write(18,'(50g14.6)') time/(3600.*24.*365.) &
     !        ,u_thermohaline/Sv &
     !        ,(T_atm(3,i),i=1,nboxes_atm)
     end if

     if ((nt.eq.1).or. &
          abs(time-int(time/tprint)*tprint) &
          .lt.dt_seaice/2.0) then
        write(16,'(" ",/," ==>> Time (seaice, years)=",g14.6)') time/year
        write(16,'('' nt(seaice)  ='',i12 &
             &,/,''  ; V_seaice_t (ocean_to_seaice_flux(:,2)*year; m/yr) ='',2g14.6 &
             &,/,''  ; V_seaice    ='',2g14.6 &
             &,/,''  ; thickness_seaice    ='',2g14.6 &
             &,/,''  ; frac_seaice ='',2g14.6 &
             &,/,''  ; ocean_to_seaice_flux(:,1) (watts/m^2) ='',2g14.6 &
             &,/,''  ; shortwave_to_seaice_heat_flux(:) (watts/m^2) ='',2g14.6 &
             &,/,''  ; --> mass balance for seaice: '' &
             &,/,''  ; a: ocean_to_seaice_flux(:,1) (m/year)           ='',2g14.6 &
             &,/,''  ; b: shortwave_to_seaice_heat_flux(:) (m/year)    ='',2g14.6 &
             &,/,''  ; c: air_to_sea_flux(:,2)*frac_seaice(:) (m/year) ='',2g14.6 &
             &,/,''  ; a+b+c =? ocean_to_seaice_flux(:,2) (m/year)     ='',2g14.6 &
             &)') &
             nt &
             ,ocean_to_seaice_flux(:,2)*year &
             ,V_seaice(2,:) &
             ,thickness_seaice(:) &
             ,frac_seaice(:) &
             ,ocean_to_seaice_flux(:,1) &
             ,shortwave_to_seaice_heat_flux(:) &
             ,ocean_to_seaice_flux(:,1)*year*(rho_0/rho_ice)/(rho_ice*ice_fusion_heat) &
             ,shortwave_to_seaice_heat_flux(:)*year*(rho_0/rho_ice)/(rho_ice*ice_fusion_heat)&
             ,air_to_sea_flux(:,2)*frac_seaice(:)*year &
             ,ocean_to_seaice_flux(:,2)*year
     end if
  end do
  ntime_steps_seaice=ntime_steps_seaice+1
end subroutine sea_ice


subroutine calc_FW_flux
  ! calculate meridional fresh water flux (f12_water), using one of 
  ! several optional parameterizations:
  implicit none
  include 'THC_4box_with_seaice.h'
  real*8 Temperature_CC, A_CC, B_CC, RH_CC, Ps_CC, epsilon_CC &
       ,qsat_CC, q_CC, dq_dT_CC, dT_dy, dT_dy_CC, dq_dy, q1, q2
  real*8 K_water_atm_dq_dy_bndry1 !XX

  if (FW_parameterization.eq.'dq_dy_bndry') then
     !  calculate dq/dy at the latitude y_bndry separating the two boxes,
     !   using an approximate Clausius-Clapeyron (CC!) equation.  Note that
     !  dq/dy(y_bndry)={dq/dT(T_bndry)}{dT/dy(y_bndry)}
     Temperature_CC=T_atm(2,1)+(T_atm(2,2)-T_atm(2,1))*L_atm(1)/(L_atm(1)+L_atm(2))
     A_CC=2.53d+11                                ! Pascall
     B_CC=5.42d+3                                 ! K
     RH_CC=0.7                                    ! relative humidity
     Ps_CC=1.0d+5                                 ! surface pressure, pascal
     epsilon_CC=0.622
     qsat_CC=A_CC*exp(-B_CC/Temperature_CC)
     q_CC=RH_CC*epsilon_CC*qsat_CC/Ps_CC
     dq_dT_CC=q_CC*(B_CC/Temperature_CC**2)
     dT_dy_CC=(T_atm(2,2)-T_atm(2,1))/(0.5*(L_atm(1)+L_atm(2)))
     dq_dy=dq_dT_CC*dT_dy_CC

     !XX set to .true. for hysteresis run:
     if (.false.) then
        if (time<tend/2) then
           K_water_atm_dq_dy_bndry1=20d13+180d13*time/(tend/2-tbegin)
        else
           K_water_atm_dq_dy_bndry1=200d13-180d13*(time-tend/2)/(tend/2-tbegin)
        end if
     else
        K_water_atm_dq_dy_bndry1=K_water_atm_dq_dy_bndry
     end if
     f12_water = -K_water_atm_dq_dy_bndry1*dq_dy
  else if (FW_parameterization.eq.'q_target') then
     A_CC=2.53d+11                                ! Pascall
     B_CC=5.42d+3                                 ! K
     RH_CC=0.7                                    ! relative humidity
     Ps_CC=1.0d+5                                 ! surface pressure, pascal
     epsilon_CC=0.622
     qsat_CC=A_CC*exp(-B_CC/T_atm(2,2))  ! NORTHERN BOX atmc humidity
     q_CC=RH_CC*epsilon_CC*qsat_CC/Ps_CC
     dT_dy_CC=(T_atm(2,2)-T_atm(2,1))/(0.5*(L_atm(1)+L_atm(2)))
     f12_water = -K_water_atm_q_target*dT_dy_CC*qsat_CC
  else if (FW_parameterization.eq.'marotzke') then
     Temperature_CC=T_atm(2,1)+(T_atm(2,2)-T_atm(2,1))*L_atm(1)/(L_atm(1)+L_atm(2))
     A_CC=2.53d+11                                ! Pascall
     B_CC=5.42d+3                                 ! K
     RH_CC=0.7                                    ! relative humidity
     Ps_CC=1.0d+5                                 ! surface pressure, pascal
     epsilon_CC=0.622
     q1=RH_CC*A_CC*epsilon_CC*exp(-B_CC/T_atm(2,1))/Ps_CC
     q2=RH_CC*A_CC*epsilon_CC*exp(-B_CC/T_atm(2,2))/Ps_CC
     dT_dy=(T_atm(2,2)-T_atm(2,1))/(0.5*(L_atm(1)+L_atm(2)))
     dq_dy=(q2-q1)/(0.5*(L_atm(1)+L_atm(2)))
     f12_water = -K_water_atm_marotzke*abs(dT_dy)**2*dq_dy
  else
     write(*,*) ' calc_FW_flux FW_parameterization=', FW_parameterization
     stop' calc_FW_flux: unknown value for variable FW_parameterization'
  end if

end subroutine calc_FW_flux




!     ;;; Local Variables: ***
!     ;;; mode:f90 ***
!     ;;; compile-command:"make" ***
!     ;;; comment-start: ";;; "  ***
!     ;;; comment-end:"***" ***
!     ;;; End: ***
