module convert_0

!#########################################################################

!THIS MODULE IMPLEMENTS CONVERSION BETWEEN MASS FRACTIONS AND MOLE FRACTIONS
!AND PROVIDES A CONVERSION FROM PRACTICAL SALINITY TO ABSOLUTE SALINITY

!#########################################################################

!IMPLEMENTATION IN FORTRAN BY D.G. WRIGHT
!FOR PUBLICATION IN OCEAN SCIENCE, AS DESCRIBED IN THE PAPERS

!FEISTEL, R., WRIGHT, D.G., JACKETT, D.R., MIYAGAWA, K., REISSMANN, J.H.,
!WAGNER, W., OVERHOFF, U., GUDER, C., FEISTEL, A., MARION, G.M.:
!NUMERICAL IMPLEMENTATION AND OCEANOGRAPHIC APPLICATION OF THE THERMODYNAMIC
!POTENTIALS OF WATER, VAPOUR, ICE, SEAWATER AND AIR. PART I: BACKGROUND AND EQUATIONS. 
!OCEAN SCIENCES, 2009, IN PREPARATION.

!WRIGHT, D.G., FEISTEL, R., JACKETT, D.R., MIYAGAWA, K., REISSMANN, J.H., 
!WAGNER, W., OVERHOFF, U., GUDER, C., FEISTEL, A., MARION, G.M.:
!NUMERICAL IMPLEMENTATION AND OCEANOGRAPHIC APPLICATION OF THE THERMODYNAMIC
!POTENTIALS OF WATER, VAPOUR, ICE, SEAWATER AND AIR. PART II: THE LIBRARY ROUTINES, 
!OCEAN SCIENCES., 2009, IN PREPARATION.

!FEISTEL, R., KRETZSCHMAR, H.-J., SPAN, R., HAGEN, E., !WRIGHT, D.G., JACKETT, D.R.:
!THERMODYNAMIC PROPERTIES OF SEA AIR.
!OCEAN SCIENCE DISCUSSION 6(2009)2193-2325.

!#########################################################################

!THIS MODULE REQUIRES THE LIBRARY MODULES
!     CONSTANTS_0, FILE CONSTANTS_0.F90

!#########################################################################

use constants_0

implicit none
private

character*16, private :: version = '9 Dec 2009'

public :: air_massfraction_air_si, air_massfraction_vap_si, &
          air_molar_mass_si, air_molfraction_air_si, &
          air_molfraction_vap_si
public :: asal_from_psal, psal_from_asal

real*8 ::  mw = molar_mass_h2o_si      !MOLAR MASS OF H2O IN KG/MOL
real*8 ::  ma = molar_mass_air_L2000   !MOLAR MASS OF AIR IN KG/MOL USED BY LEMMON et al. 2000
integer flag_del_sa

character*16, private :: version = '5 Aug 2009'

contains

!==========================================================================
function air_molar_mass_si(a_si)
!==========================================================================

!OUTPUT:  
!air_molar_mass_si = MOLAR MASS OF HUMID AIR IN KG MOL-1

!INPUT:  
!A_SI = AIR MASS FRACTION, IN KG/KG

!CHECK VALUE: air_molar_mass_si(0.5) = 2.22122197774E-02

real*8 :: air_molar_mass_si, a_si

air_molar_mass_si = errorreturn

if (a_si < 0d0 .or. a_si > 1d0) return

air_molar_mass_si = 1d0 / ((1d0 - a_si) / mw + a_si / ma)

end function

!==========================================================================
function air_molfraction_vap_si(a_si)
!==========================================================================

!OUTPUT: 
!MOLE FRACTION OF VAPOUR IN HUMID AIR IN MOL MOL-1

!INPUT: 
!A_SI = AIR MASS FRACTION, IN KG/KG

!CHECK VALUE: air_molfraction_vap_si(0.5) = 0.616483190186

real*8 air_molfraction_vap_si, a_si

air_molfraction_vap_si = errorreturn

if (a_si < 0d0 .or. a_si > 1d0) return

air_molfraction_vap_si = (1d0 - a_si) / (1d0 - a_si * (1d0 - mw / ma))

end function

!==========================================================================
function air_molfraction_air_si(a_si)
!==========================================================================

!OUTPUT: 
!MOLE FRACTION OF DRY AIR IN HUMID AIR IN MOL MOL-1

!INPUT: 
!A_SI = AIR MASS FRACTION, IN KG/KG

!CHECK VALUE: air_molfraction_air_si(0.5) = 0.383516809814

real*8 air_molfraction_air_si, a_si

air_molfraction_air_si = errorreturn


!==========================================================================
! THIS IS THE GSW ABSOLUTE SALINITY ALGORITHM IMPLEMENTED IN FORTRAN
!==========================================================================
!
! AS DEFINED IN 
!
! MCDOUGALL, T.J., JACKETT, D.R. AND MILLERO, F.J., 2009: AN ALGORITHM FOR ESTIMATING  
!            ABSOLUTE SALINITY IN THE GLOBAL OCEAN, OCEAN SCIENCE, SUBMITTED.  
!     
! 	
! DAVID JACKETT
! JANUARY 2009


!=====================================================
function asal_from_psal(sp,lon0,lat0,p_si)
!=====================================================

! CONVERT PRACTICAL SALINITY TO ABSOLUTE SALINITY 
!
! SP                  : PRACTICAL SALINITY                 [PSU]
! P_SI                : ABSOLUTE PRESSURE                  [PA]
! LONGS0              : LONGITUDE                          [DEG E]     
! LATS0               : LATITUDE                           [DEG N]
!
! RESULT              : ABSOLUTE SALINITY                  [KG/KG]

!CHECK VALUE:
!ASAL_FROM_PSAL(35.52764437773386, 201, -21, 101325+1023d4) = 0.0357
!ASAL_FROM_PSAL(35, 180, 40, 101325 + 2d7) = 3.51888478029404E-02  !NORTH pACIFIC SILICATE
!ASAL_FROM_PSAL(8, 20, 57, 101325) = 8.13338057142857E-03          !CENTRAL BALTIC SEA


!IF THE FILE "GSW_DATA.DAT" IS MISSING, SA = SR IS COMPUTED IN THE OCEAN:
!ASAL_FROM_PSAL(35.5276443777339, 201, -21, 101325 + 1023d4) = 3.56951724471082E-02
!ASAL_FROM_PSAL(35, 180, 40, 101325 + 2d7) = 0.03516504           !north Pacific (without silicate)
!ASAL_FROM_PSAL(8, 20, 57, 101325) = 8.13338057142857E-03         !central Baltic Sea


implicit none
integer flag
real*8 sa, sp, lon0, lat0, p_si, p0, asal_from_psal

p0 = (p_si - 101325d0)/1d4
sa = (35.16504d0/35d0)*sp  + gsw_delta_sa(lon0,lat0,p0)
flag = 1
asal_from_psal = 1d-3 * gsw_adjust_baltic(sa,sp,lon0,lat0,flag)

return
end function


!=====================================================
function psal_from_asal(sa_si,lon0,lat0,p_si)
!=====================================================

! CONVERT ABSOLUTE SALINITY TO PRACTICAL SALINITY 
!
! SA_SI               : ABSOLUTE SALINITY                  [KG/KG]
! P_SI                : ABSOLUTE PRESSURE                  [PA]
! LONGS0              : LONGITUDE                          [DEG E]     
! LATS0               : LATITUDE                           [DEG N]
!
! PSAL_FROM_ASAL  : PRACTICAL SALINITY                     [PSU]

!CHECK VALUE
!PSAL_FROM_ASAL(0.0357, 201, -21, 101325+1023D4) = 35.5276443777

!IF THE FILE "GSW_DATA.DAT" IS MISSING, SA = SR IS COMPUTED IN THE OCEAN:
!PSAL_FROM_ASAL(0.0357, 201, -21, 101325 + 1023d4) = 35.532449273483


implicit none
integer flag
real*8 sa_si, sa, lon0, lat0, p_si, p0, psal_from_asal

p0 = (p_si - 101325d0)/1d4
sa = 1d3 * sa_si
psal_from_asal = (35.d0/35.16504d0)*(sa - gsw_delta_sa(lon0,lat0,p0))
flag = 2
psal_from_asal = gsw_adjust_baltic(sa,psal_from_asal,lon0,lat0,flag);

return
end function


!=====================================================
function gsw_delta_sa(lon0_in,lat0,p0)
!=====================================================

! CALCULATE THE ABSOLUTE SALINITY ANOMALY
!
! P0                  : ABSOLUTE pressure                  [DBAR]
! LONGS0              : LONGITUDE                          [DEG E]     
! LATS0               : LATITUDE                           [DEG N]
!
! RESULT              : ABSOLUTE SALINITY ANOMALY          [G/KG]

implicit none

integer nx, ny, nz
parameter(nx=91,ny=44,nz=45)

integer indx0, indy0, indz, i, j, icalled
integer k, deli(4), delj(4), nmean
real*8 p0, p0_original, lon0, lat0, lon0_in, sa_upper, sa_lower, gsw_delta_sa
real*8 longs(nx), lats(ny), p(nz), del_sa(nz,ny,nx), dlon, dlat
real*8 r1, s1, t1, dsa_mean, dsa(4), ndepth(ny,nx), ndepth_max

data deli/0,1,1,0/, delj/0,0,1,1/

data icalled/0/

save icalled, longs, lats, p, ndepth, del_sa

gsw_delta_sa = 0d0

if(lat0 < -82d0 .or. lat0 > 90d0) return
lon0 = mod(lon0_in, 360d0)
if(lon0 == 0d0) lon0 = 1d-15
if(lon0 == 360d0) lon0 = 360d0 - 1d-15

if(icalled.eq.0) then
   icalled = 1
   open(10,file='gsw_data.dat',status='old',err=1)
   flag_del_sa = 1
   read(10,*) (longs(i), i=1,nx)
   read(10,*) (lats(i), i=1,ny)
   read(10,*) (p(i), i=1,nz)
   read(10,*) ((ndepth(j,i), j=1,ny), i=1,nx)
   read(10,*) (((del_sa(k,j,i), k=1,nz), j=1,ny), i=1,nx)
   close(10)
   go to 2
1  del_sa = 0d0
   flag_del_sa = 0
2  continue
end if

!SET GSW_DELTA_SA = 0 AND RETURN IF THERE IS NO DATA FILE PRESENT
if(flag_del_sa == 0) then; gsw_delta_sa = 0d0; return; endif

dlon = longs(2)-longs(1); dlat = lats(2)-lats(1)

indx0 = floor(1 + (nx-1)*(lon0-longs(1))/(longs(nx)-longs(1)))
if(indx0.eq.nx) then
   indx0 = nx-1
end if

indy0 = floor(1 + (ny-1)*(lat0-lats(1))/(lats(ny)-lats(1)))
if(indy0.eq.ny) then
   indy0 = ny-1
end if

ndepth_max = -1.d0
do k = 1,4
   if(ndepth(indy0+delj(k),indx0+deli(k)).gt.0.d0) &
      ndepth_max = max(ndepth_max,ndepth(indy0+delj(k),indx0+deli(k)))
end do

if(ndepth_max.eq.-1.d0) then
   gsw_delta_sa = 0.d0 !/0.
   return
end if    

p0_original = p0
if(p0.gt.p(int(ndepth_max))) p0 = p(int(ndepth_max))
call indx(p,nz,p0,indz)
    
r1 = (lon0-longs(indx0))/(longs(indx0+1)-longs(indx0));
s1 = (lat0-lats(indy0))/(lats(indy0+1)-lats(indy0));
t1 = (p0-p(indz))/(p(indz+1)-p(indz));

do k = 1,4
   dsa(k) = del_sa(indz,indy0+delj(k),indx0+deli(k))
end do

if(260.d0.le.lon0.and.lon0.le.291.999d0.and.3.4d0.le.lat0.and.lat0.le.19.55d0) then
   dsa = dsa_add_barrier(dsa,lon0,lat0,longs(indx0),lats(indy0),dlon,dlat)
else if(abs(sum(dsa)) >= 1d10) then 
   dsa = dsa_add_mean(dsa,lon0,lat0)
end if

sa_upper = (1.d0-s1)*(dsa(1) + r1*(dsa(2)-dsa(1))) + s1*(dsa(4) + r1*(dsa(3)-dsa(4)))

do k = 1,4
   dsa(k) = del_sa(indz+1,indy0+delj(k),indx0+deli(k))
end do

if(260.d0.le.lon0.and.lon0.le.291.999d0.and.3.4d0.le.lat0.and.lat0.le.19.55d0) then
   dsa = dsa_add_barrier(dsa,lon0,lat0,longs(indx0),lats(indy0),dlon,dlat)
else if(abs(sum(dsa)) >= 1d10) then 
   dsa = dsa_add_mean(dsa,lon0,lat0)
end if

sa_lower = (1.d0-s1)*(dsa(1) + r1*(dsa(2)-dsa(1))) + s1*(dsa(4) + r1*(dsa(3)-dsa(4)))
if(abs(sa_lower) >= 1d10) sa_lower = sa_upper
gsw_delta_sa = sa_upper + t1*(sa_lower-sa_upper)
!dbg
if(abs(gsw_delta_sa) >= 1d10) then
   !write(*,*)"gsw_delta_sa = errorreturn has been replaced by 0"
   gsw_delta_sa = 0.d0
endif
!dbg
p0 = p0_original
  
return
end function


!=====================================================
function gsw_adjust_baltic(sa,sp,longs,lats,flag)
!=====================================================

! FOR THE BALTIC SEA, OVERWRITE ABSOLUTE SALINITY WITH A VALUE
! COMPUTED ANALYTICALLY FROM PRACTICAL SALINITY, OR VICE VERSA
!
! SA                  : ABSOLUTE SALINITY                  [G/KG]
! SP                  : PRACTICAL SALINITY                 [PSU]
! LONGS               : LONGITUDE                          [DEG E]     
! LATS                : LATITUDE                           [DEG N]
! FLAG                : FLAG - 1 or 2 
!
! GSW_ADJUST_bALTIC   : ABSOLUTE SALINITY                  [G/KG]
!                         WHEN FLAG = 1
!                     : PRACTICAL SALINITY                 [PSU]
!                         WHEN FLAG = 2

implicit none
integer n2, n3, flag
real*8 SA, SP, longs, lats, gsw_adjust_Baltic
real*8 xb_left(3), xb_right(2), yb_left(3), yb_right(2), xx_left, xx_right

data xb_left/12.6d0, 7.d0, 26.d0/, yb_left/50.d0, 59.d0, 69.d0/
data xb_right/45.d0, 26.d0/, yb_right/50.d0, 69.d0/

!EXTERNAL XINTERP1

n2 = 2; n3 = 3

if(flag.eq.1) then
   gsw_adjust_Baltic = SA
elseif(flag.eq.2) then
   gsw_adjust_Baltic = SP
end if

if(xb_left(2).lt.longs .and. longs.lt.xb_right(1) .and. yb_left(1).lt.lats .and. lats.lt.yb_left(3)) then
   xx_left = xinterp1(yb_left, xb_left, n3, lats)
   xx_right = xinterp1(yb_right, xb_right, n2, lats)
   if(xx_left.le.longs .and. longs.le.xx_right) then
      if(flag.eq.1) then
         gsw_adjust_Baltic = (35.16504d0/35.d0)*SP + 0.124d0*(1.d0-SP/35.d0);
      elseif(flag.eq.2) then
         gsw_adjust_Baltic = (35.d0/35.04104d0)*(SA - 0.124d0)
      end if
   end if
end if

return
end function


!=====================================================
function dsa_add_barrier(dsa,lon0,lat0,longs,lats,dlon,dlat)
!=====================================================

! ADD A BARRIER THROUGH CENTRAL AMERICA (PANAMA) AND THEN AVERAGE
! OVER THE APPROPRIATE SIDE OF THE BARRIER
! 
! DSA                 : ABSOLUTE SALINITY ANOMALY                [G/KG]
! LONGS0              : LONGITUDES OF DATA                       [DEG E]
! LATS0               : LATITUDES OF DATA                        [DEG N]
! LONGS               : LONGITUDES OF REGULAR GRID               [DEG E]
! LATS                : LATITUDES OF REGULAR GRID                [DEG N]
! DLONGS              : LONGITUDE DIFFERENCE OF REGULAR GRID     [DEG LONGITUDE]
! DLATS               : LATITUDES DIFFERENCE OF REGULAR GRID     [DEG LATITUDE]
!
! RESULT              : ABSOLUTE SALINITY ANOMALY OF DATA        [G/KG]

implicit none

integer n4, n6

parameter(n4=4, n6=6)

integer k, nmean, above_line0, above_line(n4)

real*8 dsa_add_barrier(n4), dsa(n4)
real*8 longs_pan(n6), lats_pan(n6), lon0, lat0, r, lats_line
real*8 longs, lats, dlon, dlat, dsa_mean

data longs_pan/260.0000, 272.5900, 276.5000, 278.6500, 280.7300, 292.000/ 
data  lats_pan/ 19.5500,  13.9700,   9.6000,   8.1000,   9.3300,   3.400/ 


call indx(longs_pan,n6,lon0,k)                            !   the long0/lat0 point
r = (lon0-longs_pan(k))/(longs_pan(k+1)-longs_pan(k))
lats_line = lats_pan(k) + r*(lats_pan(k+1)-lats_pan(k))
if(lats_line.le.lat0) then
   above_line0 = 1
else
   above_line0 = 0
end if

!print *, 'above_line = ', above_line0

call indx(longs_pan,n6,longs,k)                                     !  the 1 and 4 long/lat points 
r = (longs-longs_pan(k))/(longs_pan(k+1)-longs_pan(k))
lats_line = lats_pan(k) + r*(lats_pan(k+1)-lats_pan(k))
if(lats_line.le.lats) then
   above_line(1) = 1
else
   above_line(1) = 0
end if
if(lats_line.le.lats+dlat) then
   above_line(4) = 1
else
   above_line(4) = 0
end if

call indx(longs_pan,n6,longs+dlon,k)                              !  the 2 and 3 long/lat points 
r = (longs+dlon-longs_pan(k))/(longs_pan(k+1)-longs_pan(k))
lats_line = lats_pan(k) + r*(lats_pan(k+1)-lats_pan(k))
if(lats_line.le.lats) then
   above_line(2) = 1
else
   above_line(2) = 0
end if
if(lats_line.le.lats+dlat) then
   above_line(3) = 1
else
   above_line(3) = 0
end if

nmean = 0; dsa_mean = 0.d0; dsa_add_barrier = dsa
do k = 1,n4
   if ((abs(dsa(k)) <= 100d0).and.above_line0.eq.above_line(k)) then
      nmean = nmean+1; dsa_mean = dsa_mean+dsa(k)
   end if
end do

!dbg
if(nmean == 0)then
   dsa_mean = 0d0    !errorreturn
else
   dsa_mean = dsa_mean/nmean
endif
!dbg

do k = 1,n4
   if((abs(dsa(k)) >= 1d10).or.above_line0.ne.above_line(k)) then
      dsa_add_barrier(k) = dsa_mean
   end if
end do

return
end function


!=====================================================
function dsa_add_mean(dsa,lon0,lat0)
!=====================================================

! REPLACE NANS WITH NAMEAN   
!
! DSA           : ABSOLUTE SALINITY ANOMALY                [G/KG]
!                 OF THE 4 ADJACENT NEIGHBOURS     
!
! RESULT        : NANMEAN OF THE 4 ADJACENT NEIGHBOURS     [G/KG]

implicit none

integer k, nmean
real*8 dsa(4)
real*8 dsa_add_mean(4)
real*8 dsa_mean, lon0, lat0

nmean = 0; dsa_mean = 0.d0; dsa_add_mean = dsa

do k = 1,4
   if (abs(dsa(k)) <= 100d0) then
      nmean = nmean+1; dsa_mean = dsa_mean+dsa(k)
   end if
end do

!dbg
if(nmean == 0)then
   dsa_mean = 0d0    !errorreturn
else
   dsa_mean = dsa_mean/nmean
endif
!dbg

do k = 1,4
   if(abs(dsa(k)) >= 100d0) then
      dsa_add_mean(k) = dsa_mean
   end if
end do

return
end function


!=====================================================
subroutine indx(x,n,z,k)
!=====================================================

!DESCRIPTION:	FIND THE INDEX OF A REAL NUMBER IN A
!		MONOTONICALLY INCREASING REAL ARRAY
!
!INPUT:	          X		ARRAY OF INCREASING VALUES
!                 N		LENGTH OF ARRAY
!	          Z		REAL NUMBER
!
!OUTPUT:	  K		INDEX K - IF X(K) <= Z < X(K+1), OR
!				N-1     		- IF Z = X(N)
!
!CREATED:		JUNE 1993
!

real*8 x, z
integer n, k, ku, kl, km

dimension x(n)

if(x(1).lt.z.and.z.lt.x(n)) then
   kl=1
   ku=n
   do while (ku-kl.gt.1)
      km=(ku+kl)/2
      if(z.gt.x(km))then
         kl=km
      else
         ku=km
      endif
   end do
   k=kl
   if(z.eq.x(k+1)) k = k+1
else
   if(z.eq.x(1)) then
      k = 1
   elseif(z.eq.x(n)) then
      k = n-1
   else
      write(*,*) 'ERROR in indx.f : out of range'
      write(*,*) 'z = ', z, 'n = ', n, 'x = ',x
   end if
end if

return
end subroutine


!=====================================================
function xinterp1(x,y,n,x0)
!=====================================================

! LINEARLY INTERPOLATE A REAL ARRAY   
!
! X                   : X ARRAY (MONOTONIC)               
! Y                   : Y ARRAY     
! N                   : LENGTH OF X AND Y ARRAYS
! X0                  : X VALUE TO BE INTERPOLATED
!
! RESULT              : INTERPOLATED VALUE

implicit none
integer n, k
real*8 x(n), y(n), x0, r, xinterp1

call indx(x,n,x0,k)
r = (x0-x(k))/(x(k+1)-x(k))
xinterp1 = y(k) + r*(y(k+1)-y(k))

return
end function

end module convert_0