coupling_sltn.F90

Go to the documentation of this file.

!SFX_LIC Copyright 1994-2014 CNRS, Meteo-France and Universite Paul Sabatier
!SFX_LIC This is part of the SURFEX software governed by the CeCILL-C licence
!SFX_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt  
!SFX_LIC for details. version 1.
SUBROUTINE coupling_slt_n (SLT, &
      ki,                   &!I [nbr] number of sea points 
      kslt,                 &!I [nbr] number of sea points 
      pwind,                &!I Wind velocity
      psfslt                &!O [kg/m2/sec] flux of sea salt
      ) 
  
!PURPOSE
!-------
!  Compute sea salt emission  upon Vignatti et al, 2001
!
!AUTHOR
!-------
! P. Tulet
!
!
USE modd_slt_n, ONLY : slt_t
!
USE modd_csts, ONLY : xavogadro, xpi
USE modd_slt_surf
!
USE yomhook   ,ONLY : lhook,   dr_hook
USE parkind1  ,ONLY : jprb
!
IMPLICIT NONE
!
!INPUT
!
TYPE(slt_t), INTENT(INOUT) :: slt
!
INTEGER, INTENT(IN)                :: ki             !I Number of sea points
INTEGER, INTENT(IN)                :: kslt           !I Number of sea salt emission variables
REAL, DIMENSION(KI),      INTENT(IN)  :: pwind       !I wind velocity
REAL, DIMENSION(KI,KSLT), INTENT(OUT) :: psfslt      !Out: mole particles per mole air m/s *(MWdst/MWair*rhoair)(index #1)
                                                     !Out: kg/m2/s (index #2)
                                                     !Out: moles m6/moles air m/s *(MWdst/MWair*rhoair)(index #3)
!LOCAL VARIABLES
REAL,DIMENSION(KI,3)  :: zsfslt_mde       ! sea salt flux from modes
INTEGER               :: jn, ji           !Counter for sea salt modes
REAL, DIMENSION(KI)   :: dzspeed 
INTEGER, DIMENSION(KI):: wcl
REAL                  :: zconvertfacm0_slt
REAL                  :: zconvertfacm3_slt
REAL                  :: zconvertfacm6_slt
!
!REAL, PARAMETER :: MASS1FLUX(0:40) = (/  &
!          0.000E+00, 2.483E-15, 2.591E-14, 1.022E-13, 2.707E-13, 5.761E-13,  &
!          1.068E-12, 1.800E-12, 2.829E-12, 4.215E-12, 6.023E-12, 8.317E-12, &
!          1.117E-11, 1.464E-11, 1.882E-11, 2.378E-11, 2.959E-11, 3.633E-11, &
!          4.409E-11, 5.296E-11, 6.301E-11, 7.433E-11, 8.693E-11, 1.012E-10, &
!          1.168E-10, 1.342E-10, 1.532E-10, 1.741E-10, 1.970E-10, 2.219E-10, &
!          2.489E-10, 2.781E-10, 3.097E-10, 3.437E-10, 3.803E-10, 4.195E-10, &
!          4.616E-10, 5.065E-10, 5.544E-10, 6.054E-10, 6.711E-10             /) 
!
!REAL, PARAMETER :: MASS2FLUX(0:40) = (/  &
!          0.000E+00, 2.319E-13, 2.411E-12, 9.481E-12, 2.505E-11, 5.321E-11,  &
!          9.850E-11, 1.658E-10, 2.602E-10, 3.874E-10, 5.529E-10, 7.628E-10,  &
!          1.023E-09, 1.341E-09, 1.722E-09, 2.175E-09, 2.704E-09, 3.319E-09,  &
!          4.026E-09, 4.832E-09, 5.746E-09, 6.776E-09, 7.925E-09, 9.214E-09,  &
!          1.064E-08, 1.221E-08, 1.394E-08, 1.584E-08, 1.791E-08, 2.016E-08,  &
!          2.261E-08, 2.526E-08, 2.812E-08, 3.120E-08, 3.451E-08, 3.806E-08,  &
!          4.186E-08, 4.592E-08, 5.025E-08, 5.486E-08, 6.014E-08             /) 
!
!REAL, PARAMETER :: MASS3FLUX(0:40) = (/ 0.0, &
!        1.783E-12, 1.579E-11, 5.852E-11, 1.501E-10, 3.134E-10, 5.740E-10, &
!        9.597E-10, 1.500E-09, 2.227E-09, 3.175E-09, 4.378E-09, 5.872E-09, &
!        7.698E-09, 9.897E-09, 1.250E-08, 1.556E-08, 1.912E-08, 2.323E-08, &
!        2.792E-08, 3.325E-08, 3.927E-08, 4.608E-08, 5.356E-08, 6.194E-08, &
!        7.121E-08, 8.143E-08, 9.266E-08, 1.049E-07, 1.183E-07, 1.329E-07, &
!        1.487E-07, 1.658E-07, 1.843E-07, 2.041E-07, 2.255E-07, 2.484E-07, &
!        2.729E-07, 2.991E-07, 3.270E-07, 3.517E-07 /) 


REAL, PARAMETER :: numb1flux(0:40) = (/ &
         0.000e+00, 3.004e+01, 3.245e+02, 1.306e+03, 3.505e+03, 7.542e+03,  &
          1.410e+04, 2.394e+04, 3.787e+04, 5.674e+04, 8.147e+04, 1.130e+05,  &
          1.523e+05, 2.005e+05, 2.586e+05, 3.278e+05, 4.091e+05, 5.037e+05,  &
          6.129e+05, 7.379e+05, 8.800e+05, 1.041e+06, 1.220e+06, 1.422e+06,  &
          1.646e+06, 1.893e+06, 2.166e+06, 2.466e+06, 2.794e+06, 3.152e+06,  &
          3.541e+06, 3.962e+06, 4.419e+06, 4.911e+06, 5.441e+06, 6.011e+06,  &
          6.621e+06, 7.274e+06, 7.972e+06, 8.716e+06, 8.801e+06             /) 

REAL, PARAMETER :: numb2flux(0:40) = (/  &
          0.000e+00, 1.934e+01, 2.068e+02, 8.271e+02, 2.211e+03, 4.741e+03,  &
          8.841e+03, 1.497e+04, 2.363e+04, 3.534e+04, 5.066e+04, 7.017e+04,  &
          9.447e+04, 1.242e+05, 1.600e+05, 2.025e+05, 2.525e+05, 3.106e+05,  &
          3.776e+05, 4.542e+05, 5.413e+05, 6.395e+05, 7.501e+05, 8.726e+05,  &
          1.009e+06, 1.160e+06, 1.327e+06, 1.509e+06, 1.709e+06, 1.927e+06,  &
          2.163e+06, 2.420e+06, 2.697e+06, 2.996e+06, 3.318e+06, 3.664e+06,  &
          4.034e+06, 4.430e+06, 4.852e+06, 5.303e+06, 5.740e+06             /) 

REAL, PARAMETER :: numb3flux(0:40) = (/ 0.0, &
        4.340e-01, 5.217e+00, 2.241e+01, 6.301e+01, 1.404e+02, 2.703e+02, &
        4.699e+02, 7.584e+02, 1.157e+03, 1.687e+03, 2.373e+03, 3.240e+03, &
        4.314e+03, 5.625e+03, 7.197e+03, 9.063e+03, 1.126e+04, 1.380e+04, &
        1.674e+04, 2.011e+04, 2.393e+04, 2.827e+04, 3.311e+04, 3.853e+04, &
        4.457e+04, 5.126e+04, 5.864e+04, 6.675e+04, 7.564e+04, 8.535e+04, &
        9.592e+04, 1.074e+05, 1.198e+05, 1.333e+05, 1.478e+05, 1.633e+05, &
        1.801e+05, 1.980e+05, 2.172e+05, 2.353e+05 /) 
REAL(KIND=JPRB) :: zhook_handle
!!
!!    MESONH carries the following units during transport:
!!    M0=#/molec_{air}
!!    M6=um6/molec_{air}*1.d6
!!    The surface model should have (for sea salt)
!!    M0=#/m3*[kg_{dst}/mole_{dst}/XAVOGADRO]
!!    M3=kg/m3
!!    M6=um6/m3
!!    REFERENCE
!!    ---------
!!    Tulet et al, ORILAM manuscript for transformation of modal parameters
!!    J. Geophys. Res., 110, D18201, doi:10.1029/2004JD005716
!
!Initialize output which is total flux of sea salt (kg/m2/sec). 
IF (lhook) CALL dr_hook('COUPLING_SLT_N',0,zhook_handle)
!
!Factor which is needed so that all gains normal units when leaving ground paramn
zconvertfacm0_slt = xmolarweight_slt / xavogadro !(kg_dst/mol_dst)/(molec/mol)
!Factor which is needed for moment 6, there is a factor 1.d6 transported around in M6 in MESONH
zconvertfacm6_slt = xmolarweight_slt / xavogadro*1.d6
zconvertfacm3_slt = 4./3.*xpi*xdensity_slt / 1.d18
!
psfslt(:,:)=0.d0
!
IF (cemisparam_slt.eq."Vig01")THEN
  !
  ! Vignatti et al. 2001 (in particles.cm-2.s-1)
  zsfslt_mde(:,1) =  10**(0.09  *pwind(:) + 0.283)   ! fine mode
  zsfslt_mde(:,2) =  10**(0.0422*pwind(:) + 0.288)   ! median mode
  zsfslt_mde(:,3) =  10**(0.069 *pwind(:) - 3.5)     ! coarse mode
  ! convert into  particles.m-2.s-1)
  zsfslt_mde(:,1) = max(zsfslt_mde(:,1) * 1e4, 1e-10)  
  zsfslt_mde(:,2) = max(zsfslt_mde(:,2) * 1e4, 1e-10)  
  zsfslt_mde(:,3) = max(zsfslt_mde(:,3) * 1e4, 1e-10)  
  !
ELSE ! Use Schultz et al., 2004
  !
  wcl(:) = int(pwind(:))
  wcl(:) = max(0, min(wcl(:), 39))
  !
  dzspeed(:) = max(0., min(pwind(:) - float(wcl(:)), 1.))
  !
  ! Flux given  in  particles.m-2 s-1
  !
  DO ji=1,ki
    !plm-gfortran
    zsfslt_mde(ji,1) = numb1flux(wcl(ji)) + (numb1flux(wcl(ji)+1)-numb1flux(wcl(ji)))*dzspeed(ji)
    zsfslt_mde(ji,2) = numb2flux(wcl(ji)) + (numb2flux(wcl(ji)+1)-numb2flux(wcl(ji)))*dzspeed(ji)
    zsfslt_mde(ji,3) = numb3flux(wcl(ji)) + (numb3flux(wcl(ji)+1)-numb3flux(wcl(ji)))*dzspeed(ji)
    !plm-gfortran
  ENDDO
  !
END IF
!
DO jn=1,jpmode_slt
  !
  IF (lvarsig_slt) THEN
    !
    psfslt(:,1+(jn-1)*3) = zsfslt_mde(:,jorder_slt(jn))
    psfslt(:,2+(jn-1)*3) = psfslt(:,1+(jn-1)*3) * (slt%XEMISRADIUS_SLT(jn)**3)*exp(4.5 * log(slt%XEMISSIG_SLT(jn))**2)  
    psfslt(:,3+(jn-1)*3) = psfslt(:,1+(jn-1)*3) * (slt%XEMISRADIUS_SLT(jn)**6)*exp(18. * log(slt%XEMISSIG_SLT(jn))**2)  
    !
    ! Conversion into fluxes
    psfslt(:,1+(jn-1)*3) = psfslt(:,1+(jn-1)*3) * zconvertfacm0_slt
    psfslt(:,2+(jn-1)*3) = psfslt(:,1+(jn-1)*3) * zconvertfacm3_slt
    psfslt(:,3+(jn-1)*3) = psfslt(:,3+(jn-1)*3) * zconvertfacm6_slt

  ELSE IF (lrgfix_slt) THEN
    psfslt(:,jn) =  zsfslt_mde(:,jorder_slt(jn)) * (slt%XEMISRADIUS_SLT(jn)**3)*exp(4.5 * log(slt%XEMISSIG_SLT(jn))**2) 
    ! Conversion into fluxes
    psfslt(:,jn) = psfslt(:,jn) * zconvertfacm3_slt

  ELSE
    psfslt(:,1+(jn-1)*2) = zsfslt_mde(:,jorder_slt(jn)) 
    psfslt(:,2+(jn-1)*2) = psfslt(:,1+(jn-1)*2) * (slt%XEMISRADIUS_SLT(jn)**3)*exp(4.5 * log(slt%XEMISSIG_SLT(jn))**2) 

    ! Conversion into fluxes
    psfslt(:,1+(jn-1)*2) = psfslt(:,1+(jn-1)*2) * zconvertfacm0_slt
    psfslt(:,2+(jn-1)*2) = psfslt(:,1+(jn-1)*2) * zconvertfacm3_slt

  ENDIF
END DO

IF (lhook) CALL dr_hook('COUPLING_SLT_N',1,zhook_handle)
END SUBROUTINE coupling_slt_n