MODULE OceanModel USE PhysicalFunctions , ONLY : fpvs IMPLICIT NONE ! Selecting Kinds INTEGER, PARAMETER :: r4 = SELECTED_REAL_KIND(6) ! Kind for 32-bits Real Numbers INTEGER, PARAMETER :: i4 = SELECTED_INT_KIND(9) ! Kind for 32-bits Integer Numbers INTEGER, PARAMETER :: r8 = SELECTED_REAL_KIND(15) ! Kind for 64-bits Real Numbers INTEGER, PARAMETER :: i8 = SELECTED_INT_KIND(14) ! Kind for 64-bits Integer Numbers INTEGER, PARAMETER :: r16 = SELECTED_REAL_KIND(15)! Kind for 128-bits Real Numbers !*L------------------COMDECK C_EPSLON----------------------------------- ! EPSILON IS RATIO OF MOLECULAR WEIGHTS OF WATER AND DRY AIR REAL(KIND=r8), PARAMETER :: Epsilonk = 0.62198_r8 REAL(KIND=r8), PARAMETER :: C_Virtual = 1.0_r8/Epsilonk-1.0_r8 REAL(kind=r8), PARAMETER :: con_rd =2.8705e+2_r8 ! gas constant air (J/kg/K) REAL(kind=r8), PARAMETER :: con_rv =4.6150e+2_r8 ! gas constant H2O (J/kg/K) REAL(kind=r8), PARAMETER:: EPS =con_rd/con_rv REAL(kind=r8), PARAMETER:: EPSM1 =con_rd/con_rv-1.0_r8 !*---------------------------------------------------------------------- !*L------------------COMDECK C_O_DG_C----------------------------------- ! ZERODEGC IS CONVERSION BETWEEN DEGREES CELSIUS AND KELVIN ! TFS IS TEMPERATURE AT WHICH SEA WATER FREEZES ! TM IS TEMPERATURE AT WHICH FRESH WATER FREEZES AND ICE MELTS REAL(KIND=r8), PARAMETER :: ZeroDegC = 273.15_r8 REAL(KIND=r8), PARAMETER :: TFS = 271.35_r8 REAL(KIND=r8), PARAMETER :: TM = 273.15_r8 !*---------------------------------------------------------------------- ! !----------------------------------------------------------------------- ! Local parameters and other physical constants !----------------------------------------------------------------------- REAL(KIND=r8), PARAMETER :: one_minus_epsilon = 1.0_r8 -epsilonk ! One minus the ratio of the molecular weights of water and dry ! air REAL(KIND=r8) , PARAMETER :: hltm = 2.52e6_r8 ! latent heat of vaporization (J kg^-1) REAL(KIND=r8), PARAMETER :: & & T_low = 183.15_r8 & ! Lowest temperature for which look_up table of ! saturation water vapour presssure is valid (K) &,T_high = 338.15_r8 & ! Highest temperature for which look_up table of ! saturation water vapour presssure is valid (K) &,delta_T = 0.1_r8 ! temperature increment of the look-up table ! of saturation vapour pressures. INTEGER, PARAMETER :: & & N = ((T_high - T_low + (delta_T*0.5_r8))/delta_T) + 1.0_r8 ! Gives N=1551, size of the look-up table of saturation water ! vapour pressure. INTEGER IES ! LOOP COUNTER FOR DATA STATEMENT LOOK-UP TABLE ! Local dynamic arrays------------------------------------------------- REAL(KIND=r8) ES(0:N+1) ! TABLE OF SATURATION WATER VAPOUR PRESSURE (PA) ! - SET BY DATA STATEMENT CALCULATED FROM THE ! GOFF-GRATCH FORMULAE AS TAKEN FROM LANDOLT- ! BORNSTEIN, 1987 NUMERICAL DATA AND FUNCTIONAL ! RELATIONSHIPS IN SCIENCE AND TECHNOLOGY. ! GROUP V/ VOL 4B METEOROLOGY. PHYSICAL AND ! CHEMICAL PROPERTIES OF AIR, P35 ! INTEGER,PARAMETER :: NDIV = 6 ! number of particle size divisions REAL(KIND=r8), PARAMETER :: CHARNOCK=.014_r8 INTEGER , PARAMETER ::ISeaZ0T =1 ! ! Options for marine boundary layers INTEGER, PARAMETER :: Fixed_Z0T = 0 ! Stanard flixed value of thermal roughness length over sea INTEGER, PARAMETER :: SurfDivZ0T = 1 ! Thermal roughness length over sea defined from surface ! divergence theory REAL(KIND=r8), PARAMETER :: SeaSalinityFactor=1.0_r8 ! ! Factor allowing for the ! ! effect of the salinity of ! ! sea water on the evaporative ! ! flux. LOGICAL,PARAMETER :: lq_mix_bl =.FALSE. !TRUE ! IN TRUE if mixing ratios used in boundary layer code LOGICAL,PARAMETER :: L_DUST =.FALSE.! IN switch for mineral dust LOGICAL,PARAMETER :: SFME =.FALSE. ! IN STASH flag for wind mixing energy flux LOGICAL,PARAMETER :: QSAT_opt =.FALSE. REAL(KIND=r8), ALLOCATABLE :: si(:) REAL(KIND=r8), ALLOCATABLE :: sl(:) REAL(KIND=r8), ALLOCATABLE :: del(:) CONTAINS SUBROUTINE InitOceanModel(si_in,sl_in,del_in,kMax) IMPLICIT NONE INTEGER, INTENT(IN ) :: kMax REAL(KIND=r8), INTENT(IN ) :: si_in(kMAx+1) REAL(KIND=r8), INTENT(IN ) :: sl_in(kMAx) REAL(KIND=r8), INTENT(IN ) :: del_in(kMAx) ALLOCATE(si(kMax+1)) ALLOCATE(sl(kMax)) ALLOCATE(del(kMax)) si=si_in sl=sl_in del=del_in CALL qsat_data() END SUBROUTINE InitOceanModel ! *****************************COPYRIGHT******************************* ! (C) Crown copyright Met Office. All rights reserved. ! For further details please refer to the file COPYRIGHT.txt ! which you should have received as part of this distribution. ! *****************************COPYRIGHT******************************* !!! SUBROUTINE SF_EXCH------------------------------------------------ !!! !!! Purpose: Calculate coefficients of turbulent exchange between !!! the surface and the lowest atmospheric layer, and !!! "explicit" fluxes between the surface and this layer. !!! !!! Suitable for Single Column use. !!! !!! !!! Programming standard: Unified Model Documentation Paper No 4, !!! Version 2, dated 18/1/90. !!! !!! System component covered: Part of P243. !!! !!! Project task: !!! !!! Documentation: UM Documentation Paper No 24, section P243. !!! See especially sub-section (ix). !!! !!!--------------------------------------------------------------------- ! Arguments :- SUBROUTINE SF_EXCH (& nCols , & kMAx , & CF , & QCF , & QCL , & Q , & T , & gu , & gv , & sinclt , & PS , & RADNET , & TI , & ZH ,& mskant , & iMask , & TSTAR_LAND , & TSTAR_SSI , & TSTAR_SEA ,& TSTAR_SICE ,& rmi , & rhi , & evap , & sens , & ustar , & Z0MSEA & ) IMPLICIT NONE INTEGER, INTENT(IN ) :: nCols ! IN Number of X points? INTEGER, INTENT(IN ) :: kMAx ! Factor allowing for the effect of the salinity of ! sea water on the evaporative flux. REAL(KIND=r8) , INTENT(IN ) :: CF (nCols,kMAx) ! IN Cloud fraction (decimal). REAL(KIND=r8) , INTENT(IN ) :: QCF (nCols,kMAx) ! IN Cloud ice (kg per kg air) REAL(KIND=r8) , INTENT(IN ) :: QCL (nCols,kMAx) ! IN Cloud liquid water REAL(KIND=r8) , INTENT(IN ) :: Q (nCols,kMAx) ! IN specific humidity REAL(KIND=r8) , INTENT(IN ) :: T (nCols,kMAx) ! IN temperature REAL(KIND=r8), INTENT(in ) :: gu (ncols,kMax) REAL(KIND=r8), INTENT(in ) :: gv (ncols,kMax) REAL(KIND=r8), INTENT(in ) :: sinclt(ncols) REAL(KIND=r8) , INTENT(IN ) :: PS (nCols) REAL(KIND=r8) , INTENT(IN ) :: RADNET(nCols) ! IN Net surface radiation (W/m2) positive downwards REAL(KIND=r8) , INTENT(IN ) :: TI (nCols) ! IN Temperature of sea-ice surface layer (K) REAL(KIND=r8) , INTENT(IN ) :: ZH (nCols) ! IN Depth of boundary layer (m). boundary layer (metres). ! Switch for the treatment of the thermal roughness ! length over the sea. ! INTEGER(KIND=i8) , INTENT(IN ) :: mskant (nCols) ! IN Fraction of gridbox which is sea-ice. INTEGER(KIND=i8) , INTENT(IN ) :: iMask (nCols) ! IN Land fraction on all tiles. ! ! Note, if the grid used is staggered in ! ! the vertical, Z1_UV and Z1_TQ can be ! ! different. REAL(KIND=r8) , INTENT(IN ) :: TSTAR_LAND(nCols) ! IN Land mean surface temperature (K). REAL(KIND=r8) , INTENT(IN ) :: TSTAR_SSI (nCols) ! IN Mean sea surface temperature (K). REAL(KIND=r8) , INTENT(IN ) :: TSTAR_SEA (nCols) ! IN Open sea surface temperature (K). REAL(KIND=r8) , INTENT(IN ) :: TSTAR_SICE(nCols) ! IN Sea-ice surface temperature (K). REAL(KIND=r8) , INTENT(OUT ) :: rmi (nCols) REAL(KIND=r8) , INTENT(OUT ) :: rhi (nCols) REAL(KIND=r8) , INTENT(OUT ) :: evap (nCols) REAL(KIND=r8) , INTENT(OUT ) :: sens (nCols) REAL(KIND=r8) , INTENT(OUT ) :: ustar (nCols) REAL(KIND=r8) , INTENT(INOUT) :: Z0MSEA (nCols) ! INOUT Sea-surface roughness length for momentum (m). F617. REAL(KIND=r8) :: Z1_TQ (nCols) ! IN Height of lowest tq level (m). REAL(KIND=r8) :: Z1_UV (nCols) ! IN Height of lowest uv level (m). REAL(KIND=r8) :: PSTAR (nCols) ! IN Surface pressure (Pascals). REAL(KIND=r8) :: VSHR_SSI (nCols) ! IN Mag. of mean sea sfc-to-lowest-level wind shear REAL(KIND=r8) :: TL_1(nCols) ! IN Liquid/frozen water temperature for lowest atmospheric layer (K). REAL(KIND=r8) :: BQ_1(nCols) ! IN A buoyancy parameter for lowest atm ! ! level ("beta-q twiddle"). REAL(KIND=r8) :: BT_1(nCols) ! IN A buoyancy parameter for lowest atm ! ! level ("beta-T twiddle"). REAL(KIND=r8) :: QW_1(nCols) ! IN Total water content of lowest ! ! atmospheric layer (kg per kg air). REAL(KIND=r8) :: RHOSTAR(nCols) ! OUT Surface air density REAL(KIND=r8) :: FLANDG (nCols) ! IN Land fraction on all tiles. REAL(KIND=r8) :: ICE_FRACT (nCols) ! IN Fraction of gridbox which is sea-ice. ! Define local storage. ! (a) Workspace. REAL(KIND=r8) :: QS1(nCols) ! Sat. specific humidity ! ! qsat(TL_1,PSTAR) REAL(KIND=r8) :: QSTAR_SEA(nCols) ! Surface saturated sp humidity REAL(KIND=r8) :: QSTAR_ICE(nCols)! Surface saturated sp humidity ! REAL(KIND=r8) :: PSTAR_LAND(LAND_PTS)! Surface pressure for land points. REAL(KIND=r8) :: Z0H_SEA(nCols)! Roughness length for heat and ! ! moisture transport INTEGER :: NSICE ! Number of sea-ice points. INTEGER :: SICE_INDEX (nCols,2)! Index of sea-ice points REAL(KIND=r8) :: Z1_TQ_SEA(nCols)! Height of lowest model level ! ! relative to sea. REAL(KIND=r8) :: Z0_ICE(nCols) ! Roughness length. REAL(KIND=r8) :: Z0_MIZ(nCols) ! Roughness length. REAL(KIND=r8) :: RIB_SEA(nCols)! Bulk Richardson number REAL(KIND=r8) :: RIB_ICE(nCols)! Bulk Richardson number REAL(KIND=r8) :: DB_SEA(nCols) ! Buoyancy difference for sea points REAL(KIND=r8) :: DB_ICE(nCols) ! Buoyancy difference for sea ice REAL(KIND=r8) :: BT (nCols,kMAx) REAL(KIND=r8) :: BQ (nCols,kMAx) REAL(KIND=r8) :: P (nCols,kMAx) REAL(KIND=r8) :: BT_CLD (nCols,kMAx) REAL(KIND=r8) :: BQ_CLD (nCols,kMAx) REAL(KIND=r8) :: BT_GB (nCols,kMAx) REAL(KIND=r8) :: BQ_GB (nCols,kMAx) REAL(KIND=r8) :: A_QS (nCols,kMAx) REAL(KIND=r8) :: A_DQSDT(nCols,kMAx) REAL(KIND=r8) :: DQSDT (nCols,kMAx) REAL(KIND=r8) :: QW(nCols, kMAx)! OUT Total water content REAL(KIND=r8) :: TL(nCols, kMAx)! OUT Ice/liquid water temperature REAL(KIND=r8) :: ALPHA1_SICE(nCols) REAL(KIND=r8) :: ASHTF(nCols) REAL(KIND=r8) :: E_SEA(nCols) REAL(KIND=r8) :: FTL_ICE(nCols) REAL(KIND=r8) :: FTL_1(nCols) REAL(KIND=r8) :: H_SEA(nCols) REAL(KIND=r8) :: RHO_ARESIST(nCols) ! ! OUT RHOSTAR*CD_STD*VSHR for SCYCLE ! REAL(KIND=r8) :: RHO_ARESIST_LAND(nCols) ! ! Land mean of rho_aresist_tile REAL(KIND=r8) :: ARESIST(nCols) ! ! OUT 1/(CD_STD*VSHR) for SCYCLE REAL(KIND=r8) :: RESIST_B(nCols) ! ! OUT (1/CH-1/CD_STD)/VSHR for SCYCLE ! REAL(KIND=r8) :: RHOKM_LAND(nCols,nCols) ! ! OUT For land momentum. NB: This is output REAL(KIND=r8) :: RHOKH_1_SICE(nCols) ! ! OUT Surface exchange coefficient for sea ! ! or sea-ice. ! OUT "Explicit" surface flux of TL = H/CP. ! ! (sensible heat / CP). grid-box mean REAL(KIND=r8) :: FQW_1(nCols) ! OUT "Explicit" surface flux of QW (i.e. ! ! evaporation), on P-grid (kg/m2/s). ! ! for whole grid-box REAL(KIND=r8) :: FQW_ICE(nCols) ! ! OUT GBM FQW_1 for sea-ice. REAL(KIND=r8) :: RHOKM_SSI(nCols) ! ! OUT For mean sea mom. NB: This is output ! ! on UV-grid, but with the first and ! ! last set to "missing data". REAL(KIND=r8) :: RHOKPM_SICE(nCols) ! ! OUT Mixing coefficient for sea-ice. REAL(KIND=r8) :: CD_SSI(nCols) ! ! OUT Bulk transfer coefficient for ! ! momentum over sea mean. REAL(KIND=r8) :: CH_SSI(nCols) ! ! OUT Bulk transfer coefficient for heat ! ! and/or moisture over sea mean. REAL(KIND=r8) :: CD_LAND(nCols) ! Bulk transfer coefficient for ! ! momentum over land. REAL(KIND=r8) :: CD(nCols) ! OUT Bulk transfer coefficient for ! ! momentum. REAL(KIND=r8) :: CH(nCols) ! OUT Bulk transfer coefficient for heat ! ! and/or moisture. REAL(KIND=r8) :: CD_STD_DUST(nCols) !OUT Bulk transfer coef. for ! ! momentum, excluding orographic effects ! Workspace for sea and sea-ice leads REAL(KIND=r8) :: CD_ICE(nCols) ! Drag coefficient REAL(KIND=r8) :: CD_SEA(nCols) ! Drag coefficient REAL(KIND=r8) :: CH_SEA(nCols) ! Transfer coefficient for heat and ! ! moisture REAL(KIND=r8) :: CD_MIZ(nCols) ! Drag coefficient REAL(KIND=r8) :: CH_ICE(nCols) ! Transfer coefficient for heat and ! ! moisture REAL(KIND=r8) :: CH_MIZ(nCols) ! Transfer coefficient for heat and ! ! moisture REAL(KIND=r8) :: RECIP_L_MO_SEA(nCols) ! ! OUT Reciprocal of the Monin-Obukhov ! ! length for sea points (m^-1) REAL(KIND=r8) :: RECIP_L_MO_MIZ(nCols) ! ! Reciprocal of the Monin-Obukhov ! ! length for marginal sea ice (m^-1). REAL(KIND=r8) :: RECIP_L_MO_ICE(nCols) ! ! Reciprocal of the Monin-Obukhov ! ! length for sea ice (m^-1). REAL(KIND=r8) :: V_S_SEA(nCols) ! Surface layer scaling velocity ! ! for sea points (m/s). REAL(KIND=r8) :: V_S_ICE(nCols) ! Surface layer scaling velocity ! ! for sea ice (m/s). REAL(KIND=r8) :: V_S_MIZ(nCols) ! Surface layer scaling velocity ! ! for marginal sea ice (m/s). REAL(KIND=r8) :: Q1_SD(nCols) ! OUT Standard deviation of turbulent ! ! fluctuations of surface layer ! ! specific humidity (kg/kg). REAL(KIND=r8) :: T1_SD(nCols) ! OUT Standard deviation of turbulent ! ! fluctuations of surface layer ! ! temperature (K). REAL(KIND=r8) :: speedm(nCols) REAL(KIND=r8) :: FME(nCols) ! OUT Wind mixing energy flux (Watts/sq m). REAL(KIND=r8) :: VSHR(nCols) REAL(KIND=r8) :: TAU ! Magnitude of surface wind stress over sea. REAL(KIND=r8) :: TOL_USTR_N ! Tolerance for USTR_N REAL(KIND=r8) :: TOL_USTR_L ! Tolerance for USTR_L (see below) REAL(KIND=r8) :: USTR_L ! Low-wind estimate of friction velocity REAL(KIND=r8) :: USTR_N ! Neutral estimate of friction velocity REAL(KIND=r8) :: R_B_DUST(nCols,NDIV) !OUT surf layer res for dust REAL(KIND=r8) :: qss(nCols) INTEGER :: I,K INTEGER :: JITS,IDIV !!---------------------------------------------------------------------- !!!-----------COMDECK C_ROUGH FOR SUBROUTINE SF_EXCH---------- ! Sea ice parameters ! Z0FSEA = roughness length for free convective heat and moisture ! transport over the sea (m). ! DUMMY VARIABLE - Only used in 7A boundary layer scheme ! Z0HSEA = roughness length for heat and moisture transport ! over the sea (m). ! Z0MIZ = roughness length for heat, moisture and momentum over ! the Marginal Ice Zone (m). ! Z0SICE = roughness length for heat, moisture and momentum over ! sea-ice (m). REAL(KIND=r8), PARAMETER :: Z0HSEA=1.14e-5_r8 REAL(KIND=r8), PARAMETER :: Z0MIZ=1.5e-5_r8 REAL(KIND=r8), PARAMETER :: Z0SICE=0.0002e0_r8! !!*---------------------------------------------------------------------- !*L------------------COMDECK C_G---------------------------------------- ! G IS MEAN ACCEL DUE TO GRAVITY AT EARTH'S SURFACE REAL(KIND=r8), PARAMETER :: G = 9.80665_r8 REAL(KIND=r8), PARAMETER :: CP = 1005.0_r8 REAL(KIND=r8), PARAMETER :: R = 287.05_r8 !!---------------------------------------------------------------------- ! C_VKMAN start REAL(KIND=r8),PARAMETER:: VKMAN=0.4_r8 ! Von Karman's constant !----------------------------------------------------------------------- !! Calculate total water content, QW and Liquid water temperature, TL !----------------------------------------------------------------------- ! Derived local parameters. !*---------------------------------------------------------------------- ! C_LHEAT start ! latent heat of condensation of water at 0degc REAL(KIND=r8),PARAMETER:: LC=2.501E6_r8 ! latent heat of fusion at 0degc REAL(KIND=r8),PARAMETER:: LF=0.334E6_r8 REAL(KIND=r8), PARAMETER :: tice=271.16_r8 REAL(KIND=r8),PARAMETER:: RHOSEA = 1000.0_r8 ! density of sea water (kg/m3) REAL(KIND=r8),PARAMETER :: LCRCP=LC/CP ! Evaporation-to-dT conversion factor. REAL(KIND=r8),PARAMETER :: LS=LF+LC ! Latent heat of sublimation. REAL(KIND=r8),PARAMETER :: LSRCP=LS/CP ! Sublimation-to-dT conversion factor. REAL(KIND=r8) :: USHEAR,VSHEAR,VSHR2 REAL(KIND=r8) :: rbyg rbyg =R/G*del(1)*0.5_r8 rmi =0.0_r8 rhi =0.0_r8 evap =0.0_r8 sens =0.0_r8 ustar=0.0_r8 DO K=1,kMAx DO I=1,nCols IF(mskant(i) == 1_i8)THEN P(i,k) =ps(i)*sl(k)*100.0_r8 QW(I,K) = Q(I,K) + QCL(I,K) + QCF(I,K) ! P243.10 TL(I,K) = T(I,K) - LCRCP*QCL(I,K) - LSRCP*QCF(I,K)! P243.9 END IF ENDDO ENDDO ICE_FRACT=0.0_r8 DO i = 1, ncols IF(mskant(i) == 1_i8)THEN IF (TSTAR_SEA(i) < 0.0_r8 .AND. ABS(TSTAR_SEA(i)) < tice+0.01_r8) THEN ICE_FRACT(i)=1.0_r8 ELSE IF (TSTAR_SEA(i) < 0.0_r8 .AND. ABS(TSTAR_SEA(i)) > tice+0.01_r8) THEN ICE_FRACT(i)=0.0_r8 END IF END IF END DO FLANDG=0.0_r8 DO I=1,nCols IF(mskant(i) == 1_i8)THEN IF(IMASK(i) > 1_i8)FLANDG(i)=1.0_r8 QW_1(i)=QW(I,1) TL_1(i)=TL(I,1) Z1_TQ(i)=rbyg*T(i,1) Z1_UV(i)=Z1_TQ(i) Pstar(i) =ps(i)*si(1)*100.0_r8 speedm(i)=SQRT(gu(i,1)**2+gv(i,1)**2)*sinclt(i) USHEAR = ((gu(I,1) - 0.01_r8*ABS(gu(I,1)))*sinclt(i)) VSHEAR = ((gv(I,1) - 0.01_r8*ABS(gv(I,1)))*sinclt(i)) !USHEAR = ((gu(I,1) - ((gu(I,1)/MAX(ABS(gu(I,1)),0.1_r8))*sqrt(ABS(gu(I,1)))))*sinclt(i))/Z1_UV(i) !VSHEAR = ((gv(I,1) - ((gv(I,1)/MAX(ABS(gv(I,1)),0.1_r8))*sqrt(ABS(gv(I,1)))))*sinclt(i))/Z1_UV(i) !USHEAR = VKMAN*(gu(I,1) - ((gu(I,1)/MAX(ABS(gu(I,1)),0.1_r8))*sqrt(ABS(gu(I,1)))))* & ! sinclt(i)/(log(Z1_UV(i)/Z0MSEA(i))) !VSHEAR = VKMAN*(gv(I,1) - ((gv(I,1)/MAX(ABS(gv(I,1)),0.1_r8))*sqrt(ABS(gv(I,1)))))* & ! sinclt(i)/(log(Z1_UV(i)/Z0MSEA(i))) VSHR2 = MAX (1.0E-6_r8 , USHEAR*USHEAR + VSHEAR*VSHEAR) VSHR_SSI(I) = SQRT(VSHR2) END IF END DO CALL BOUY_TQ (& nCols , & !INTEGER, INTENT(IN ) :: nCols kMAx , & !INTEGER, INTENT(IN ) :: kMAx LQ_MIX_BL , & !LOGICAL, INTENT(IN ) :: LQ_MIX_BL mskant (1:nCols ) , & P (1:nCols,1:kMAx), & !REAL, INTENT(IN ) :: P (nCols,rows,kMAx) T (1:nCols,1:kMAx), & !REAL, INTENT(IN ) :: T (nCols,rows,kMAx) Q (1:nCols,1:kMAx), & !REAL, INTENT(IN ) :: Q (nCols,rows,kMAx) QCF (1:nCols,1:kMAx), & !REAL(KIND=r8)(KIND=r8), INTENT(IN ) :: QCF (nCols,rows,kMAx) QCL (1:nCols,1:kMAx), & !REAL(KIND=r8)(KIND=r8), INTENT(IN ) :: QCL (nCols,rows,kMAx) CF (1:nCols,1:kMAx), & !REAL(KIND=r8)(KIND=r8), INTENT(IN ) :: CF (nCols,rows,kMAx) BT (1:nCols,1:kMAx), & !REAL(KIND=r8)(KIND=r8), INTENT(out ) :: BT (nCols,rows,kMAx) BQ (1:nCols,1:kMAx), & !REAL(KIND=r8)(KIND=r8), INTENT(out ) :: BQ (nCols,rows,kMAx) BT_CLD (1:nCols,1:kMAx), & !REAL(KIND=r8)(KIND=r8), INTENT(out ) :: BT_CLD (nCols,rows,kMAx) BQ_CLD (1:nCols,1:kMAx), & !REAL(KIND=r8)(KIND=r8), INTENT(out ) :: BQ_CLD (nCols,rows,kMAx) BT_GB (1:nCols,1:kMAx), & !REAL(KIND=r8)(KIND=r8), INTENT(out ) :: BT_GB (nCols,rows,kMAx) BQ_GB (1:nCols,1:kMAx), & !REAL(KIND=r8)(KIND=r8), INTENT(out ) :: BQ_GB (nCols,rows,kMAx) A_QS (1:nCols,1:kMAx), & !REAL(KIND=r8), INTENT(out ) :: A_QS (nCols,rows,kMAx) A_DQSDT (1:nCols,1:kMAx), & !REAL(KIND=r8), INTENT(out ) :: A_DQSDT(nCols,rows,kMAx) DQSDT (1:nCols,1:kMAx) & !REAL(KIND=r8), INTENT(out ) :: DQSDT (nCols,rows,kMAx) ) DO I=1,nCols IF(mskant(i) == 1_i8)THEN BQ_1(i)=BQ_GB (i,1) BT_1(i)=BT_GB (I,1) END IF END DO ! IF (LTIMER) THEN ! DEPENDS ON: timer ! CALL TIMER('SFEXCH ',3) ! ENDIF !--------- &,VSHR_SSI(nCols,ROWS) & ! ! IN Mag. of mean sea sfc-to-lowest-level ! ! wind shear !-------------------------------------------------------------- !! 1. Index array for sea-ice !----------------------------------------------------------------------- NSICE = 0 DO I=1,nCols IF(mskant(i) == 1_i8)THEN IF ( ICE_FRACT(I) > 0.0_r8 .AND. FLANDG(I) < 1.0_r8 ) THEN NSICE = NSICE + 1 SICE_INDEX(NSICE,1) = I ENDIF END IF ENDDO !----------------------------------------------------------------------- !! 1.1 Calculate height of lowest model level relative to sea. !----------------------------------------------------------------------- DO I=1,nCols IF(mskant(i) == 1_i8)THEN Z1_TQ_SEA(I)=Z1_TQ(I) !IF(L_CTILE.AND.FLANDG(I) > 0.0.AND.FLANDG(I) < 1.0) Z1_TQ_SEA(I)=Z1_TQ(I)+Z_LAND(I) !IF(FLANDG(I) > 0.0.AND.FLANDG(I) < 1.0) Z1_TQ_SEA(I)=Z1_TQ(I)+Z_LAND(I) END IF ENDDO !----------------------------------------------------------------------- !! 2. Calculate QSAT values required later. !----------------------------------------------------------------------- DO I=1,nCols IF(mskant(i) == 1_i8)THEN RHOSTAR(I) = PSTAR(I) /( R*(FLANDG(I)*TSTAR_LAND(I)+(1.0_r8-FLANDG(I))*TSTAR_SSI(I)) ) ! ... surface air density from ideal gas equation END IF ENDDO ! DEPENDS ON: qsat_mix IF(QSAT_opt)THEN CALL QSAT_mix( QS1 (1:nCols) , &!REAL(KIND=r8), intent(out) :: QmixS(npnts) ! Output Saturation mixing ratio or saturationeing processed by qSAT scheme. TL_1 (1:nCols) , &!REAL(KIND=r8), intent(in) :: T(npnts) ! Temperature (K). PSTAR (1:nCols) , &!REAL(KIND=r8), intent(in) :: P(npnts) ! Pressure (Pa). nCols , &!Integer, intent(in) :: npnts !Points (=horizontal dimensions) being processed by qSAT scheme. lq_mix_bl ,&!logical, intent(in) :: lq_mix .true. return qsat as a mixing ratio mskant (1:nCols) ) ! .false. return qsat as a specific humidity ELSE DO i=1,nCols IF(mskant(i) == 1_i8)THEN qss(i) = fpvs(TL_1(i)) QS1(I) = EPS * qss(I) / (PSTAR(I) + EPSM1 * qss(I)) END IF END DO END IF ! DEPENDS ON: qsat_mix IF(QSAT_opt)THEN CALL QSAT_mix( QSTAR_SEA (1:nCols), &!REAL(KIND=r8), intent(out) :: QmixS(npnts) ! Output Saturation mixing ratio or saturationeing processed by qSAT scheme. TSTAR_SEA (1:nCols), &!REAL(KIND=r8), intent(in) :: T(npnts)! Temperature (K). PSTAR (1:nCols), &!REAL(KIND=r8), intent(in) :: P(npnts)! Pressure (Pa). nCols , &!Integer, intent(in) :: npnts !Points (=horizontal dimensions) being processed by qSAT scheme. lq_mix_bl ,&!logical, intent(in) :: lq_mix .true. return qsat as a mixing ratio mskant (1:nCols) ) ! .false. return qsat as a specific humidity ELSE DO i=1,nCols IF(mskant(i) == 1_i8)THEN qss(i) = fpvs(TSTAR_SEA(i)) QSTAR_SEA(I) = EPS * qss(I) / (PSTAR(I) + EPSM1 * qss(I)) END IF END DO END IF ! DEPENDS ON: qsat_mix IF(QSAT_opt)THEN CALL QSAT_mix (QSTAR_ICE (1:nCols), &!REAL(KIND=r8), intent(out) :: QmixS(npnts) ! Output Saturation mixing ratio or saturationeing processed by qSAT scheme. TSTAR_SICE(1:nCols), &!REAL(KIND=r8), intent(in) :: T(npnts)! Temperature (K). PSTAR (1:nCols), &!REAL(KIND=r8), intent(in) :: P(npnts)! Pressure (Pa). nCols , &!Integer, intent(in) :: npnts !Points (=horizontal dimensions) being processed by qSAT scheme. lq_mix_bl ,&!logical, intent(in) :: lq_mix .true. return qsat as a mixing ratio mskant (1:nCols) ) ! .false. return qsat as a specific humidity ELSE DO i=1,nCols IF(mskant(i) == 1_i8)THEN qss(i) = fpvs(TSTAR_SICE(i)) QSTAR_ICE(I) = EPS * qss(I) / (PSTAR(I) + EPSM1 * qss(I)) END IF END DO END IF !----------------------------------------------------------------------- !! 3. Calculation of transfer coefficients and surface layer stability !----------------------------------------------------------------------- !----------------------------------------------------------------------- !! 3.1 Calculate neutral roughness lengths !----------------------------------------------------------------------- ! Sea, sea-ice leads, sea-ice and marginal ice zone DO I=1,nCols IF(mskant(i) == 1_i8)THEN Z0_MIZ(I) = Z0MIZ Z0_ICE(I) = Z0SICE RIB_SEA(I) = 0.0_r8 RIB_ICE(I) = 0.0_r8 DB_SEA(I) = 0.0_r8 DB_ICE(I) = 0.0_r8 END IF ENDDO IF (ISeaZ0T == SurfDivZ0T) THEN ! ! Composite formulation for thermal roughness lengths, ! incoporating the smooth aerodynamic limit for low ! wind speeds and a value based on surface divergence ! theory for higher wind speeds. ! ! The friction velocity is diagnosed in the surface ! transfer scheme, using z0m from the previous time-step. ! z0[T,q] is also required but depends on u_* in this ! scheme. For practical purposes, it is sufficient to ! infer it from u_* determined from z0m, but in general ! a given value of z0m corresponds to two values of ! u_*, so we need to know whether we are on the low or ! high wind speed side of the minimum value of z0m. ! If we are on the high side, z0[T,q] will be inversely ! proportional to z0m, but on the low side it may follow ! this relationship, or be aerodynamically smooth. In ! the smooth case we iterate u_* from the expression for ! the roughness length and then take the maximum of the ! smooth and high-wind expressions for z0[T,q]. An ! iteration for the low-wind value of u_*, USTR_L is ! carried out. This will converge to the correct limit only ! on the low-wind side of the minimum, and the standard ! criterion that the gradient of a fixed-point iteration ! should be less than 1 in modulus gievs a more precise ! boundary, TOL_USTR_L. For consistency with earlier versions ! of the modset, hard-wired values are retained for the ! operational value of Charnock's parameter. An additional ! check is required, since z0m can be large at low or at ! high wind-speeds. This is less precise and a fixed ! value of 0.07 is used to test USTR_N, which was determined ! by inspection of a graph of z0m against u_*: it is ! unlikely that this will need to be changed unless ! Charnock's constant is greatly altered. ! IF (CHARNOCK == 0.018_r8) THEN TOL_USTR_L = 0.055_r8 TOL_USTR_N = 0.07_r8 ELSE TOL_USTR_L = 0.75_r8*(1.54E-6_r8*G/(2.0_r8*CHARNOCK))**0.33333_r8 TOL_USTR_N = 0.07_r8 ENDIF ! DO I=1, nCols IF(mskant(i) == 1_i8)THEN Z0MSEA(I) = MIN(Z0MSEA(I),0.1_r8 * Z1_UV(I)) Z0H_SEA(I) = Z0MSEA(I) ustar(i)=VSHR_SSI(I) ! We need to infer u_* from Z0M. IF (VSHR_SSI(I) > 0.0_r8) THEN ! Compute u_* using neutral stratification. ! stratification. USTR_N = VKMAN * VSHR_SSI(I) / & & LOG(Z1_UV(I) / Z0MSEA(I) ) ustar(i)=USTR_N ! Compute u_* using low wind approximation. USTR_L = 1.54E-06_r8 / Z0MSEA(I) - 1.0E-05_r8 ! Since Z0M could be large for low and high u_*, we use ! USTR_N as an additional check on the regime. IF ( (USTR_N < TOL_USTR_N) .AND. & & (USTR_L < TOL_USTR_L) ) THEN ! Iterate u_* for low winds. DO JITS=1, 5 USTR_L=1.54E-06_r8/(Z0MSEA(I)-(CHARNOCK/G)*USTR_L**2) & & -1.0E-05_r8 ENDDO ! Take the maximum of the smooth and high-wind values. ! A lower limit is imposed on the friction velocity to ! allow for the exceptional case of very low winds: the ! value of 10^-5 is the same as the limit for the momentum ! roughness length. ustar(i)=USTR_L Z0H_SEA(I) = MAX( 2.52E-6_r8/(USTR_L+1.0E-05_r8),2.56E-9_r8/Z0MSEA(I) ) ELSE ! Take the high-wind value, but limit it to the molecular ! mean free path (we should not hit this limit ! in practice). Z0H_SEA(I) = MAX( 2.56E-9_r8/Z0MSEA(I), 7.0E-08_r8 ) ENDIF ELSE IF(ustar(i)==0.0_r8 )ustar(i)=0.000001_r8 ENDIF Z0H_SEA(I) = MIN(Z0H_SEA(I),0.1_r8 * Z1_UV(I)) END IF ENDDO ! ELSE IF (ISeaZ0T == Fixed_Z0T) THEN ! ! Use a fixed thermal roughness length. Z0MSEA(1:nCols) = MIN(Z0MSEA(1:nCols),0.1_r8 * Z1_UV(1:nCols)) Z0H_SEA(1:nCols) = Z0MSEA(1:nCols) !Z0H_SEA(1:nCols) = Z0HSEA ! ENDIF ! !----------------------------------------------------------------------- !! 3.2 Calculate bulk Richardson number for the lowest model level. !----------------------------------------------------------------------- ! Sea, sea-ice and sea-ice leads ! DEPENDS ON: sf_rib_sea CALL SF_RIB_SEA (& nCols , &!INTEGER, INTENT(IN ) :: nCols ! IN Number of X points? mskant (1:nCols) , &!INTEGER, INTENT(IN ) :: nCols ! IN Number of X points? FLANDG (1:nCols) , &!REAL(KIND=r8) , INTENT(IN ) :: FLANDG (nCols,ROWS) ! IN Land fraction on all pts. NSICE , &!INTEGER, INTENT(IN ) :: NSICE ! IN Number of sea-ice points. SICE_INDEX (1:nCols,1:2) , &!INTEGER, INTENT(IN ) :: SICE_INDEX (nCols*ROWS,2) ! IN Index of sea-ice points. BQ_1 (1:nCols) , &!REAL(KIND=r8) , INTENT(IN ) :: BQ_1 (nCols,ROWS) ! IN A buoyancy parameter for lowest atm level. ("beta-q twiddle"). BT_1 (1:nCols) , &!REAL(KIND=r8) , INTENT(IN ) :: BT_1 (nCols,ROWS) ! IN A buoyancy parameter for lowest atm level. ("beta-T twiddle"). QSTAR_ICE (1:nCols) , &!REAL(KIND=r8) , INTENT(IN ) :: QSTAR_ICE (nCols,ROWS) ! IN Surface saturated sp humidity over sea-ice. QSTAR_SEA (1:nCols) , &!REAL(KIND=r8) , INTENT(IN ) :: QSTAR_SEA (nCols,ROWS) ! IN Surface saturated sp humidity over sea and sea-ice leads. QW_1 (1:nCols) , &!REAL(KIND=r8) , INTENT(IN ) :: QW_1 (nCols,ROWS) ! IN Total water content of lowest atmospheric layer (kg per kg air) TL_1 (1:nCols) , &!REAL(KIND=r8) , INTENT(IN ) :: TL_1 (nCols,ROWS) ! IN Liquid/frozen water temperature for lowest atmospheric layer (K). TSTAR_SICE (1:nCols) , &!REAL(KIND=r8) , INTENT(IN ) :: TSTAR_SICE (nCols,ROWS) ! IN Surface temperature of sea-ice (K). TSTAR_SEA (1:nCols) , &!REAL(KIND=r8) , INTENT(IN ) :: TSTAR_SEA (nCols,ROWS) ! IN Surface temperature of sea and sea-ice leads (K). VSHR_SSI (1:nCols) , &!REAL(KIND=r8) , INTENT(IN ) :: VSHR (nCols,ROWS) ! IN Magnitude of surface- to-lowest-level wind shear. Z0_ICE (1:nCols) , &!REAL(KIND=r8) , INTENT(IN ) :: Z0H_ICE (nCols,ROWS) ! IN Roughness length for heat and moisture transport over sea-ice (m). Z0H_SEA (1:nCols) , &!REAL(KIND=r8) , INTENT(IN ) :: Z0H_SEA (nCols,ROWS) ! IN Roughness length for heat and moisture transport over sea or sea-ice Z0_ICE (1:nCols) , &!REAL(KIND=r8) , INTENT(IN ) :: Z0M_ICE (nCols,ROWS) ! IN Roughness length for momentum over sea-ice (m). Z0MSEA (1:nCols) , &!REAL(KIND=r8) , INTENT(IN ) :: Z0M_SEA (nCols,ROWS) ! IN Roughness length for momentum over sea or sea-ice leads (m). Z1_TQ_SEA (1:nCols) , &!REAL(KIND=r8) , INTENT(IN ) :: Z1_TQ (nCols,ROWS) ! IN Height of lowest TQ level (m). Z1_UV (1:nCols) , &!REAL(KIND=r8) , INTENT(IN ) :: Z1_UV (nCols,ROWS) ! IN Height of lowest UV level (m). RIB_SEA (1:nCols) , &!REAL(KIND=r8) , INTENT(OUT ) :: RIB_SEA (nCols,ROWS) ! OUT Bulk Richardson number for lowest layer over sea or sea-ice leads. RIB_ICE (1:nCols) , &!REAL(KIND=r8) , INTENT(OUT ) :: RIB_ICE (nCols,ROWS) ! OUT Bulk Richardson number for lowest layer over sea-ice. DB_SEA (1:nCols) , &!REAL(KIND=r8) , INTENT(OUT ) :: DB_SEA (nCols,ROWS) ! OUT Buoyancy difference between DB_ICE (1:nCols) &!REAL(KIND=r8) , INTENT(OUT ) :: DB_ICE (nCols,ROWS) ! OUT Buoyancy difference between ) !----------------------------------------------------------------------- !! 3.4 Calculate CD, CH via routine FCDCH. !! Note that these are returned as the dimensionless surface !! exchange coefficients. !----------------------------------------------------------------------- ! Sea-ice ! DEPENDS ON: fcdch_sea CALL FCDCH_SEA( & nCols , &!INTEGER, INTENT(IN ) :: nCols ! IN Number of X points? mskant (1:nCols), &!INTEGER, INTENT(IN ) :: nCols ! IN Number of X points? FLANDG (1:nCols), &!REAL(KIND=r8) , INTENT(IN ) :: FLANDG (nCols,ROWS) ! IN Land fraction DB_ICE (1:nCols), &!REAL(KIND=r8) , INTENT(IN ) :: DB (nCols,ROWS) ! VSHR_SSI (1:nCols), &!REAL(KIND=r8) , INTENT(IN ) :: VSHR (nCols,ROWS) Z0_ICE (1:nCols), &!REAL(KIND=r8) , INTENT(IN ) :: Z0M (nCols,ROWS) Z0_ICE (1:nCols), &!REAL(KIND=r8) , INTENT(IN ) :: Z0H (nCols,ROWS) ZH (1:nCols), &!REAL(KIND=r8) , INTENT(IN ) :: ZH (nCols,i) Z1_UV (1:nCols), &!REAL(KIND=r8) , INTENT(IN ) :: Z1_UV (nCols,ROWS) Z1_TQ_SEA (1:nCols), &!REAL(KIND=r8) , INTENT(IN ) :: Z1_TQ (nCols,ROWS) CD_ICE (1:nCols), &!REAL(KIND=r8) , INTENT(OUT ) :: CDV (nCols,ROWS) CH_ICE (1:nCols), &!REAL(KIND=r8) , INTENT(OUT ) :: CHV (nCols,ROWS) V_S_ICE (1:nCols), &!REAL(KIND=r8) , INTENT(OUT ) :: V_S (nCols,ROWS) RECIP_L_MO_ICE(1:nCols) &!REAL(KIND=r8) , INTENT(OUT ) :: RECIP_L_MO (nCols,ROWS) ) ! Marginal Ice Zone ! DEPENDS ON: fcdch_sea CALL FCDCH_SEA(& nCols , &!INTEGER, INTENT(IN ) :: nCols ! IN Number of X points? mskant (1:nCols), &!INTEGER, INTENT(IN ) :: nCols ! IN Number of X points? FLANDG (1:nCols), &!REAL(KIND=r8) , INTENT(IN ) :: FLANDG (nCols,ROWS) ! IN Land fraction DB_ICE (1:nCols), &!REAL(KIND=r8) , INTENT(IN ) :: DB (nCols,ROWS) ! VSHR_SSI (1:nCols), &!REAL(KIND=r8) , INTENT(IN ) :: VSHR (nCols,ROWS) Z0_MIZ (1:nCols), &!REAL(KIND=r8) , INTENT(IN ) :: Z0M (nCols,ROWS) Z0_MIZ (1:nCols), &!REAL(KIND=r8) , INTENT(IN ) :: Z0H (nCols,ROWS) ZH (1:nCols), &!REAL(KIND=r8) , INTENT(IN ) :: ZH (nCols,ROWS) Z1_UV (1:nCols), &!REAL(KIND=r8) , INTENT(IN ) :: Z1_UV (nCols,ROWS) Z1_TQ_SEA (1:nCols), &!REAL(KIND=r8) , INTENT(IN ) :: Z1_TQ (nCols,ROWS) CD_MIZ (1:nCols), &!REAL(KIND=r8) , INTENT(OUT ) :: CDV (nCols,ROWS) CH_MIZ (1:nCols), &!REAL(KIND=r8) , INTENT(OUT ) :: CHV (nCols,ROWS) V_S_MIZ (1:nCols), &!REAL(KIND=r8) , INTENT(OUT ) :: V_S (nCols,ROWS) RECIP_L_MO_MIZ(1:nCols) &!REAL(KIND=r8) , INTENT(OUT ) :: RECIP_L_MO (nCols,ROWS) ) ! Sea and sea-ice leads ! DEPENDS ON: fcdch_sea CALL FCDCH_SEA(& nCols , & !INTEGER, INTENT(IN ) :: nCols! IN Number of X points? mskant (1:nCols), &!INTEGER, INTENT(IN ) :: nCols ! IN Number of X points? FLANDG (1:nCols), & !REAL(KIND=r8) , INTENT(IN ) :: FLANDG(nCols,ROWS) ! IN Land fraction DB_SEA (1:nCols), & !REAL(KIND=r8) , INTENT(IN ) :: DB (nCols,ROWS) ! VSHR_SSI (1:nCols), & !REAL(KIND=r8) , INTENT(IN ) :: VSHR (nCols,ROWS) Z0MSEA (1:nCols), & !REAL(KIND=r8) , INTENT(IN ) :: Z0M (nCols,ROWS) Z0H_SEA (1:nCols), & !REAL(KIND=r8) , INTENT(IN ) :: Z0H (nCols,ROWS) ZH (1:nCols), & !REAL(KIND=r8) , INTENT(IN ) :: ZH (nCols,ROWS) Z1_UV (1:nCols), & !REAL(KIND=r8) , INTENT(IN ) :: Z1_UV(nCols,ROWS) Z1_TQ_SEA (1:nCols), & !REAL(KIND=r8) , INTENT(IN ) :: Z1_TQ(nCols,ROWS) CD_SEA (1:nCols), & !REAL(KIND=r8) , INTENT(OUT ) :: CDV (nCols,ROWS) CH_SEA (1:nCols), & !REAL(KIND=r8) , INTENT(OUT ) :: CHV (nCols,ROWS) V_S_SEA (1:nCols), & !REAL(KIND=r8) , INTENT(OUT ) :: V_S(nCols,ROWS) RECIP_L_MO_SEA(1:nCols) & !REAL(KIND=r8) , INTENT(OUT ) :: RECIP_L_MO(nCols,ROWS) ) !----------------------------------------------------------------------- ! Calculate gridbox-means of transfer coefficients. !----------------------------------------------------------------------- DO I=1,nCols IF(mskant(i) == 1_i8)THEN CD_LAND(I) = 0.0_r8 CD_SSI(I) = 0.0_r8 CH_SSI(I) = 0.0_r8 CD(I) = 0.0_r8 CH(I) = 0.0_r8 IF (L_DUST) CD_STD_DUST(I) = 0.0_r8 ! Sea and sea-ice IF (FLANDG(I) < 1.0_r8 ) THEN IF ( ICE_FRACT(I) < 0.7_r8 ) THEN CD_SSI(I) = ( ICE_FRACT(I)*CD_MIZ(I) + (0.7_r8-ICE_FRACT(I))*CD_SEA(I) ) / 0.7_r8 ! P2430.5 CH_SSI(I) = ( ICE_FRACT(I)*CH_MIZ(I) + (0.7_r8-ICE_FRACT(I))*CH_SEA(I) ) / 0.7_r8 ! P2430.4 ELSE CD_SSI(I) = ( (1.0_r8-ICE_FRACT(I))*CD_MIZ(I) + (ICE_FRACT(I)-0.7_r8)*CD_ICE(I) ) / 0.3_r8 ! P2430.7 CH_SSI(I) = ( (1.0_r8-ICE_FRACT(I))*CH_MIZ(I) + (ICE_FRACT(I)-0.7_r8)*CH_ICE(I) ) / 0.3_r8 ! P2430.7 ENDIF CD(I)=(1.0_r8-FLANDG(I))*CD_SSI(I) CH(I)=(1.0_r8-FLANDG(I))*CH_SSI(I) IF (L_DUST) CD_STD_DUST(I)=(1.0_r8-FLANDG(I))*CD_SSI(I) ENDIF END IF ENDDO ! Land tiles ! DO N=1,NTILES !CDIR NODEP ! DO K=1,TILE_PTS(N) ! L = TILE_INDEX(K,N) ! J=(LAND_INDEX(L)-1)/nCols + 1 ! I = LAND_INDEX(L) - (J-1)*nCols ! CD_LAND(I,J) = CD_LAND(I,J) + TILE_FRAC(L,N)*CD_TILE(L,N) ! CD(I,J) = CD(I,J) + FLANDG(I,J)*TILE_FRAC(L,N)*CD_TILE(L,N) ! CH(I,J) = CH(I,J) + FLANDG(I,J)*TILE_FRAC(L,N)*CH_TILE(L,N) ! IF (L_DUST) CD_STD_DUST(I,J) = CD_STD_DUST(I,J) + & ! & FLANDG(I,J)*TILE_FRAC(L,N)*CD_STD(L,N)! ! ! ENDDO ! ENDDO !----------------------------------------------------------------------- !! 4.3 Calculate the surface exchange coefficients RHOK(*) and ! resistances for use in Sulphur Cycle ! (Note that CD_STD, CH and VSHR should never = 0) ! RHOSTAR * CD * VSHR stored for diagnostic output before ! horizontal interpolation. !----------------------------------------------------------------------- DO I=1,nCols IF(mskant(i) == 1_i8)THEN RHO_ARESIST(I) = 0.0_r8 ! RHO_ARESIST_LAND(I)=0.0 ARESIST(I) = 0.0_r8 RESIST_B(I) = 0.0_r8 ! RHOKM_LAND(I) = 0. RHOKM_SSI(I) = 0.0_r8 ! Sea and sea-ice IF ( FLANDG(I) < 1.0_r8 ) THEN RHOKM_SSI(I) = RHOSTAR(I)*CD_SSI(I)*VSHR_SSI(I) ! ! P243.124 RHOKH_1_SICE(I) = RHOSTAR(I)*CH_SSI(I)*VSHR_SSI(I) ! ! P243.125 RHO_ARESIST(I) = RHOSTAR(I)*CD_SSI(I)*VSHR_SSI(I) ARESIST (I) = 1.0_r8 / (CD_SSI(I) * VSHR_SSI(I)) RESIST_B(I)= (CD_SSI(I)/CH_SSI(I) - 1.0_r8) * ARESIST(I) ENDIF END IF ENDDO !----------------------------------------------------------------------- !! Calculate surface layer resistance for mineral dust !----------------------------------------------------------------------- IF (L_DUST) THEN DO I = 1,nCols IF(mskant(i) == 1_i8)THEN DO IDIV=1,NDIV R_B_DUST(I,IDIV)=0.0_r8 ENDDO !IDIV VSHR(I)=(1.0_r8 - FLANDG(I)) * VSHR_SSI(I) !+ FLANDG(I) * VSHR_LAND(I) END IF ENDDO !I ! DEPENDS ON: dustresb CALL DUSTRESB (& nCols , &!INTEGER,INTENT(IN ) :: nCols !IN mskant, & PSTAR , &!REAL(KIND=r8), INTENT(IN ) :: PSTAR(nCols,ROWS) !IN surface pressure TSTAR_SEA , &!REAL(KIND=r8), INTENT(IN ) :: TSTAR(nCols,ROWS) !IN surface temperature RHOSTAR , &!REAL(KIND=r8), INTENT(IN ) :: RHOSTAR(nCols,ROWS)!IN surface air density VSHR , &!REAL(KIND=r8), INTENT(IN ) :: VSHR(nCols,ROWS) !IN surface to lowest lev windspeed difference CD_STD_DUST , &!REAL(KIND=r8), INTENT(IN ) :: CD_STD_DUST(nCols,ROWS) !IN surface transfer coeffient for R_B_DUST &! REAL(KIND=r8), INTENT(OUT ) :: R_B_DUST(nCols,ROWS,NDIV) !OUT surface layer resistance for ) ENDIF !(L_DUST) !----------------------------------------------------------------------- !! Calculate local and gridbox-average surface fluxes of heat and !! moisture. !----------------------------------------------------------------------- DO I=1,nCols IF(mskant(i) == 1_i8)THEN FTL_1(I) = 0.0_r8 FQW_1(I) = 0.0_r8 END IF ENDDO ! Sea and sea-ice ! DEPENDS ON: sf_flux_sea CALL SF_FLUX_SEA ( & nCols , &!INTEGER, INTENT(IN ) :: nCols ! IN Number of X points? NSICE , &!INTEGER, INTENT(IN ) :: NSICE ! IN Number of sea-ice points. mskant , & SICE_INDEX (1:nCols,1:2) , &!INTEGER, INTENT(IN ) :: SICE_INDEX(nCols*ROWS,2)! IN Index of sea-ice points FLANDG (1:nCols) , &!REAL(KIND=r8) , INTENT(IN ) :: FLANDG (nCols,ROWS) ICE_FRACT (1:nCols) , &!REAL(KIND=r8) , INTENT(IN ) :: ICE_FRACT (nCols,ROWS) ! IN Fraction of gridbox which is sea-ice. QS1 (1:nCols) , &!REAL(KIND=r8) , INTENT(IN ) :: QS1 (nCols,ROWS) ! IN Sat. specific humidity qsat(TL_1,PSTAR) QSTAR_ICE (1:nCols) , &!REAL(KIND=r8) , INTENT(IN ) :: QSTAR_ICE (nCols,ROWS) ! IN Surface qsat for sea-ice. QSTAR_SEA (1:nCols) , &!REAL(KIND=r8) , INTENT(IN ) :: QSTAR_SEA (nCols,ROWS) ! IN Surface qsat for sea or sea-ice leads. QW_1 (1:nCols) , &!REAL(KIND=r8) , INTENT(IN ) :: QW_1 (nCols,ROWS) ! IN Total water content of lowest RADNET (1:nCols) , &!REAL(KIND=r8) , INTENT(IN ) :: RADNET (nCols,ROWS) ! IN Net surface radiation (W/m2)positive downwards RHOKH_1_SICE (1:nCols) , &!REAL(KIND=r8) , INTENT(IN ) :: RHOKH_1 (nCols,ROWS) ! IN Surface exchange coefficient. TI (1:nCols) , &!REAL(KIND=r8) , INTENT(IN ) :: TI (nCols,ROWS) ! IN Temperature of sea-ice surfacelayer (K) TL_1 (1:nCols) , &!REAL(KIND=r8) , INTENT(IN ) :: TL_1 (nCols,ROWS) ! IN Liquid/frozen water temperature for lowest atmospheric layer (K) TSTAR_SICE (1:nCols) , &!REAL(KIND=r8) , INTENT(IN ) :: TSTAR_SICE(nCols,ROWS) ! IN Sea-ice surface temperature (K) TSTAR_SEA (1:nCols) , &!REAL(KIND=r8) , INTENT(IN ) :: TSTAR_SEA (nCols,ROWS) ! IN Sea surface temperature (K). Z0_ICE (1:nCols) , &!REAL(KIND=r8) , INTENT(IN ) :: Z0H_ICE (nCols,ROWS) ! IN Sea-ice heat and moisture Z0_ICE (1:nCols) , &!REAL(KIND=r8) , INTENT(IN ) :: Z0M_ICE (nCols,ROWS) ! IN Sea-ice momentum roughness Z0H_SEA (1:nCols) , &!REAL(KIND=r8) , INTENT(IN ) :: Z0H_SEA (nCols,ROWS) ! IN Sea and lead heat and moisture Z0MSEA (1:nCols) , &!REAL(KIND=r8) , INTENT(IN ) :: Z0M_SEA (nCols,ROWS) ! IN Sea and lead momentum roughness Z1_TQ_SEA (1:nCols) , &!REAL(KIND=r8) , INTENT(IN ) :: Z1_TQ (nCols,ROWS) ! IN Height of lowest atmospheric SeaSalinityFactor , &!REAL(KIND=r8) , Intent(IN ) :: SeaSalinityFactor ALPHA1_SICE (1:nCols) , &!REAL(KIND=r8) , INTENT(OUT ) :: ALPHA1 (nCols,ROWS) ASHTF (1:nCols) , &!REAL(KIND=r8) , INTENT(OUT ) :: ASHTF (nCols,ROWS) E_SEA (1:nCols) , &!REAL(KIND=r8) , INTENT(OUT ) :: E_SEA (nCols,ROWS) FQW_ICE (1:nCols) , &!REAL(KIND=r8) , INTENT(OUT ) :: FQW_ICE (nCols,ROWS) FQW_1 (1:nCols) , &!REAL(KIND=r8) , INTENT(OUT ) :: FQW_1 (nCols,ROWS) FTL_ICE (1:nCols) , &!REAL(KIND=r8) , INTENT(OUT ) :: FTL_ICE (nCols,ROWS) FTL_1 (1:nCols) , &!REAL(KIND=r8) , INTENT(OUT ) :: FTL_1 (nCols,ROWS) H_SEA (1:nCols) , &!REAL(KIND=r8) , INTENT(OUT ) :: H_SEA (nCols,ROWS) RHOKPM_SICE (1:nCols) &!REAL(KIND=r8) , INTENT(OUT ) :: RHOKPM (nCols,ROWS) ! OUT Modified surface exchange ) DO I=1,nCols IF(mskant(i) == 1_i8)THEN rmi (i)=CD_SSI(I)*VSHR_SSI(I) rhi (i)=CH_SSI(I)*VSHR_SSI(I) evap (i)=FQW_1(i)*hltm sens (i)=FTL_1(i)*CP END IF ENDDO !----------------------------------------------------------------------- !! 4.4 Calculate the standard deviations of layer 1 turbulent !! fluctuations of temperature and humidity using approximate !! formulae from first order closure. !----------------------------------------------------------------------- DO I=1,nCols IF(mskant(i) == 1_i8)THEN Q1_SD(I) = 0.0_r8 T1_SD(I) = 0.0_r8 END IF ENDDO ! Sea and sea-ice ! DEPENDS ON: stdev1_sea CALL STDEV1_SEA (& nCols , & !INTEGER, INTENT(IN ) :: nCols ! IN Number of X points? mskant , & FLANDG (1:nCols), & !REAL(KIND=r8) , INTENT(IN ) :: FLANDG (nCols,ROWS) ! IN Land fraction. BQ_1 (1:nCols), & !REAL(KIND=r8) , INTENT(IN ) :: BQ_1 (nCols,ROWS) ! IN Buoyancy parameter. BT_1 (1:nCols), & !REAL(KIND=r8) , INTENT(IN ) :: BT_1 (nCols,ROWS) ! IN Buoyancy parameter. FQW_1 (1:nCols), & !REAL(KIND=r8) , INTENT(IN ) :: FQW_1 (nCols,ROWS) ! IN Surface flux of QW. FTL_1 (1:nCols), & !REAL(KIND=r8) , INTENT(IN ) :: FTL_1 (nCols,ROWS) ! IN Surface flux of TL. ICE_FRACT (1:nCols), & !REAL(KIND=r8) , INTENT(IN ) :: ICE_FRACT(nCols,ROWS) ! IN Fraction of gridbox which is sea-ice. RHOKM_SSI (1:nCols), & !REAL(KIND=r8) , INTENT(IN ) :: RHOKM_1 (nCols,ROWS) ! IN Surface momentum exchange RHOSTAR (1:nCols), & !REAL(KIND=r8) , INTENT(IN ) :: RHOSTAR (nCols,ROWS) ! IN Surface air density. VSHR_SSI (1:nCols), & !REAL(KIND=r8) , INTENT(IN ) :: VSHR (nCols,ROWS) ! IN Magnitude of surface- Z0MSEA (1:nCols), & !REAL(KIND=r8) , INTENT(IN ) :: Z0MSEA (nCols,ROWS) ! IN Sea roughness length. Z0_ICE (1:nCols), & !REAL(KIND=r8) , INTENT(IN ) :: Z0_ICE (nCols,ROWS) ! IN Sea-ice roughness length. Z1_TQ_SEA (1:nCols), & !REAL(KIND=r8) , INTENT(IN ) :: Z1_TQ (nCols,ROWS) ! IN Height of lowest tq level. Q1_SD (1:nCols), & !REAL(KIND=r8) , INTENT(out ) :: Q1_SD (nCols,ROWS) ! OUT Standard deviation of T1_SD (1:nCols) ) !REAL(KIND=r8) , INTENT(out ) :: T1_SD (nCols,ROWS) ! OUT Standard deviation of !----------------------------------------------------------------------- !! 4.6 For sea points, calculate the wind mixing energy flux and the !! sea-surface roughness length on the P-grid, using time-level n !! quantities. !----------------------------------------------------------------------- DO I=1,nCols IF(mskant(i) == 1_i8)THEN IF (SFME) FME(I) = 0.0_r8 IF (FLANDG(I) < 1.0_r8) THEN TAU = RHOKM_SSI(I) * VSHR_SSI(I) ! P243.130 IF (ICE_FRACT(I) > 0.0_r8) & & TAU = RHOSTAR(I) * CD_SEA(I) * VSHR_SSI(I) * VSHR_SSI(I) IF (SFME)FME(I) = (1.0_r8-ICE_FRACT(I)) * TAU * SQRT(TAU/RHOSEA) ! ! P243.96 ! Limit Z0MSEA to 0.154m for TAU very small Z0MSEA(I) = 1.54E-6_r8 / (SQRT(TAU/RHOSTAR(I)) + 1.0E-5_r8) & & + (CHARNOCK/G) * (TAU / RHOSTAR(I)) Z0MSEA(I) = MAX ( Z0HSEA , Z0MSEA(I) ) ! ... P243.B6 (Charnock formula) ! TAU/RHOSTAR is "mod VS squared", see eqn P243.131 ENDIF END IF ENDDO END SUBROUTINE SF_EXCH ! *****************************COPYRIGHT******************************* ! (C) Crown copyright Met Office. All rights reserved. ! For further details please refer to the file COPYRIGHT.txt ! which you should have received as part of this distribution. ! *****************************COPYRIGHT******************************* !----------------------------------------------------------------------- ! ! Subroutines SF_FLUX_LAND and SF_FLUX_SEA to calculate explicit surface ! fluxes of heat and moisture ! ! ! Model Modification history ! version date ! 5.2 15/11/00 New Deck M. Best ! 5.3 25/04/01 Add coastal tiling. Nic Gedney ! 6.1 01/09/04 Calculate potential evaporation related variables. ! Nic Gedney ! 6.2 07/11/05 Allow for the salinity of sea water in ! the evaporative flux. ! J. M. Edwards ! ! Programming standard: ! !----------------------------------------------------------------------- ! SUBROUTINE SF_FLUX_LAND------------------------------------------- ! ! Calculate explicit surface fluxes of heat and moisture over ! land tiles ! ! ------------------------------------------------------------------ ! SUBROUTINE SF_FLUX_SEA-------------------------------------------- ! ! Calculate explicit surface fluxes of heat and moisture over sea ! and sea-ice ! ! ------------------------------------------------------------------ SUBROUTINE SF_FLUX_SEA (& nCols , &!INTEGER, INTENT(IN ) :: nCols ! IN Number of X points? NSICE , &!INTEGER, INTENT(IN ) :: NSICE ! IN Number of sea-ice points. mskant, & SICE_INDEX , &!INTEGER, INTENT(IN ) :: SICE_INDEX(nCols*ROWS,2)! IN Index of sea-ice points FLANDG , &!REAL(KIND=r8), INTENT(IN ) :: FLANDG(nCols,ROWS) ICE_FRACT , &!REAL(KIND=r8), INTENT(IN ) :: ICE_FRACT(nCols,ROWS) ! IN Fraction of gridbox which is sea-ice. QS1 , &!REAL(KIND=r8), INTENT(IN ) ::QS1(nCols,ROWS)! IN Sat. specific humidity qsat(TL_1,PSTAR) QSTAR_ICE , &!REAL(KIND=r8), INTENT(IN ) ::QSTAR_ICE(nCols,ROWS)! IN Surface qsat for sea-ice. QSTAR_SEA , &!REAL(KIND=r8), INTENT(IN ) ::QSTAR_SEA(nCols,ROWS)! IN Surface qsat for sea or sea-ice leads. QW_1 , &!REAL(KIND=r8), INTENT(IN ) ::QW_1(nCols,ROWS)! IN Total water content of lowest RADNET , &!REAL(KIND=r8), INTENT(IN ) ::RADNET(nCols,ROWS)! IN Net surface radiation (W/m2)positive downwards RHOKH_1 , &!REAL(KIND=r8), INTENT(IN ) ::RHOKH_1(nCols,ROWS)! IN Surface exchange coefficient. TI , &!REAL(KIND=r8), INTENT(IN ) ::TI(nCols,ROWS) ! IN Temperature of sea-ice surfacelayer (K) TL_1 , &!REAL(KIND=r8), INTENT(IN ) ::TL_1(nCols,ROWS) ! IN Liquid/frozen water temperature for lowest atmospheric layer (K) TSTAR_SICE , &!REAL(KIND=r8), INTENT(IN ) ::TSTAR_SICE(nCols,ROWS) ! IN Sea-ice surface temperature (K) TSTAR_SEA , &!REAL(KIND=r8), INTENT(IN ) ::TSTAR_SEA(nCols,ROWS) ! IN Sea surface temperature (K). Z0H_ICE , &!REAL(KIND=r8), INTENT(IN ) ::Z0H_ICE(nCols,ROWS)! IN Sea-ice heat and moisture Z0M_ICE , &!REAL(KIND=r8), INTENT(IN ) ::Z0M_ICE(nCols,ROWS)! IN Sea-ice momentum roughness Z0H_SEA , &!REAL(KIND=r8), INTENT(IN ) ::Z0H_SEA(nCols,ROWS) ! IN Sea and lead heat and moisture Z0M_SEA , &!REAL(KIND=r8), INTENT(IN ) ::Z0M_SEA(nCols,ROWS) ! IN Sea and lead momentum roughness Z1_TQ , &!REAL(KIND=r8), INTENT(IN ) ::Z1_TQ(nCols,ROWS)! IN Height of lowest atmospheric SeaSalinityFactor, &!REAL(KIND=r8), Intent(IN) :: SeaSalinityFactor ALPHA1 , &!REAL(KIND=r8) , INTENT(OUT ):: ALPHA1(nCols,ROWS) ASHTF , &!REAL(KIND=r8) , INTENT(OUT ):: ASHTF(nCols,ROWS) E_SEA , &!REAL(KIND=r8) , INTENT(OUT ):: E_SEA(nCols,ROWS) FQW_ICE , &!REAL(KIND=r8) , INTENT(OUT ):: FQW_ICE(nCols,ROWS) FQW_1 , &!REAL(KIND=r8) , INTENT(OUT ):: FQW_1(nCols,ROWS) FTL_ICE , &!REAL(KIND=r8) , INTENT(OUT ):: FTL_ICE(nCols,ROWS) FTL_1 , &!REAL(KIND=r8) , INTENT(OUT ):: FTL_1(nCols,ROWS) H_SEA , &!REAL(KIND=r8) , INTENT(OUT ):: H_SEA(nCols,ROWS) RHOKPM &!REAL(KIND=r8) , INTENT(OUT ):: RHOKPM(nCols,ROWS) ! OUT Modified surface exchange ) IMPLICIT NONE INTEGER, INTENT(IN ) :: nCols ! IN Number of X points? INTEGER, INTENT(IN ) :: NSICE ! IN Number of sea-ice points. INTEGER, INTENT(IN ) :: SICE_INDEX(nCols,2)! IN Index of sea-ice points INTEGER(KIND=i8), INTENT(IN ) :: mskant(nCols) REAL(KIND=r8), INTENT(IN ) :: FLANDG(nCols) REAL(KIND=r8), INTENT(IN ) :: ICE_FRACT(nCols) ! IN Fraction of gridbox which is sea-ice. REAL(KIND=r8), INTENT(IN ) ::QS1(nCols)! IN Sat. specific humidity qsat(TL_1,PSTAR) REAL(KIND=r8), INTENT(IN ) ::QSTAR_ICE(nCols)! IN Surface qsat for sea-ice. REAL(KIND=r8), INTENT(IN ) ::QSTAR_SEA(nCols)! IN Surface qsat for sea or sea-ice leads. REAL(KIND=r8), INTENT(IN ) ::QW_1(nCols)! IN Total water content of lowest ! ! atmospheric layer ! ! (kg per kg air). REAL(KIND=r8), INTENT(IN ) ::RADNET(nCols)! IN Net surface radiation (W/m2) ! ! positive downwards REAL(KIND=r8), INTENT(IN ) ::RHOKH_1(nCols)! IN Surface exchange coefficient. REAL(KIND=r8), INTENT(IN ) ::TI(nCols) ! IN Temperature of sea-ice surface ! ! layer (K) REAL(KIND=r8), INTENT(IN ) ::TL_1(nCols) ! IN Liquid/frozen water temperature ! ! for lowest atmospheric layer (K) REAL(KIND=r8), INTENT(IN ) ::TSTAR_SICE(nCols) ! IN Sea-ice surface temperature (K) REAL(KIND=r8), INTENT(IN ) ::TSTAR_SEA(nCols) ! IN Sea surface temperature (K). REAL(KIND=r8), INTENT(IN ) ::Z0H_ICE(nCols)! IN Sea-ice heat and moisture ! ! roughness length (m). REAL(KIND=r8), INTENT(IN ) ::Z0M_ICE(nCols)! IN Sea-ice momentum roughness ! ! length (m). REAL(KIND=r8), INTENT(IN ) ::Z0H_SEA(nCols) ! IN Sea and lead heat and moisture ! ! roughness length (m). REAL(KIND=r8), INTENT(IN ) ::Z0M_SEA(nCols) ! IN Sea and lead momentum roughness ! ! length. REAL(KIND=r8), INTENT(IN ) ::Z1_TQ(nCols)! IN Height of lowest atmospheric ! ! level (m). ! ! REAL(KIND=r8), INTENT(IN) :: SeaSalinityFactor ! ! Factor allowing for the ! ! effect of the salinity of ! ! sea water on the evaporative ! ! flux. REAL(KIND=r8) , INTENT(OUT ):: ALPHA1(nCols) ! OUT Gradient of saturated specific ! ! humidity with respect to ! ! temperature between the bottom ! ! model layer and the surface. REAL(KIND=r8) , INTENT(OUT ):: ASHTF(nCols) ! OUT Coefficient to calculate ! ! surface heat flux into sea-ice ! ! (W/m2/K). REAL(KIND=r8) , INTENT(OUT ):: E_SEA(nCols) ! OUT Evaporation from sea times ! ! leads fraction (kg/m2/s). REAL(KIND=r8) , INTENT(OUT ):: FQW_ICE(nCols) ! OUT Surface flux of QW for sea-ice. REAL(KIND=r8) , INTENT(OUT ):: FQW_1(nCols) ! OUT GBM surface flux of QW ! ! (kg/m2/s). REAL(KIND=r8) , INTENT(OUT ):: FTL_ICE(nCols) ! OUT Surface flux of TL for sea-ice. REAL(KIND=r8) , INTENT(OUT ):: FTL_1(nCols) ! OUT GBM surface flux of TL. REAL(KIND=r8) , INTENT(OUT ):: H_SEA(nCols)! OUT Surface sensible heat flux over ! ! sea times leads fraction (W/m2) REAL(KIND=r8) , INTENT(OUT ):: RHOKPM(nCols) ! OUT Modified surface exchange ! ! coefficient. !*L------------------COMDECK C_EPSLON----------------------------------- ! EPSILON IS RATIO OF MOLECULAR WEIGHTS OF WATER AND DRY AIR REAL(KIND=r8), PARAMETER :: Epsilon = 0.62198_r8 REAL(KIND=r8), PARAMETER :: C_Virtual = 1.0_r8/Epsilon-1.0_r8 !*---------------------------------------------------------------------- !*L------------------COMDECK C_O_DG_C----------------------------------- ! ZERODEGC IS CONVERSION BETWEEN DEGREES CELSIUS AND KELVIN ! TFS IS TEMPERATURE AT WHICH SEA WATER FREEZES ! TM IS TEMPERATURE AT WHICH FRESH WATER FREEZES AND ICE MELTS REAL(KIND=r8), PARAMETER :: ZeroDegC = 273.15_r8 REAL(KIND=r8), PARAMETER :: TFS = 271.35_r8 REAL(KIND=r8), PARAMETER :: TM = 273.15_r8 !*---------------------------------------------------------------------- !*L------------------COMDECK C_G---------------------------------------- ! G IS MEAN ACCEL DUE TO GRAVITY AT EARTH'S SURFACE REAL(KIND=r8), PARAMETER :: G = 9.80665_r8 !*---------------------------------------------------------------------- ! C_KAPPAI start ! Thermal conductivity of sea-ice (W per m per K). REAL(KIND=r8),PARAMETER:: KAPPAI=2.09_r8 ! Effective thickness of sea-ice surface layer (m). REAL(KIND=r8),PARAMETER:: DE = 0.1_r8 ! C_KAPPAI end ! C_LHEAT start ! latent heat of condensation of water at 0degc REAL(KIND=r8),PARAMETER:: LC=2.501E6_r8 ! latent heat of fusion at 0degc REAL(KIND=r8),PARAMETER:: LF=0.334E6_r8 ! C_LHEAT end !*L------------------COMDECK C_R_CP------------------------------------- ! History: ! Version Date Comment. ! 5.0 07/05/99 Add variable P_zero for consistency with ! conversion to C-P 'C' dynamics grid. R. Rawlins ! 5.1 07/03/00 Fixed/Free format conversion P. Selwood ! R IS GAS CONSTANT FOR DRY AIR ! CP IS SPECIFIC HEAT OF DRY AIR AT CONSTANT PRESSURE ! PREF IS REFERENCE SURFACE PRESSURE REAL(KIND=r8), PARAMETER :: R = 287.05_r8 REAL(KIND=r8), PARAMETER :: CP = 1005.0_r8 REAL(KIND=r8), PARAMETER :: Kappa = R/CP REAL(KIND=r8), PARAMETER :: Pref = 100000.0_r8 ! Reference surface pressure = PREF REAL(KIND=r8), PARAMETER :: P_zero = Pref REAL(KIND=r8), PARAMETER :: sclht = 6.8E+03_r8 ! Scale Height H !*---------------------------------------------------------------------- ! CSIGMA start ! Stefan-Boltzmann constant (W/m**2/K**4). REAL(KIND=r8), PARAMETER :: SBCON=5.67E-8_r8 ! CSIGMA end ! Derived local parameters. REAL(KIND=r8) GRCP,LS PARAMETER ( & & GRCP=G/CP & &,LS=LF+LC & ! Latent heat of sublimation. & ) ! Scalars INTEGER & & I & ! Horizontal field index. &,K ! Sea-ice field index. REAL(KIND=r8) & & DQ1 & ! (qsat(TL_1,PSTAR)-QW_1) + g/cp*alpha1*Z1 &,D_T & ! Temporary in calculation of alpha1. &,RAD_REDUC ! Radiation term required for surface flux ! ! calcs. ! EXTERNAL TIMER ! IF (LTIMER) THEN ! DEPENDS ON: timer ! CALL TIMER('SF_FLUX ',3) ! ENDIF DO I=1,nCols IF(mskant(i) == 1_i8)THEN ALPHA1(I) = 0.0_r8 E_SEA(I) = 0.0_r8 H_SEA(I) = 0.0_r8 FQW_ICE(I) = 0.0_r8 FTL_ICE(I) = 0.0_r8 RHOKPM(I) = 0.0_r8 END IF ENDDO !---------------------------------------------------------------------- !! 1 Calculate gradient of saturated specific humidity for use in !! calculation of surface fluxes - only required for sea-ice points !---------------------------------------------------------------------- DO K=1,NSICE I = SICE_INDEX(K,1) IF(mskant(i) == 1_i8)THEN D_T = TSTAR_SICE(I) - TL_1(I) IF (D_T > 0.05_r8 .OR. D_T < -0.05_r8) THEN ALPHA1(I) = (QSTAR_ICE(I) - QS1(I)) / D_T ELSEIF (TL_1(I) > TM) THEN ALPHA1(I) = EPSILON*LC*QS1(I)*( 1.0_r8+C_VIRTUAL*QS1(I) ) & & /(R*TL_1(I)*TL_1(I)) ELSE ALPHA1(I) = EPSILON*LS*QS1(I)*( 1.0_r8+C_VIRTUAL*QS1(I) ) & & /(R*TL_1(I)*TL_1(I)) ENDIF ASHTF(I) = 2 * KAPPAI / DE + 4*SBCON*TI(I)**3 END IF ENDDO DO I=1,nCols IF(mskant(i) == 1_i8)THEN IF (FLANDG(I) < 1.0_r8 ) THEN E_SEA(I) = - (1.0_r8 - ICE_FRACT(I)) * RHOKH_1(I)*(QW_1(I) - SeaSalinityFactor*QSTAR_SEA(I)) H_SEA(I) = - (1.0_r8 - ICE_FRACT(I)) *CP* RHOKH_1(I) * ( TL_1(I) - TSTAR_SEA(I) & & + GRCP*(Z1_TQ(I) + Z0M_SEA(I) - Z0H_SEA(I)) ) IF ( ICE_FRACT(I) > 0.0_r8 ) THEN ! Sea-ice RHOKPM(I) = RHOKH_1(I) / ( ASHTF(I) + RHOKH_1(I)*(LS*ALPHA1(I) + CP) ) RAD_REDUC = RADNET(I) - ICE_FRACT(I) * ASHTF(I) * & & ( TL_1(I) - TI(I) + & & GRCP*(Z1_TQ(I) + Z0M_ICE(I) - Z0H_ICE(I)) ) DQ1 = QS1(I) - QW_1(I) + & & GRCP*ALPHA1(I)* & & (Z1_TQ(I) + Z0M_ICE(I) - Z0H_ICE(I)) FQW_ICE(I) = RHOKPM(I) * ( ALPHA1(I)*RAD_REDUC + & & (CP*RHOKH_1(I) + & & ASHTF(I))*DQ1*ICE_FRACT(I) ) FTL_ICE(I) = RHOKPM(I) * ( RAD_REDUC - & & ICE_FRACT(I)*LS*RHOKH_1(I)*DQ1 ) ENDIF FTL_1(I) = (1.0_r8-FLANDG(I))*(FTL_ICE(I) + H_SEA(I) / CP) FQW_1(I) = (1.0_r8-FLANDG(I))*(FQW_ICE(I) + E_SEA(I)) ENDIF END IF ENDDO ! IF (LTIMER) THEN ! DEPENDS ON: timer ! CALL TIMER('SF_FLUX ',4) ! ENDIF RETURN END SUBROUTINE SF_FLUX_SEA ! ! *****************************COPYRIGHT******************************* ! (C) Crown copyright Met Office. All rights reserved. ! For further details please refer to the file COPYRIGHT.txt ! which you should have received as part of this distribution. ! *****************************COPYRIGHT******************************* ! subroutine DUSTRESB ! ! ! Purpose: ! To calculate the surface layer resistance for mineral dust ! ! Called by sfexch ! ! Current owners of code: S.Woodward ! ! History: ! Version Date Comment ! ------- ---- ------- ! ! 5.5 12/02/03 Original code S Woodward ! 6.2 11/05/06 Fix for bit-non-reproducibility A. Malcolm ! 6.4 15/01/07 Malcolm McVean's portability fix S.Woodward ! ! Code Description: ! Language: FORTRAN77 + common extensions ! This code is written to UMDP3 v6 programming standards ! ! Documentation: "Modelling the atmospheric lifecycle..." ! Woodward, JGR106, D16, pp18155-18166 !--------------------------------------------------------------------- ! SUBROUTINE DUSTRESB( & nCols , &!INTEGER,INTENT(IN ) :: nCols !IN mskant , & PSTAR , &!REAL(KIND=r8), INTENT(IN ) :: PSTAR(nCols,ROWS) !IN surface pressure TSTAR , &!REAL(KIND=r8), INTENT(IN ) :: TSTAR(nCols,ROWS) !IN surface temperature RHOSTAR , &!REAL(KIND=r8), INTENT(IN ) :: RHOSTAR(nCols,ROWS)!IN surface air density VSHR , &!REAL(KIND=r8), INTENT(IN ) :: VSHR(nCols,ROWS) !IN surface to lowest lev windspeed difference CD_STD_DUST, &!REAL(KIND=r8), INTENT(IN ) :: CD_STD_DUST(nCols,ROWS) !IN surface transfer coeffient for R_B_DUST &! REAL(KIND=r8), INTENT(OUT ) :: R_B_DUST(nCols,ROWS,NDIV) !OUT surface layer resistance for ) IMPLICIT NONE !C_DUST_NDIV............................................................. ! Description: Contains parameters for mineral dust code ! Current Code Owner: Stephanie Woodward ! ! History: ! Version Date Comment ! ------- ---- ------- ! 5.5 12/02/03 Original Code. Stephanie Woodward ! ! Declarations: ! INTEGER NDIV ! number of particle size divisions PARAMETER (NDIV = 6) !..................................................................... INTEGER,INTENT(IN ) :: nCols !IN INTEGER(KIND=i8),INTENT(IN ) :: mskant(nCols) REAL(KIND=r8), INTENT(IN ) :: PSTAR(nCols) !IN surface pressure REAL(KIND=r8), INTENT(IN ) :: TSTAR(nCols) !IN surface temperature REAL(KIND=r8), INTENT(IN ) :: RHOSTAR(nCols)!IN surface air density REAL(KIND=r8), INTENT(IN ) :: VSHR(nCols) !IN surface to lowest lev windspeed difference REAL(KIND=r8), INTENT(IN ) :: CD_STD_DUST(nCols) !IN surface transfer coeffient for ! ! momentum, excluding orographic ! ! form drag REAL(KIND=r8), INTENT(OUT ) :: R_B_DUST(nCols,NDIV) !OUT surface layer resistance for ! ! mineral dust ! local variables INTEGER & & IDIV & !loop counter, dust divisions &,I & !loop counter &,LEV1 !number of levels for vstokes calculation REAL(KIND=r8) :: ETAA(nCols) & !dynamic viscosity of air &,VSTOKES1(nCols) & !gravitational settling velocity, lev1 &,NSTOKES(nCols) & !stokes number = VstokesVshrVshr/nu g &,NSCHMIDT(nCols) & !schmidt number = nu/diffusivit &,CCF(nCols) & !Cunningham correction factor &,STOKES_EXP & !stokes term in R_B_DUST equation &,SMALLP !small +ve number, negligible compared to 1 !*L------------------COMDECK C_G---------------------------------------- ! G IS MEAN ACCEL DUE TO GRAVITY AT EARTH'S SURFACE REAL(KIND=r8), PARAMETER :: G = 9.80665_r8 !*---------------------------------------------------------------------- !*L------------------COMDECK C_PI--------------------------------------- !LL !LL 4.0 19/09/95 New value for PI. Old value incorrect !LL from 12th decimal place. D. Robinson !LL 5.1 7/03/00 Fixed/Free format P.Selwood !LL ! Pi REAL(KIND=r8), PARAMETER :: Pi = 3.14159265358979323846_r8 ! Conversion factor degrees to radians REAL(KIND=r8), PARAMETER :: Pi_Over_180 = Pi/180.0_r8 ! Conversion factor radians to degrees REAL(KIND=r8), PARAMETER :: Recip_Pi_Over_180 = 180.0_r8/Pi !*---------------------------------------------------------------------- !*L------------------COMDECK C_O_DG_C----------------------------------- ! ZERODEGC IS CONVERSION BETWEEN DEGREES CELSIUS AND KELVIN ! TFS IS TEMPERATURE AT WHICH SEA WATER FREEZES ! TM IS TEMPERATURE AT WHICH FRESH WATER FREEZES AND ICE MELTS REAL(KIND=r8), PARAMETER :: ZeroDegC = 273.15_r8 REAL(KIND=r8), PARAMETER :: TFS = 271.35_r8 REAL(KIND=r8), PARAMETER :: TM = 273.15_r8 !*---------------------------------------------------------------------- !*L------------------COMDECK C_R_CP------------------------------------- ! History: ! Version Date Comment. ! 5.0 07/05/99 Add variable P_zero for consistency with ! conversion to C-P 'C' dynamics grid. R. Rawlins ! 5.1 07/03/00 Fixed/Free format conversion P. Selwood ! R IS GAS CONSTANT FOR DRY AIR ! CP IS SPECIFIC HEAT OF DRY AIR AT CONSTANT PRESSURE ! PREF IS REFERENCE SURFACE PRESSURE REAL(KIND=r8), PARAMETER :: R = 287.05_r8 REAL(KIND=r8), PARAMETER :: CP = 1005.0_r8 REAL(KIND=r8), PARAMETER :: Kappa = R/CP REAL(KIND=r8), PARAMETER :: Pref = 100000.0_r8 ! Reference surface pressure = PREF REAL(KIND=r8), PARAMETER :: P_zero = Pref REAL(KIND=r8), PARAMETER :: sclht = 6.8E+03_r8 ! Scale Height H !*---------------------------------------------------------------------- !C_DUSTGEN.............................................................. ! Description: Contains parameters for mineral dust generation ! Current Code Owner: Stephanie Woodward ! ! History: ! Version Date Comment ! ------- ---- ------- ! 5.5 12/02/03 Original Code. Stephanie Woodward ! 6.2 14/02/06 Alternative emissions terms for HadGEM1A. ! Stephanie Woodward ! 6.4 03/01/07 Alternative U*t terms. Stephanie Woodward ! Parameters for mineral dust generation ! REAL(KIND=r8), PARAMETER :: HORIZ_C = 2.61_r8 ! C in horizontal flux calc. REAL(KIND=r8), PARAMETER :: VERT_A = 13.4_r8 ! A in vertical flux calc REAL(KIND=r8), PARAMETER :: VERT_B = -6.0_r8 ! B in vertical flux calc REAL(KIND=r8),PARAMETER :: RHOP = 2.65E+3_r8 ! density of a dust particle REAL(KIND=r8),PARAMETER :: Z0B = 0.0003_r8 !roughness length for bare soil REAL(KIND=r8) Z0S(NDIV) ! smooth roughness len (calc.d from part size) REAL(KIND=r8) DREP(NDIV) ! representative particle diameter REAL(KIND=r8) DMAX(NDIV) ! max diameter of particles in each div. REAL(KIND=r8) DMIN(NDIV) ! min diameter of particles in each div. ! note that by using two arrays here we can set ! up overlapping divisions, however this means ! we have to be careful to make them consistent ! DATA Z0S/ .374894E-08_r8, .118552E-07_r8, .374894E-07_r8, .118552E-06_r8, & & .374894E-06_r8, .118552E-05_r8/ DATA DREP/ .112468E-06_r8, .355656E-06_r8, .112468E-05_r8, .355656E-05_r8, & & .112468E-04_r8, .355656E-04_r8/ DATA DMAX/2.0E-7_r8,6.32456E-7_r8,2.0E-6_r8, & & 6.32456E-6_r8,2.0E-5_r8,6.32456E-5_r8/ DATA DMIN/6.32456E-8_r8,2.0E-7_r8,6.32456E-7_r8, & & 2.0E-6_r8,6.32456E-6_r8,2.0E-5_r8/ !....................................................................... !C_DUSTGRAV............................................................. ! Description: ! Contains parameters for mineral dust gravitational settling ! Current Code Owner: Stephanie Woodward ! ! History: ! Version Date Comment ! ------- ---- ------- ! 1 12/02/03 Original Code. Stephanie Woodward ! REAL(KIND=r8),PARAMETER :: ACCF=1.257_r8 ! Cunningham correction factor term A REAL(KIND=r8),PARAMETER :: BCCF=0.4_r8 ! Cunningham correction factor term B REAL(KIND=r8),PARAMETER :: CCCF=-1.1_r8 ! Cunningham correction factor term C !....................................................................... !-------------------COMDECK C_SULCHM-------------------------------- ! Parameters for Sulphur Cycle Chemistry REAL(KIND=r8) & & EVAPTAU, & ! timescale for dissolved SO4 to evaporate & NUCTAU, & ! timescale for accumulation mode particles ! to nucleate once they enter a cloud. & DIFFUSE_AIT, & ! diffusion coefficient of Aitken particles & K_SO2OH_HI, & ! high pressure reaction rate limit & K_DMS_OH, & ! reaction rate for DMS+OH cc/mcl/s & K4_CH3SO2_O3, & ! Rate coeff for CH3SO2+O3 -> CH3SO3+O2 & K5_CH3SO3_HO2, & ! Rate coeff for CH3SO3+HO2 -> MSA+O2 & RMM_O3, & ! relative molecular mass O3 & BRAT_SO2, & ! branching ratio for SO2 in DMS oxidn & BRAT_MSA, & ! branching ratio for MSA in DMS oxidn & AVOGADRO, & ! no. of molecules in 1 mole & RMM_S, & ! relative molecular mass S kg/mole & RMM_H2O2, & ! relative molecular mass H2O2 kg/mole & RMM_HO2, & ! relative molecular mass HO2 kg/mole & RMM_AIR, & ! relative molecular mass dry air & RMM_W, & ! relative molecular mass water & RELM_S_H2O2, & ! rel atomic mass sulphur/RMM_H2O2 & RELM_S_2N, & ! rel atomic mass Sulphur/2*Nitrogen & PARH, & ! power of temp dependence of K_SO2OH_LO & K1, & ! parameters for calcn of K_SO2OH_LO & T1, & ! & FC, & ! parameters for interpolation between & FAC1, & ! LO and HI reaction rate limits & K2,K3,K4, & ! parameters for calcn of K_HO2_HO2 & T2,T3,T4, & ! & CLOUDTAU, & ! air parcel lifetime in cloud & CHEMTAU, & ! chem lifetime in cloud before oxidn & O3_MIN, & ! min mmr of O3 required for oxidn & THOLD ! threshold for cloud liquid water ! ! PARAMETER ( & & EVAPTAU = 300.0_r8, & ! secs (=5 mins) & NUCTAU = 30.0_r8, & ! secs & DIFFUSE_AIT = 1.7134E-9_r8, & ! sq m/s & K_SO2OH_HI = 2.0E-12_r8, & ! cc/mcl/s from STOCHEM model & K_DMS_OH = 9.1E-12_r8, & ! cc/mcl/s & K4_CH3SO2_O3 = 1.0E-14_r8, & ! cc/mcl/s & K5_CH3SO3_HO2 = 4.0E-11_r8, & & RMM_O3 = 4.8E-2_r8, & ! kg/mole & BRAT_SO2 = 0.9_r8, & & BRAT_MSA = 1.0_r8-BRAT_SO2, & & AVOGADRO = 6.022E23_r8, & ! per mole & RMM_S = 3.20E-2_r8, & ! kg/mole & RMM_H2O2 = 3.40E-2_r8, & ! kg/mole & RMM_HO2 = 3.30E-2_r8, & ! kg/mole & RMM_AIR = 2.896E-2_r8, & ! kg/mole & RMM_W = 1.8E-2_r8, & ! kg/mole & RELM_S_H2O2 = 3.206_r8/3.40_r8, & & RELM_S_2N = 3.206_r8/2.80_r8, & & PARH = 3.3_r8, & & K1 = 4.0E-31_r8, & ! (cc/mcl)2/s from STOCHEM & T1 = 300.0_r8, & ! K & FC = 0.45_r8, & ! from STOCHEM model & FAC1 = 1.1904_r8, & ! 0.75-1.27*LOG10(FC) from STOCHEM & K2 = 2.2E-13_r8, & ! cc/mcl/s & K3 = 1.9E-33_r8, & ! (cc/mcl)2/s & K4 = 1.4E-21_r8, & ! cc/mcl & T2 = 600.0_r8, & ! K & T3 = 890.0_r8, & ! K & T4 = 2200.0_r8, & ! K & CLOUDTAU = 1.08E4_r8, & ! secs (=3 hours) & CHEMTAU = 9.0E2_r8, & ! secs (=15 mins) & O3_MIN = 1.6E-8_r8, & !(kg/kg, equiv. 10ppbv) & THOLD = 1.0E-8_r8 & ! kg/kg & ) ! REAL(KIND=r8) RAD_AIT, & ! median radius of Aitken mode particles & DIAM_AIT, & ! " diameter " & RAD_ACC, & ! median radius of acccumulation mode & DIAM_ACC, & ! " diameter " & CHI, & ! mole fraction of S in particle & RHO_SO4, & ! density of SO4 particle & SIGMA, & ! standard devn of particle size distn ! for accumulation mode & E_PARM, & ! param relating size distns of Ait & Acc & NUM_STAR ! threshold concn of accu mode particles ! below which PSI=1 ! PARAMETER ( & & RAD_AIT = 6.5E-9_r8, & ! m & DIAM_AIT = 2.0_r8*RAD_AIT, & & RAD_ACC = 95.0E-9_r8, & ! m & DIAM_ACC = 2.0_r8*RAD_ACC, & & CHI = 32.0_r8/132.0_r8, & & RHO_SO4 = 1769.0_r8, & ! kg/m3 & SIGMA = 1.4_r8, & & E_PARM = 0.9398_r8, & & NUM_STAR = 1.0E6_r8 & ! m-3 & ) ! REAL(KIND=r8) BOLTZMANN !Boltzmanns constant. REAL(KIND=r8) MFP_REF !Reference value of mean free path. REAL(KIND=r8) TREF_MFP !Reference temperature for mean free path. REAL(KIND=r8) PREF_MFP !Reference pressure for mean free path. REAL(KIND=r8) SIGMA_AIT !Geometric standard deviation of the Aitken ! mode distribution. ! PARAMETER (BOLTZMANN = 1.3804E-23_r8) ! J K-1 PARAMETER (MFP_REF = 6.6E-8_r8 & ! m & , TREF_MFP = 293.15_r8 & ! K & , PREF_MFP = 1.01325E5_r8) ! Pa PARAMETER (SIGMA_AIT = 1.30_r8) ! !*--------------------------------------------------------------------- ! ! EXTERNAL VGRAV ! !... epsilon() is defined as almost negligible, so eps/100 is negligible ! SMALLP = EPSILON(1.0_r8) / 100.0_r8 ! !...calc stokes number, schmidt number and finally resistance ! LEV1=1 DO IDIV=1,NDIV ! DEPENDS ON: vgrav CALL VGRAV( & & nCols,mskant,LEV1,DREP(IDIV),RHOP,PSTAR,TSTAR, & & VSTOKES1,CCF,ETAA & & ) !CDIR NOVECTOR DO I= 1,nCols IF(mskant(i) == 1_i8)THEN NSCHMIDT(I)=3.0_r8*PI*ETAA(I)*ETAA(I)*DREP(IDIV)/ & & (RHOSTAR(I)*BOLTZMANN*TSTAR(I)*CCF(I)) NSTOKES(I)=VSTOKES1(I)*CD_STD_DUST(I)*RHOSTAR(I)* & & VSHR(I)*VSHR(I)/(ETAA(I)*G) ! Avoid underflow in Stokes term by setting to zero if ! negligible compared to Schmidt term, i.e., if NSTOKES ! is too small. IF ( 3.0_r8 / NSTOKES(I) < & - LOG10( SMALLP *NSCHMIDT(I)**(-2.0_r8/3.0_r8) ) ) THEN STOKES_EXP = 10.0_r8**(-3.0_r8/NSTOKES(I)) ELSE STOKES_EXP = 0.0_r8 ENDIF R_B_DUST(I,IDIV)=1.0_r8/( SQRT(CD_STD_DUST(I)) * & & (NSCHMIDT(I)**(-2.0_r8/3.0_r8)+STOKES_EXP) ) END IF ENDDO !nCols ENDDO !NDIV RETURN END SUBROUTINE DUSTRESB ! *****************************COPYRIGHT******************************* ! (C) Crown copyright Met Office. All rights reserved. ! For further details please refer to the file COPYRIGHT.txt ! which you should have received as part of this distribution. ! *****************************COPYRIGHT******************************* ! Subroutine VGRAV ------------------------------------------------- ! ! Purpose: To calculate the gravitational sedimentation velocity of ! tracer particles according to Stoke's law, including the ! Cunningham correction factor. ! ! Current owners of code: S Woodward, M Woodage ! ! History: ! Version Date Comment ! ------- ---- ------- ! 4.4 03/10/97 Original code S Woodward, M Woodage ! 5.5 12/02/03 Updated for vn 5.5 S Woodward ! ! ! Code description: ! Language: FORTRAN77 + extensions ! Programming standard: UMDP 3 Vn 6 ! ! System components covered: ! ! System task: ! !Documentation: Ref. Pruppacher & Klett ! Microphysics of clouds & ppn 1978,1980 edns. ! !----------------------------------------------------------------------- ! SUBROUTINE VGRAV( & & nCols,mskant,NLEVS,DIAM,RHOP,P,T, & & VSTOKES,CCF,ETAA) ! IMPLICIT NONE ! ! INTEGER nCols !IN row length INTEGER NLEVS !IN number of levels ! INTEGER(KIND=i8) :: mskant(nCols) REAL(KIND=r8) DIAM !IN particle diameter REAL(KIND=r8) RHOP !IN particles density REAL(KIND=r8) P(nCols,NLEVS)!IN pressure REAL(KIND=r8) T(nCols,NLEVS)!IN temperature ! REAL(KIND=r8) VSTOKES(nCols,NLEVS) !OUT sedimentation velocity REAL(KIND=r8) ETAA(nCols,NLEVS)!OUT viscosity of air REAL(KIND=r8) CCF(nCols,NLEVS) !OUT cunningham correction factor ! ! !*L------------------COMDECK C_G---------------------------------------- ! G IS MEAN ACCEL DUE TO GRAVITY AT EARTH'S SURFACE REAL(KIND=r8), PARAMETER :: G = 9.80665_r8 !*---------------------------------------------------------------------- !*L------------------COMDECK C_O_DG_C----------------------------------- ! ZERODEGC IS CONVERSION BETWEEN DEGREES CELSIUS AND KELVIN ! TFS IS TEMPERATURE AT WHICH SEA WATER FREEZES ! TM IS TEMPERATURE AT WHICH FRESH WATER FREEZES AND ICE MELTS REAL(KIND=r8), PARAMETER :: ZeroDegC = 273.15_r8 REAL(KIND=r8), PARAMETER :: TFS = 271.35_r8 REAL(KIND=r8), PARAMETER :: TM = 273.15_r8 !*---------------------------------------------------------------------- !-------------------COMDECK C_SULCHM-------------------------------- ! Parameters for Sulphur Cycle Chemistry REAL(KIND=r8) & & EVAPTAU, & ! timescale for dissolved SO4 to evaporate & NUCTAU, & ! timescale for accumulation mode particles ! to nucleate once they enter a cloud. & DIFFUSE_AIT, & ! diffusion coefficient of Aitken particles & K_SO2OH_HI, & ! high pressure reaction rate limit & K_DMS_OH, & ! reaction rate for DMS+OH cc/mcl/s & K4_CH3SO2_O3, & ! Rate coeff for CH3SO2+O3 -> CH3SO3+O2 & K5_CH3SO3_HO2, & ! Rate coeff for CH3SO3+HO2 -> MSA+O2 & RMM_O3, & ! relative molecular mass O3 & BRAT_SO2, & ! branching ratio for SO2 in DMS oxidn & BRAT_MSA, & ! branching ratio for MSA in DMS oxidn & AVOGADRO, & ! no. of molecules in 1 mole & RMM_S, & ! relative molecular mass S kg/mole & RMM_H2O2, & ! relative molecular mass H2O2 kg/mole & RMM_HO2, & ! relative molecular mass HO2 kg/mole & RMM_AIR, & ! relative molecular mass dry air & RMM_W, & ! relative molecular mass water & RELM_S_H2O2, & ! rel atomic mass sulphur/RMM_H2O2 & RELM_S_2N, & ! rel atomic mass Sulphur/2*Nitrogen & PARH, & ! power of temp dependence of K_SO2OH_LO & K1, & ! parameters for calcn of K_SO2OH_LO & T1, & ! & FC, & ! parameters for interpolation between & FAC1, & ! LO and HI reaction rate limits & K2,K3,K4, & ! parameters for calcn of K_HO2_HO2 & T2,T3,T4, & ! & CLOUDTAU, & ! air parcel lifetime in cloud & CHEMTAU, & ! chem lifetime in cloud before oxidn & O3_MIN, & ! min mmr of O3 required for oxidn & THOLD ! threshold for cloud liquid water ! ! PARAMETER ( & & EVAPTAU = 300.0_r8, & ! secs (=5 mins) & NUCTAU = 30.0_r8, & ! secs & DIFFUSE_AIT = 1.7134E-9_r8, & ! sq m/s & K_SO2OH_HI = 2.0E-12_r8, & ! cc/mcl/s from STOCHEM model & K_DMS_OH = 9.1E-12_r8, & ! cc/mcl/s & K4_CH3SO2_O3 = 1.0E-14_r8, & ! cc/mcl/s & K5_CH3SO3_HO2 = 4.0E-11_r8, & & RMM_O3 = 4.8E-2_r8, & ! kg/mole & BRAT_SO2 = 0.9_r8, & & BRAT_MSA = 1.0_r8-BRAT_SO2, & & AVOGADRO = 6.022E23_r8, & ! per mole & RMM_S = 3.20E-2_r8, & ! kg/mole & RMM_H2O2 = 3.40E-2_r8, & ! kg/mole & RMM_HO2 = 3.30E-2_r8, & ! kg/mole & RMM_AIR = 2.896E-2_r8, & ! kg/mole & RMM_W = 1.8E-2_r8, & ! kg/mole & RELM_S_H2O2 = 3.206_r8/3.40_r8, & & RELM_S_2N = 3.206_r8/2.80_r8, & & PARH = 3.3_r8, & & K1 = 4.0E-31_r8, & ! (cc/mcl)2/s from STOCHEM & T1 = 300.0_r8, & ! K & FC = 0.45_r8, & ! from STOCHEM model & FAC1 = 1.1904_r8, & ! 0.75-1.27*LOG10(FC) from STOCHEM & K2 = 2.2E-13_r8, & ! cc/mcl/s & K3 = 1.9E-33_r8, & ! (cc/mcl)2/s & K4 = 1.4E-21_r8, & ! cc/mcl & T2 = 600.0_r8, & ! K & T3 = 890.0_r8, & ! K & T4 = 2200.0_r8, & ! K & CLOUDTAU = 1.08E4_r8, & ! secs (=3 hours) & CHEMTAU = 9.0E2_r8, & ! secs (=15 mins) & O3_MIN = 1.6E-8_r8, & !(kg/kg, equiv. 10ppbv) & THOLD = 1.0E-8_r8 & ! kg/kg & ) ! REAL(KIND=r8) RAD_AIT, & ! median radius of Aitken mode particles & DIAM_AIT, & ! " diameter " & RAD_ACC, & ! median radius of acccumulation mode & DIAM_ACC, & ! " diameter " & CHI, & ! mole fraction of S in particle & RHO_SO4, & ! density of SO4 particle & SIGMA, & ! standard devn of particle size distn ! for accumulation mode & E_PARM, & ! param relating size distns of Ait & Acc & NUM_STAR ! threshold concn of accu mode particles ! below which PSI=1 ! PARAMETER ( & & RAD_AIT = 6.5E-9_r8, & ! m & DIAM_AIT = 2.0_r8*RAD_AIT, & & RAD_ACC = 95.0E-9_r8, & ! m & DIAM_ACC = 2.0_r8*RAD_ACC, & & CHI = 32.0_r8/132.0_r8, & & RHO_SO4 = 1769.0_r8, & ! kg/m3 & SIGMA = 1.4_r8, & & E_PARM = 0.9398_r8, & & NUM_STAR = 1.0E6_r8 & ! m-3 & ) ! REAL(KIND=r8) BOLTZMANN !Boltzmanns constant. REAL(KIND=r8) MFP_REF !Reference value of mean free path. REAL(KIND=r8) TREF_MFP !Reference temperature for mean free path. REAL(KIND=r8) PREF_MFP !Reference pressure for mean free path. REAL(KIND=r8) SIGMA_AIT !Geometric standard deviation of the Aitken ! mode distribution. ! PARAMETER (BOLTZMANN = 1.3804E-23_r8) ! J K-1 PARAMETER (MFP_REF = 6.6E-8_r8 & ! m & , TREF_MFP = 293.15_r8 & ! K & , PREF_MFP = 1.01325E5_r8) ! Pa PARAMETER (SIGMA_AIT = 1.30_r8) ! !*--------------------------------------------------------------------- !C_DUSTGRAV............................................................. ! Description: ! Contains parameters for mineral dust gravitational settling ! Current Code Owner: Stephanie Woodward ! ! History: ! Version Date Comment ! ------- ---- ------- ! 1 12/02/03 Original Code. Stephanie Woodward ! REAL(KIND=r8),PARAMETER :: ACCF=1.257_r8 ! Cunningham correction factor term A REAL(KIND=r8),PARAMETER :: BCCF=0.4_r8 ! Cunningham correction factor term B REAL(KIND=r8),PARAMETER :: CCCF=-1.1_r8 ! Cunningham correction factor term C !....................................................................... ! ! local variables ! INTEGER ILEV !LOC loop counter for levels INTEGER I !LOC loop counter ! REAL(KIND=r8) TC(nCols) !LOC temperature in deg C REAL(KIND=r8) LAMDAA(nCols)!LOC mean free path of particle REAL(KIND=r8) ALPHACCF(nCols)!LOC ! ! Calculate viscosity of air (Pruppacher & Klett p.323) DO ILEV=1,NLEVS DO I = 1,nCols IF(mskant(i) == 1_i8)THEN TC(I)=T(I,ILEV)-ZERODEGC IF (TC(I) >= 0.0_r8) THEN ETAA(I,ILEV)=(1.718_r8+0.0049_r8*TC(I))*1.E-5_r8 ELSE ETAA(I,ILEV)=(1.718_r8+0.0049_r8*TC(I)-1.2E-5_r8*TC(I)*TC(I))*1.E-5_r8 ENDIF ! END IF ENDDO !nCols ENDDO !NLEVS ! DO ILEV=1,NLEVS DO I=1,nCols IF(mskant(i) == 1_i8)THEN ! ! Calculate mean free path of particle (Pruppacher & Klett p.323) LAMDAA(I)=MFP_REF*PREF_MFP*T(I,ILEV)/ & & (P(I,ILEV)*TREF_MFP) ! Calculate Cunningham correction factor(Pruppacher & Klett p.361) ALPHACCF(I)=ACCF+BCCF*EXP(CCCF*DIAM*.5_r8/LAMDAA(I)) CCF(I,ILEV)=(1.0_r8+ALPHACCF(I)*LAMDAA(I)/(0.5_r8*DIAM)) ! Calculate sedimentation velocity (Pruppacher & Klett p.362) VSTOKES(I,ILEV)=CCF(I,ILEV)*(DIAM*DIAM*G*RHOP)/ & & (18.0_r8*ETAA(I,ILEV)) ! END IF ENDDO !nCols ENDDO !NLEV ! RETURN END SUBROUTINE VGRAV ! *****************************COPYRIGHT******************************* ! (C) Crown copyright Met Office. All rights reserved. ! For further details please refer to the file COPYRIGHT.txt ! which you should have received as part of this distribution. ! *****************************COPYRIGHT******************************* !!! SUBROUTINES SFL_INT_SEA AND SFL_INT_LAND-------------------------- !!! !!! Purpose: To calculate interpolation coefficients for 10m winds !!! and 1.5m temperature/specific humidity diagnostics. !!! !!! Model Modification history: !!! version Date !!! !!! 5.2 15/11/00 New Deck M. Best ! 5.3 25/04/01 Add coastal tiling. Nic Gedney !!! !!! Programming standard: !!! !!! Logical component covered: Part of P243. !!! !!! System Task: !!! !!! External Documentation: UMDP No.24 !!! !!!--------------------------------------------------------------------- !*L Arguments :- SUBROUTINE SFL_INT_SEA ( & & nCols,mskant,FLANDG & &,VSHR,CD,CH,Z0M,Z0H & &,RECIP_L_MO,V_S & &,SU10,SV10,ST1P5,SQ1P5,CDR10M,CHR1P5M) IMPLICIT NONE INTEGER :: nCols ! IN Number of X points? INTEGER(KIND=i8) :: mskant(nCols) REAL(KIND=r8) & & FLANDG(nCols) & ! ! IN Land fraction. &,Z0M(nCols) & ! IN Roughness length for momentum (m). &,Z0H(nCols) & ! IN Roughness length for heat and ! ! moisture (m). &,VSHR(nCols) & ! IN Wind speed difference between the ! ! surface and the lowest wind level in ! ! the atmosphere (m/s). &,CD(nCols) & ! IN Surface drag coefficient. &,CH(nCols) & ! IN Surface transfer coefficient for heat ! ! and moisture. &,RECIP_L_MO(nCols) & ! ! IN Reciprocal of the Monin-Obukhov ! ! length (m) &,V_S(nCols) ! IN Surface layer scaling velocity ! ! including orographic form drag (m/s). LOGICAL & & SU10 & ! IN 10m U-wind diagnostic flag &,SV10 & ! IN 10m V-wind diagnostic flag &,ST1P5 & ! IN screen temp diagnostic flag &,SQ1P5 ! IN screen specific humidity ! ! diagnostic flag ! Output variables ! REAL(KIND=r8) :: CDR10M(nCols) & ! ! OUT interpolation coefficicent for ! ! 10m wind &,CHR1P5M(nCols) ! OUT Interpolation coefficient for ! ! 1.5m temperature !* !*L--------------------------------------------------------------------- ! External routines called :- ! EXTERNAL PHI_M_H_SEA ! EXTERNAL TIMER !* !*L--------------------------------------------------------------------- ! Local and other symbolic constants :- ! C_VKMAN start REAL(KIND=r8),PARAMETER:: VKMAN=0.4_r8 ! Von Karman's constant ! C_VKMAN end REAL(KIND=r8) Z_OBS_TQ,Z_OBS_WIND PARAMETER ( & & Z_OBS_TQ = 1.5_r8 & ! Height of screen observations of temperature ! ! and humidity. &,Z_OBS_WIND = 10.0_r8 & ! Height of surface wind observations. &) ! ! Define local storage. ! ! (a) Local work arrays. ! REAL(KIND=r8) & & Z_WIND(nCols) & ! Height of wind observations. &,Z_TEMP(nCols) & ! Height of temperature and humidity ! ! observations. &,PHI_M_OBS(nCols) & ! Monin-Obukhov stability function for ! ! momentum integrated to the wind ! ! observation height. &,PHI_H_OBS(nCols)! Monin-Obukhov stability function for ! ! scalars integrated to their ! ! observation height. ! ! (b) Scalars. ! INTEGER & & I ! Loop counter (horizontal field index). !* ! IF (LTIMER) THEN ! DEPENDS ON: timer ! CALL TIMER('SFL_INT ',3) ! ENDIF ! !----------------------------------------------------------------------- !! 1. If diagnostics required calculate M-O stability functions at !! observation heights. !----------------------------------------------------------------------- IF (SU10 .OR. SV10 .OR. ST1P5 .OR. SQ1P5) THEN DO I=1,nCols IF ( FLANDG(I) < 1.0_r8 ) THEN Z_WIND(I) = Z_OBS_WIND Z_TEMP(I) = Z_OBS_TQ + Z0H(I) - Z0M(I) END IF ENDDO ! DEPENDS ON: phi_m_h_sea CALL PHI_M_H_SEA (nCols,mskant,FLANDG, & & RECIP_L_MO,Z_WIND,Z_TEMP,Z0M,Z0H, & & PHI_M_OBS,PHI_H_OBS) ENDIF !----------------------------------------------------------------------- !! 2. If diagnostics required calculate interpolation coefficient !! for 1.5m screen temperature and specific humidity. !----------------------------------------------------------------------- ! IF (ST1P5 .OR. SQ1P5) THEN DO I=1,nCols IF ( FLANDG(I) < 1.0_r8 ) THEN CHR1P5M(I) = CH(I) * VSHR(I) * & & PHI_H_OBS(I)/(VKMAN*V_S(I)) ENDIF ENDDO ENDIF ! !----------------------------------------------------------------------- !! 3. If diagnostics required calculate interpolation coefficient !! for 10m winds. !----------------------------------------------------------------------- ! IF ( SU10 .OR. SV10 ) THEN DO I=1,nCols IF ( FLANDG(I) < 1.0_r8 ) THEN CDR10M(I) = (1.0_r8-FLANDG(I)) * CD(I) * VSHR(I) * & & PHI_M_OBS(I)/(VKMAN*V_S(I)) ENDIF ENDDO ENDIF ! ! IF (LTIMER) THEN ! DEPENDS ON: timer ! CALL TIMER('SFL_INT ',4) ! ENDIF RETURN END SUBROUTINE SFL_INT_SEA !!! !!!--------------------------------------------------------------------- !*L Arguments :- ! *****************************COPYRIGHT******************************* ! (C) Crown copyright Met Office. All rights reserved. ! For further details please refer to the file COPYRIGHT.txt ! which you should have received as part of this distribution. ! *****************************COPYRIGHT******************************* ! !!! SUBROUTINES STDEV1_SEA and STDEV1_LAND ---------------------------- !!! !!! Purpose: Calculate the standard deviations of layer 1 turbulent !!! fluctuations of temperature and humidity using approximate !!! formulae from first order closure. !!! !!! Model Modification history !!! version date !!! 5.2 15/11/00 New Deck M. Best ! 5.3 25/04/01 Add coastal tiling. Nic Gedney !!! !!! Programming standard: !!! !!! ------------------------------------------------------------------- ! !!! SUBROUTINE STDEV1_SEA --------------------------------------------- !!! Layer 1 standard deviations for sea and sea-ice !!! ------------------------------------------------------------------- SUBROUTINE STDEV1_SEA (& nCols , &!INTEGER, INTENT(IN ) :: nCols! IN Number of X points? mskant , & FLANDG , &!REAL(KIND=r8) , INTENT(IN ) :: FLANDG(nCols,ROWS)! IN Land fraction. BQ_1 , &!REAL(KIND=r8) , INTENT(IN ) :: BQ_1(nCols,ROWS) ! IN Buoyancy parameter. BT_1 , &!REAL(KIND=r8) , INTENT(IN ) :: BT_1(nCols,ROWS) ! IN Buoyancy parameter. FQW_1 , &!REAL(KIND=r8) , INTENT(IN ) :: FQW_1(nCols,ROWS) ! IN Surface flux of QW. FTL_1 , &!REAL(KIND=r8) , INTENT(IN ) :: FTL_1(nCols,ROWS) ! IN Surface flux of TL. ICE_FRACT , &!REAL(KIND=r8) , INTENT(IN ) :: ICE_FRACT(nCols,ROWS) ! IN Fraction of gridbox which is sea-ice. RHOKM_1 , &!REAL(KIND=r8) , INTENT(IN ) :: RHOKM_1(nCols,ROWS) ! IN Surface momentum exchange RHOSTAR , &!REAL(KIND=r8) , INTENT(IN ) :: RHOSTAR(nCols,ROWS) ! IN Surface air density. VSHR , &!REAL(KIND=r8) , INTENT(IN ) :: VSHR(nCols,ROWS) ! IN Magnitude of surface- Z0MSEA , &!REAL(KIND=r8) , INTENT(IN ) :: Z0MSEA(nCols,ROWS) ! IN Sea roughness length. Z0_ICE , &!REAL(KIND=r8) , INTENT(IN ) :: Z0_ICE(nCols,ROWS) ! IN Sea-ice roughness length. Z1_TQ , &!REAL(KIND=r8) , INTENT(IN ) :: Z1_TQ(nCols,ROWS) ! IN Height of lowest tq level. Q1_SD , &!REAL(KIND=r8) , INTENT(out ) :: Q1_SD(nCols,ROWS) ! OUT Standard deviation of T1_SD &!REAL(KIND=r8) , INTENT(out ) :: T1_SD(nCols,ROWS) ! OUT Standard deviation of ) IMPLICIT NONE INTEGER, INTENT(IN ) :: nCols! IN Number of X points? INTEGER(KIND=i8), INTENT(IN ) :: mskant(nCols) REAL(KIND=r8) , INTENT(IN ) :: FLANDG(nCols)! IN Land fraction. REAL(KIND=r8) , INTENT(IN ) :: BQ_1(nCols) ! IN Buoyancy parameter. REAL(KIND=r8) , INTENT(IN ) :: BT_1(nCols) ! IN Buoyancy parameter. REAL(KIND=r8) , INTENT(IN ) :: FQW_1(nCols) ! IN Surface flux of QW. REAL(KIND=r8) , INTENT(IN ) :: FTL_1(nCols) ! IN Surface flux of TL. REAL(KIND=r8) , INTENT(IN ) :: ICE_FRACT(nCols) ! IN Fraction of gridbox which is sea-ice. REAL(KIND=r8) , INTENT(IN ) :: RHOKM_1(nCols) ! IN Surface momentum exchange ! ! coefficient. REAL(KIND=r8) , INTENT(IN ) :: RHOSTAR(nCols) ! IN Surface air density. REAL(KIND=r8) , INTENT(IN ) :: VSHR(nCols) ! IN Magnitude of surface- ! ! to-lowest-level wind shear. REAL(KIND=r8) , INTENT(IN ) :: Z0MSEA(nCols) ! IN Sea roughness length. REAL(KIND=r8) , INTENT(IN ) :: Z0_ICE(nCols) ! IN Sea-ice roughness length. REAL(KIND=r8) , INTENT(IN ) :: Z1_TQ(nCols) ! IN Height of lowest tq level. REAL(KIND=r8) , INTENT(out ) :: Q1_SD(nCols) ! OUT Standard deviation of ! ! turbulent fluctuations of ! ! surface layer specific ! ! humidity (kg/kg). REAL(KIND=r8) , INTENT(out ) :: T1_SD(nCols) ! OUT Standard deviation of ! ! turbulent fluctuations of ! ! surface layer temperature (K). ! External routines called :- ! EXTERNAL TIMER !*L------------------COMDECK C_G---------------------------------------- ! G IS MEAN ACCEL DUE TO GRAVITY AT EARTH'S SURFACE REAL(KIND=r8), PARAMETER :: G = 9.80665_r8 !*---------------------------------------------------------------------- ! Workspace -------------------------------------------------------- INTEGER & & I ! Loop counter (horizontal field index). REAL(KIND=r8) :: VS & ! Surface layer friction velocity &,VSF1_CUBED & ! Cube of surface layer free convective ! ! scaling velocity &,WS1 & ! Turbulent velocity scale for surface ! ! layer &,Z0 ! Roughness length ! IF (LTIMER) THEN ! DEPENDS ON: timer ! CALL TIMER('STDEV1 ',3) ! ENDIF DO I=1,nCols IF(mskant(i) == 1_i8)THEN IF ( FLANDG(I) < 1.0_r8 ) THEN Z0 = Z0MSEA(I) IF ( ICE_FRACT(I) > 0.0_r8 ) Z0 = Z0_ICE(I) VS = SQRT ( RHOKM_1(I)/RHOSTAR(I) * VSHR(I) ) VSF1_CUBED = 1.25_r8*G*(Z1_TQ(I) + Z0) * & & ( BT_1(I)*FTL_1(I) + & & BQ_1(I)*FQW_1(I) )/RHOSTAR(I) IF ( VSF1_CUBED > 0.0_r8 ) THEN WS1 = ( VSF1_CUBED + VS*VS*VS ) ** (1.0_r8/3.0_r8) T1_SD(I) = MAX ( 0.0_r8 , & & (1.0_r8-FLANDG(I))*1.93_r8*FTL_1(I) / (RHOSTAR(I)*WS1) ) Q1_SD(I) = MAX ( 0.0_r8 , & & (1.0_r8-FLANDG(I))*1.93_r8*FQW_1(I) / (RHOSTAR(I)*WS1) ) ENDIF ENDIF END IF ENDDO ! IF (LTIMER) THEN ! DEPENDS ON: timer ! CALL TIMER('STDEV1 ',4) ! ENDIF RETURN END SUBROUTINE STDEV1_SEA !!! SUBROUTINE STDEV1_LAND -------------------------------------------- !!! Layer 1 standard deviations for land tiles !!! ------------------------------------------------------------------- ! *****************************COPYRIGHT******************************* ! (C) Crown copyright Met Office. All rights reserved. ! For further details please refer to the file COPYRIGHT.txt ! which you should have received as part of this distribution. ! *****************************COPYRIGHT******************************* !!! SUBROUTINES FCDCH_SEA AND FCDCH_LAND----------------------------- !!! !!! Purpose: Calculate surface transfer coefficients at one or more !!! gridpoints. !!! !!! Model Modification history: !!! version Date !!! !!! 5.2 15/11/00 New Deck M. Best ! 5.3 25/04/01 Add coastal tiling. Nic Gedney ! 5.5 12/02/03 Include code for mineral dust scheme. S Woodward ! 6.1 08/12/03 Add !CDIR NODEP to force vectorisation. R Barnes !!! !!! Programming standard: !!! !!! System component covered: Part of P243. !!! !!! Documentation: UM Documentation Paper No 24, section P243. !!! ! SUBROUTINE FCDCH_SEA--------------------------------------------- ! ! Transfer coefficients for sea, sea-ice and leads ! ! ----------------------------------------------------------------- !! Arguments:--------------------------------------------------------- SUBROUTINE FCDCH_SEA(& nCols , & !INTEGER, INTENT(IN ) :: nCols ! IN Number of X points? mskant , & FLANDG , & !REAL(KIND=r8) , INTENT(IN ) :: FLANDG(nCols,ROWS) ! IN Land fraction DB , & !REAL(KIND=r8) , INTENT(IN ) :: DB (nCols,ROWS) ! VSHR , & !REAL(KIND=r8) , INTENT(IN ) :: VSHR (nCols,ROWS) Z0M , & !REAL(KIND=r8) , INTENT(IN ) :: Z0M (nCols,ROWS) Z0H , & !REAL(KIND=r8) , INTENT(IN ) :: Z0H (nCols,ROWS) ZH , & !REAL(KIND=r8) , INTENT(IN ) :: ZH (nCols,ROWS) Z1_UV , & !REAL(KIND=r8) , INTENT(IN ) :: Z1_UV(nCols,ROWS) Z1_TQ , & !REAL(KIND=r8) , INTENT(IN ) :: Z1_TQ(nCols,ROWS) CDV , & !REAL(KIND=r8) , INTENT(OUT ) :: CDV (nCols,ROWS) CHV , & !REAL(KIND=r8) , INTENT(OUT ) :: CHV (nCols,ROWS) V_S , & !REAL(KIND=r8) , INTENT(OUT ) :: V_S(nCols,ROWS) RECIP_L_MO & !REAL(KIND=r8) , INTENT(OUT ) :: RECIP_L_MO(nCols,ROWS) ) IMPLICIT NONE INTEGER, INTENT(IN ) :: nCols ! IN Number of X points? INTEGER(KIND=i8), INTENT(IN ) :: mskant(nCols) REAL(KIND=r8) , INTENT(IN ) :: FLANDG(nCols) ! IN Land fraction REAL(KIND=r8) , INTENT(IN ) :: DB (nCols) ! IN Buoyancy difference between surface ! and lowest temperature and humidity ! level in the atmosphere (m/s^2). REAL(KIND=r8) , INTENT(IN ) :: VSHR (nCols) ! IN Wind speed difference between the ! surface and the lowest wind level in ! the atmosphere (m/s). REAL(KIND=r8) , INTENT(IN ) :: Z0M (nCols) ! IN Roughness length for momentum ! transport (m). REAL(KIND=r8) , INTENT(IN ) :: Z0H (nCols) ! IN Roughness length for heat and ! moisture (m). REAL(KIND=r8) , INTENT(IN ) :: ZH (nCols) ! IN Depth of boundary layer (m). REAL(KIND=r8) , INTENT(IN ) :: Z1_UV(nCols) ! IN Height of lowest wind level (m). REAL(KIND=r8) , INTENT(IN ) :: Z1_TQ(nCols) ! IN Height of lowest temperature and ! humidity level (m). REAL(KIND=r8) , INTENT(OUT ) :: CDV (nCols) ! OUT Surface transfer coefficient for ! momentum including orographic form drag (m/s). REAL(KIND=r8) , INTENT(OUT ) :: CHV (nCols) ! OUT Surface transfer coefficient for ! heat, moisture & other scalars (m/s). REAL(KIND=r8) , INTENT(OUT ) :: V_S(nCols) ! OUT Surface layer scaling velocity ! ! including orographic form drag (m/s). REAL(KIND=r8) , INTENT(OUT ) :: RECIP_L_MO(nCols) ! OUT Reciprocal of the Monin-Obukhov length ! (m^-1). !*L Workspace usage---------------------------------------------------- ! ! Local work arrays. ! REAL(KIND=r8) & & PHI_M(nCols) & ! Monin-Obukhov stability function for ! ! momentum integrated to the model's ! ! lowest wind level. &,PHI_H(nCols) ! Monin-Obukhov stability function for ! ! scalars integrated to the model's lowest ! ! temperature and humidity level. ! !*---------------------------------------------------------------------- ! EXTERNAL PHI_M_H_SEA ! EXTERNAL TIMER !*---------------------------------------------------------------------- ! Common and local constants. ! C_VKMAN start REAL(KIND=r8),PARAMETER:: VKMAN=0.4_r8 ! Von Karman's constant ! C_VKMAN end REAL(KIND=r8),PARAMETER:: CTRL_TUNABLE=1.0_r8 ! Tunable parameter in the surface layer scaling ! ! velocity formula (multiplying the turbulent ! ! convective scaling velocity). REAL(KIND=r8) BETA,THIRD PARAMETER ( & & BETA=0.009_r8, & ! Tunable parameter in the surface layer scaling ! ! velocity formula (multiplying the turbulent ! ! convective scaling velocity). & THIRD=1.0_r8/3.0_r8 & ! One third. &) INTEGER, PARAMETER :: N_ITS=8 ! 5 Number of iterations for Monin-Obukhov length ! ! and stability functions. ! ! Define local variables ! INTEGER I ! Loop counters; horizontal field index. INTEGER IT ! Iteration loop counter. REAL(KIND=r8) & & B_FLUX & ! Surface bouyancy flux over air density. &,U_S & ! Surface friction velocity (effective value). &,W_S ! Surface turbulent convective scaling velocity. ! IF (LTIMER) THEN ! DEPENDS ON: timer ! CALL TIMER('FCDCH ',3) ! ENDIF ! !----------------------------------------------------------------------- !! 1. Set initial values for the iteration. !----------------------------------------------------------------------- ! DO I=1,nCols IF(mskant(i) == 1_i8)THEN IF ( FLANDG(I) < 1.0_r8 ) THEN IF (DB(I) < 0.0_r8 .AND. VSHR(I) < 2.0_r8) THEN !----------------------------------------------------------------------- ! Start the iteration from the convective limit. !----------------------------------------------------------------------- RECIP_L_MO(I) = -VKMAN/(BETA*BETA*BETA*ZH(I)) ELSE !----------------------------------------------------------------------- ! Start the iteration from neutral values. !----------------------------------------------------------------------- RECIP_L_MO(I) = 0.0_r8 ENDIF ENDIF ! SEA_MASK END IF ENDDO ! DEPENDS ON: phi_m_h_sea CALL PHI_M_H_SEA (nCols,mskant,FLANDG, & & RECIP_L_MO,Z1_UV,Z1_TQ,Z0M,Z0H, & & PHI_M,PHI_H) ! DO I=1,nCols IF(mskant(i) == 1_i8)THEN IF ( FLANDG(I) < 1.0_r8 ) THEN IF (DB(I) < 0.0_r8 .AND. VSHR(I) < 2.0_r8) THEN !----------------------------------------------------------------------- ! Start the iteration from the convective limit. !----------------------------------------------------------------------- V_S(I) = BETA * & & SQRT( BETA * ( VKMAN / PHI_H(I) ) * & & ZH(I) * (-DB(I)) ) ELSE !----------------------------------------------------------------------- ! Start the iteration from neutral values. !----------------------------------------------------------------------- V_S(I) = ( VKMAN / PHI_M(I) ) * VSHR(I) ENDIF CHV(I) = ( VKMAN / PHI_H(I) ) * V_S(I) CDV(I) = ( VKMAN / PHI_M(I) ) * V_S(I) ENDIF ! SEA_MASK END IF ENDDO !----------------------------------------------------------------------- !! 2. Iterate to obtain sucessively better approximations for CD & CH. !----------------------------------------------------------------------- DO IT = 1,N_ITS ! DO I=1,nCols IF(mskant(i) == 1_i8)THEN IF ( FLANDG(I) < 1.0_r8 ) THEN B_FLUX = -CHV(I) * DB(I) U_S = SQRT( CDV(I) * VSHR(I) ) IF (DB(I) < 0.0_r8) THEN W_S = (ZH(I) * B_FLUX)**THIRD V_S(I) = SQRT(U_S*U_S + BETA*BETA*W_S*W_S) ELSE V_S(I) = U_S ENDIF RECIP_L_MO(I) = -VKMAN * B_FLUX / & & (V_S(I)*V_S(I)*V_S(I)) ENDIF ! SEA_MASK END IF ENDDO ! DEPENDS ON: phi_m_h_sea CALL PHI_M_H_SEA (nCols,mskant,FLANDG, & & RECIP_L_MO,Z1_UV,Z1_TQ,Z0M,Z0H, & & PHI_M,PHI_H) ! DO I=1,nCols IF(mskant(i) == 1_i8)THEN IF ( FLANDG(I) < 1.0_r8 ) THEN CHV(I) = ( VKMAN / PHI_H(I) ) * V_S(I) CDV(I) = ( VKMAN / PHI_M(I) ) * V_S(I) ENDIF ! SEA_MASK END IF ENDDO ENDDO ! Iteration loop !----------------------------------------------------------------------- !! Set CD's and CH's to be dimensionless paremters !----------------------------------------------------------------------- DO I=1,nCols IF(mskant(i) == 1_i8)THEN IF ( FLANDG(I) < 1.0_r8 ) THEN CDV(I) = (CDV(I) / VSHR(I))*CTRL_TUNABLE CHV(I) = (CHV(I) / VSHR(I))*CTRL_TUNABLE ENDIF ! SEA_MASK END IF ENDDO ! IF (LTIMER) THEN ! DEPENDS ON: timer ! CALL TIMER('FCDCH ',4) ! ENDIF RETURN END SUBROUTINE FCDCH_SEA !! Arguments:--------------------------------------------------------- ! *****************************COPYRIGHT******************************* ! (C) Crown copyright Met Office. All rights reserved. ! For further details please refer to the file COPYRIGHT.txt ! which you should have received as part of this distribution. ! *****************************COPYRIGHT******************************* !!! SUBROUTINES PHI_M_H_SEA AND PHI_M_H_LAND ------------------------ !!! !!! Purpose: Calculate the integrated froms of the Monin-Obukhov !!! stability functions for surface exchanges. !!! !!! !!! Model Modification history: !!! version Date !!! !!! 5.2 15/11/00 New Deck M. Best ! 5.3 25/04/01 Add coastal tiling. Nic Gedney ! 6.0 19/08/03 NEC SX-6 optimisation - force vectorisation of loop ! in phi_m_h_land. R Barnes & J-C Rioual. !!! !!! Programming standard: !!! !!! System component covered: Part of P243. !!! !!! Documentation: UM Documentation Paper No 24. !!! !*L Arguments:--------------------------------------------------------- SUBROUTINE PHI_M_H_SEA( & & nCols,mskant,FLANDG, & & RECIP_L_MO,Z_UV,Z_TQ,Z0M,Z0H,PHI_M,PHI_H) IMPLICIT NONE INTEGER & & nCols ! IN Number of X points? INTEGER(KIND=i8), INTENT(IN ) :: mskant(nCols) REAL(KIND=r8) & & FLANDG(nCols) & ! ! IN Land fraction &,RECIP_L_MO(nCols) & ! ! IN Reciprocal of the Monin-Obukhov length (m^-1). &,Z_UV(nCols) & ! ! IN Height of wind level above roughness height (m &,Z_TQ(nCols) & ! ! IN Height of temperature, moisture and scalar lev ! ! above the roughness height (m). &,Z0M(nCols) & ! IN Roughness length for momentum (m). &,Z0H(nCols) ! ! IN Roughness length for heat/moisture/scalars (m) ! REAL(KIND=r8) & & PHI_M(nCols) & ! ! OUT Stability function for momentum. &,PHI_H(nCols) ! ! OUT Stability function for heat/moisture/scalars. ! !*L Workspace usage---------------------------------------------------- ! No work areas are required. ! !*---------------------------------------------------------------------- !*L External subprograms called: ! EXTERNAL TIMER !*---------------------------------------------------------------------- ! Common and local physical constants. ! REAL(KIND=r8) A,B,C,D,C_OVER_D PARAMETER ( & & A=1.0_r8 & ! &,B=2.0_r8/3.0_r8 & ! Constants used in the Beljaars and &,C=5.0_r8 & ! Holtslag stable stability functions &,D=0.35_r8 & ! &,C_OVER_D=C/D & ! &) ! ! Define local variables. ! INTEGER I ! Loop counter; horizontal field index. ! REAL(KIND=r8) & & PHI_MN & ! Neutral value of stability function for momentum &,PHI_HN & ! Neutral value of stability function for scalars. &,ZETA_UV & ! Temporary in calculation of PHI_M. &,ZETA_0M & ! Temporary in calculation of PHI_M. &,ZETA_TQ & ! Temporary in calculation of PHI_H. &,ZETA_0H & ! Temporary in calculation of PHI_H. &,X_UV_SQ & ! Temporary in calculation of PHI_M. &,X_0M_SQ & ! Temporary in calculation of PHI_M. &,X_UV & ! Temporary in calculation of PHI_M. &,X_0M & ! Temporary in calculation of PHI_M. &,Y_TQ & ! Temporary in calculation of PHI_H. &,Y_0H & ! Temporary in calculation of PHI_H. &,PHI_H_FZ1 & ! Temporary in calculation of PHI_H. &,PHI_H_FZ0 ! Temporary in calculation of PHI_H. ! IF (LTIMER) THEN ! DEPENDS ON: timer ! CALL TIMER('PHI_M_H ',3) ! ENDIF DO I=1,nCols IF(mskant(i) == 1_i8)THEN IF ( FLANDG(I) < 1.0_r8 ) THEN !----------------------------------------------------------------------- !! 1. Calculate neutral values of PHI_M and PHI_H. !----------------------------------------------------------------------- ! PHI_MN = LOG( (Z_UV(I) + Z0M(I)) / Z0M(I) ) PHI_HN = LOG( (Z_TQ(I) + Z0M(I)) / Z0H(I) ) ! !----------------------------------------------------------------------- !! 2. Calculate stability parameters. !----------------------------------------------------------------------- ! ZETA_UV = (Z_UV(I) + Z0M(I)) * RECIP_L_MO(I) ZETA_TQ = (Z_TQ(I) + Z0M(I)) * RECIP_L_MO(I) ZETA_0M = Z0M(I) * RECIP_L_MO(I) ZETA_0H = Z0H(I) * RECIP_L_MO(I) ! !----------------------------------------------------------------------- !! 3. Calculate PHI_M and PHI_H for neutral and stable conditions. !! Formulation of Beljaars and Holtslag (1991). !----------------------------------------------------------------------- ! IF (RECIP_L_MO(I) >= 0.0_r8) THEN PHI_M(I) = PHI_MN & & + A * (ZETA_UV - ZETA_0M) & & + B * ( (ZETA_UV - C_OVER_D) * EXP(-D*ZETA_UV) & & -(ZETA_0M - C_OVER_D) * EXP(-D*ZETA_0M) ) PHI_H_FZ1 = SQRT(1.0_r8 + (2.0_r8/3.0_r8)*A*ZETA_TQ) PHI_H_FZ0 = SQRT(1.0_r8 + (2.0_r8/3.0_r8)*A*ZETA_0H) PHI_H(I) = PHI_HN + & & PHI_H_FZ1*PHI_H_FZ1*PHI_H_FZ1 & & - PHI_H_FZ0*PHI_H_FZ0*PHI_H_FZ0 & & + B * ( (ZETA_TQ - C_OVER_D) * EXP(-D*ZETA_TQ) & & -(ZETA_0H - C_OVER_D) * EXP(-D*ZETA_0H) ) ! !----------------------------------------------------------------------- !! 4. Calculate PHI_M and PHI_H for unstable conditions. !----------------------------------------------------------------------- ! ELSE X_UV_SQ = SQRT(1.0_r8 - 16.0_r8*ZETA_UV) X_0M_SQ = SQRT(1.0_r8 - 16.0_r8*ZETA_0M) X_UV = SQRT(X_UV_SQ) X_0M = SQRT(X_0M_SQ) PHI_M(I) = PHI_MN - 2.0_r8*LOG( (1.0_r8+X_UV) / (1.0_r8+X_0M) ) & & - LOG( (1.0_r8+X_UV_SQ) / (1.0_r8+X_0M_SQ) ) & & + 2.0_r8*( ATAN(X_UV) - ATAN(X_0M) ) Y_TQ = SQRT(1.0_r8 - 16.0_r8*ZETA_TQ) Y_0H = SQRT(1.0_r8 - 16.0_r8*ZETA_0H) PHI_H(I) = PHI_HN - 2.0_r8*LOG( (1.0_r8+Y_TQ) / (1.0_r8+Y_0H) ) ENDIF ENDIF ! SEA_MASK END IF ENDDO ! IF (LTIMER) THEN ! DEPENDS ON: timer ! CALL TIMER('PHI_M_H ',4) ! ENDIF RETURN END SUBROUTINE PHI_M_H_SEA !!! !*L Arguments:--------------------------------------------------------- ! *****************************COPYRIGHT******************************* ! (C) Crown copyright Met Office. All rights reserved. ! For further details please refer to the file COPYRIGHT.txt ! which you should have received as part of this distribution. ! *****************************COPYRIGHT******************************* ! !!! SUBROUTINES SF_RIB_LAND and SF_RIB_SEA --------------------------- !!! !!! Purpose: Calculate bulk Richardson number for surface layer !!! !!! !!! Model Modification history !!! version date !!! 5.2 15/11/00 New Deck M. Best ! 5.3 25/04/01 Add coastal tiling. Nic Gedney !!! !!! Programming standard: !!! !!! ------------------------------------------------------------------ ! SUBROUTINE SF_RIB_LAND-------------------------------------------- ! ! Calculate RIB for land tiles ! ! ------------------------------------------------------------------ ! SUBROUTINE SF_RIB_SEA--------------------------------------------- ! ! Calculate RIB for sea, sea-ice and sea-ice leads ! ! ------------------------------------------------------------------ SUBROUTINE SF_RIB_SEA ( & nCols, &!INTEGER, INTENT(IN ) :: nCols ! IN Number of X points? mskant,& FLANDG , &!REAL(KIND=r8) , INTENT(IN ) :: FLANDG(nCols,ROWS) ! IN Land fraction on all pts. NSICE , &!INTEGER, INTENT(IN ) :: NSICE ! IN Number of sea-ice points. SICE_INDEX, &!INTEGER, INTENT(IN ) :: SICE_INDEX(nCols*ROWS,2) ! IN Index of sea-ice points. BQ_1 , &!REAL(KIND=r8) , INTENT(IN ) :: BQ_1(nCols,ROWS) ! IN A buoyancy parameter for lowest atm level. ("beta-q twiddle"). BT_1 , &!REAL(KIND=r8) , INTENT(IN ) :: BT_1(nCols,ROWS)! IN A buoyancy parameter for lowest atm level. ("beta-T twiddle"). QSTAR_ICE , &!REAL(KIND=r8) , INTENT(IN ) :: QSTAR_ICE(nCols,ROWS)! IN Surface saturated sp humidity over sea-ice. QSTAR_SEA , &!REAL(KIND=r8) , INTENT(IN ) :: QSTAR_SEA(nCols,ROWS)! IN Surface saturated sp humidity over sea and sea-ice leads. QW_1 , &!REAL(KIND=r8) , INTENT(IN ) :: QW_1(nCols,ROWS)! IN Total water content of lowest atmospheric layer (kg per kg air) TL_1 , &!REAL(KIND=r8) , INTENT(IN ) :: TL_1(nCols,ROWS)! IN Liquid/frozen water temperature for lowest atmospheric layer (K). TSTAR_SICE, &!REAL(KIND=r8) , INTENT(IN ) :: TSTAR_SICE(nCols,ROWS)! IN Surface temperature of sea-ice (K). TSTAR_SEA , &!REAL(KIND=r8) , INTENT(IN ) :: TSTAR_SEA(nCols,ROWS)! IN Surface temperature of sea and sea-ice leads (K). VSHR , &!REAL(KIND=r8) , INTENT(IN ) :: VSHR(nCols,ROWS)! IN Magnitude of surface- to-lowest-level wind shear. Z0H_ICE , &!REAL(KIND=r8) , INTENT(IN ) :: Z0H_ICE(nCols,ROWS)! IN Roughness length for heat and moisture transport over sea-ice (m). Z0H_SEA , &!REAL(KIND=r8) , INTENT(IN ) :: Z0H_SEA(nCols,ROWS)! IN Roughness length for heat and moisture transport over sea or sea-ice leads (m). Z0M_ICE , &!REAL(KIND=r8) , INTENT(IN ) :: Z0M_ICE(nCols,ROWS)! IN Roughness length for momentum over sea-ice (m). Z0M_SEA , &!REAL(KIND=r8) , INTENT(IN ) :: Z0M_SEA(nCols,ROWS)! IN Roughness length for momentum over sea or sea-ice leads (m). Z1_TQ , &!REAL(KIND=r8) , INTENT(IN ) :: Z1_TQ(nCols,ROWS)! IN Height of lowest TQ level (m). Z1_UV , &!REAL(KIND=r8) , INTENT(IN ) :: Z1_UV(nCols,ROWS) ! IN Height of lowest UV level (m). RIB_SEA , &!REAL(KIND=r8) , INTENT(OUT ) :: RIB_SEA(nCols,ROWS)! OUT Bulk Richardson number for lowest layer over sea or sea-ice leads. RIB_ICE , &!REAL(KIND=r8) , INTENT(OUT ) :: RIB_ICE(nCols,ROWS)! OUT Bulk Richardson number for lowest layer over sea-ice. DB_SEA , &!REAL(KIND=r8) , INTENT(OUT ) :: DB_SEA(nCols,ROWS)! OUT Buoyancy difference between DB_ICE &!REAL(KIND=r8) , INTENT(OUT ) :: DB_ICE(nCols,ROWS) ! OUT Buoyancy difference between ) IMPLICIT NONE INTEGER, INTENT(IN ) :: nCols ! IN Number of X points? INTEGER, INTENT(IN ) :: NSICE ! IN Number of sea-ice points. INTEGER, INTENT(IN ) :: SICE_INDEX(nCols,2) ! IN Index of sea-ice points. INTEGER(KIND=i8), INTENT(IN ) :: mskant(nCols) REAL(KIND=r8) , INTENT(IN ) :: FLANDG(nCols) ! IN Land fraction on all pts. REAL(KIND=r8) , INTENT(IN ) :: BQ_1(nCols) ! IN A buoyancy parameter for lowest atm level. ("beta-q twiddle"). REAL(KIND=r8) , INTENT(IN ) :: BT_1(nCols)! IN A buoyancy parameter for lowest atm level. ("beta-T twiddle"). REAL(KIND=r8) , INTENT(IN ) :: QSTAR_ICE(nCols)! IN Surface saturated sp humidity over sea-ice. REAL(KIND=r8) , INTENT(IN ) :: QSTAR_SEA(nCols)! IN Surface saturated sp humidity over sea and sea-ice leads. REAL(KIND=r8) , INTENT(IN ) :: QW_1(nCols)! IN Total water content of lowest atmospheric layer (kg per kg air) REAL(KIND=r8) , INTENT(IN ) :: TL_1(nCols)! IN Liquid/frozen water temperature for lowest atmospheric layer (K). REAL(KIND=r8) , INTENT(IN ) :: TSTAR_SICE(nCols)! IN Surface temperature of sea-ice (K). REAL(KIND=r8) , INTENT(IN ) :: TSTAR_SEA(nCols)! IN Surface temperature of sea and sea-ice leads (K). REAL(KIND=r8) , INTENT(IN ) :: VSHR(nCols)! IN Magnitude of surface- to-lowest-level wind shear. REAL(KIND=r8) , INTENT(IN ) :: Z0H_ICE(nCols)! IN Roughness length for heat and moisture transport over sea-ice (m). REAL(KIND=r8) , INTENT(IN ) :: Z0H_SEA(nCols)! IN Roughness length for heat and moisture transport over sea or sea-ice leads (m). REAL(KIND=r8) , INTENT(IN ) :: Z0M_ICE(nCols)! IN Roughness length for momentum over sea-ice (m). REAL(KIND=r8) , INTENT(IN ) :: Z0M_SEA(nCols)! IN Roughness length for momentum over sea or sea-ice leads (m). REAL(KIND=r8) , INTENT(IN ) :: Z1_TQ(nCols)! IN Height of lowest TQ level (m). REAL(KIND=r8) , INTENT(IN ) :: Z1_UV(nCols) ! IN Height of lowest UV level (m). REAL(KIND=r8) , INTENT(OUT ) :: RIB_SEA(nCols)! OUT Bulk Richardson number for lowest layer over sea or sea-ice leads. REAL(KIND=r8) , INTENT(OUT ) :: RIB_ICE(nCols)! OUT Bulk Richardson number for lowest layer over sea-ice. REAL(KIND=r8) , INTENT(OUT ) :: DB_SEA(nCols)! OUT Buoyancy difference between ! ! surface and lowest atmospheric ! ! level over sea or sea-ice leads. REAL(KIND=r8) , INTENT(OUT ) :: DB_ICE(nCols) ! OUT Buoyancy difference between ! ! surface and lowest atmospheric ! ! level over sea-ice. ! External routines called :- ! EXTERNAL TIMER ! Symbolic constants ----------------------------------------------- !*L------------------COMDECK C_O_DG_C----------------------------------- ! ZERODEGC IS CONVERSION BETWEEN DEGREES CELSIUS AND KELVIN ! TFS IS TEMPERATURE AT WHICH SEA WATER FREEZES ! TM IS TEMPERATURE AT WHICH FRESH WATER FREEZES AND ICE MELTS REAL(KIND=r8), PARAMETER :: ZeroDegC = 273.15_r8 REAL(KIND=r8), PARAMETER :: TFS = 271.35_r8 REAL(KIND=r8), PARAMETER :: TM = 273.15_r8 !*---------------------------------------------------------------------- !*L------------------COMDECK C_G---------------------------------------- ! G IS MEAN ACCEL DUE TO GRAVITY AT EARTH'S SURFACE REAL(KIND=r8), PARAMETER :: G = 9.80665_r8 !*---------------------------------------------------------------------- !*L------------------COMDECK C_R_CP------------------------------------- ! History: ! Version Date Comment. ! 5.0 07/05/99 Add variable P_zero for consistency with ! conversion to C-P 'C' dynamics grid. R. Rawlins ! 5.1 07/03/00 Fixed/Free format conversion P. Selwood ! R IS GAS CONSTANT FOR DRY AIR ! CP IS SPECIFIC HEAT OF DRY AIR AT CONSTANT PRESSURE ! PREF IS REFERENCE SURFACE PRESSURE REAL(KIND=r8), PARAMETER :: R = 287.05_r8 REAL(KIND=r8), PARAMETER :: CP = 1005.0_r8 REAL(KIND=r8), PARAMETER :: Kappa = R/CP REAL(KIND=r8), PARAMETER :: Pref = 100000.0_r8 ! Reference surface pressure = PREF REAL(KIND=r8), PARAMETER :: P_zero = Pref REAL(KIND=r8), PARAMETER :: sclht = 6.8E+03_r8 ! Scale Height H !*---------------------------------------------------------------------- ! Workspace -------------------------------------------------------- INTEGER & & I & ! Horizontal field index. &,SI !Sea-ice field index. REAL(KIND=r8) & & DQ & ! Sp humidity difference between surface ! ! and lowest atmospheric level (Q1 - Q*). &,DTEMP ! Modified temperature difference between ! ! surface and lowest atmospheric level. ! IF (LTIMER) THEN ! DEPENDS ON: timer ! CALL TIMER('SF_RIB ',3) ! ENDIF DO I=1,nCols IF(mskant(i) == 1_i8)THEN IF ( FLANDG(I) < 1.0_r8 ) THEN ! Sea and sea-ice leads DTEMP = TL_1(I) - TSTAR_SEA(I) & ! P243.118 & + (G/CP)*(Z1_TQ(I) + Z0M_SEA(I) - Z0H_SEA(I)) DQ = QW_1(I) - QSTAR_SEA(I) ! P243.119 DB_SEA(I) = G*( BT_1(I)*DTEMP + BQ_1(I)*DQ ) RIB_SEA(I) = Z1_UV(I)*DB_SEA(I) / & & ( VSHR(I)*VSHR(I) ) ENDIF END IF ENDDO DO SI=1,NSICE IF(mskant(i) == 1_i8)THEN I = SICE_INDEX(SI,1) ! Sea-ice DTEMP = TL_1(I) - TSTAR_SICE(I) & & + (G/CP)*(Z1_TQ(I) + Z0M_ICE(I) - Z0H_ICE(I)) DQ = QW_1(I) - QSTAR_ICE(I) DB_ICE(I) = G*( BT_1(I)*DTEMP + BQ_1(I)*DQ ) RIB_ICE(I) = Z1_UV(I)*DB_ICE(I) / & & ( VSHR(I) * VSHR(I) ) END IF ENDDO ! IF (LTIMER) THEN ! DEPENDS ON: timer ! !CALL TIMER('SF_RIB ',4) ! ENDIF RETURN END SUBROUTINE SF_RIB_SEA ! *****************************COPYRIGHT******************************* ! (C) Crown copyright Met Office. All rights reserved. ! For further details please refer to the file COPYRIGHT.txt ! which you should have received as part of this distribution. ! *****************************COPYRIGHT******************************* ! ! SUBROUTINE BOUY_TQ ! PURPOSE: To calculate buoyancy parameters on p,T,q-levels ! ! METHOD: ! ! HISTORY: ! DATE VERSION COMMENT ! ---- ------- ------- ! new deck !!! 4.5 Jul. 98 Kill the IBM specific lines. (JCThil) !LL 5.2 27/09/00 change from QSAT_2D to QSAT A.Malcolm !!! 5.3 Feb. 01 Calculate grid-box mean parameters (APLock) ! 6.1 17/05/04 Change Q, QCL, QCF dims to enable substepping. ! M. Diamantakis ! ! CODE DESCRIPTION: ! LANGUAGE: FORTRAN 77 + CRAY EXTENSIONS ! THIS CODE IS WRITTEN TO UMDP 3 PROGRAMMING STANDARDS. ! SUBROUTINE BOUY_TQ (& nCols, & !INTEGER, INTENT(IN ) :: nCols kMAx , & !INTEGER, INTENT(IN ) :: kMAx LQ_MIX_BL , & !LOGICAL, INTENT(IN ) :: LQ_MIX_BL mskant , & P , & !REAL(KIND=r8), INTENT(IN ) :: P(nCols,rows,kMAx) T , & !REAL(KIND=r8), INTENT(IN ) ::T(nCols,rows,kMAx) Q , & !REAL(KIND=r8), INTENT(IN ) ::Q(nCols,rows,kMAx) QCF , & !REAL(KIND=r8), INTENT(IN ) ::QCF(nCols,rows,kMAx) QCL , & !REAL(KIND=r8), INTENT(IN ) ::QCL(nCols,rows,kMAx) CF , & !REAL(KIND=r8), INTENT(IN ) ::CF(nCols, rows, kMAx) BT , & !REAL(KIND=r8), INTENT(out ) :: BT(nCols,rows,kMAx) BQ , & !REAL(KIND=r8), INTENT(out ) :: BQ(nCols,rows,kMAx) BT_CLD , & !REAL(KIND=r8), INTENT(out ) :: BT_CLD(nCols,rows,kMAx) BQ_CLD , & !REAL(KIND=r8), INTENT(out ) :: BQ_CLD(nCols,rows,kMAx) BT_GB , & !REAL(KIND=r8), INTENT(out ) :: BT_GB(nCols,rows,kMAx) BQ_GB , & !REAL(KIND=r8), INTENT(out ) :: BQ_GB(nCols,rows,kMAx) A_QS , & !REAL(KIND=r8), INTENT(out ) :: A_QS(nCols,rows,kMAx) A_DQSDT , & !REAL(KIND=r8), INTENT(out ) :: A_DQSDT(nCols,rows,kMAx) DQSDT & !REAL(KIND=r8), INTENT(out ) :: DQSDT(nCols,rows,kMAx) ) IMPLICIT NONE ! ARGUMENTS WITH INTENT IN. IE: INPUT VARIABLES. LOGICAL, INTENT(IN ) :: LQ_MIX_BL ! IN switch for using mixing ratios INTEGER, INTENT(IN ) :: nCols INTEGER, INTENT(IN ) :: kMAx ! IN No. of atmospheric levels for which ! boundary layer fluxes are calculated. ! Assumed <=30 for dimensioning GAMMA() ! in common deck C_GAMMA INTEGER(KIND=i8), INTENT(IN ) :: mskant(nCols) REAL(KIND=r8), INTENT(IN ) :: P(nCols,kMAx) ! IN Pressure at pressure points. REAL(KIND=r8), INTENT(IN ) ::T(nCols,kMAx) ! IN Temperature (K). At P points REAL(KIND=r8), INTENT(IN ) ::Q(nCols,kMAx) ! IN Sp humidity (kg water per kg air). REAL(KIND=r8), INTENT(IN ) ::QCL(nCols,kMAx) ! IN Cloud liq water (kg per kg air). REAL(KIND=r8), INTENT(IN ) ::QCF(nCols,kMAx) ! IN Cloud liq water (kg per kg air). REAL(KIND=r8), INTENT(IN ) ::CF(nCols, kMAx)! IN Cloud fraction (decimal). ! ARGUMENTS WITH INTENT OUT. IE: OUTPUT VARIABLES. REAL(KIND=r8), INTENT(out ) :: BQ(nCols,kMAx) ! OUT A buoyancy parameter for clear air REAL(KIND=r8), INTENT(out ) :: BT(nCols,kMAx) ! OUT A buoyancy parameter for clear air REAL(KIND=r8), INTENT(out ) :: BQ_CLD(nCols,kMAx) ! ! OUT A buoyancy parameter for cloudy air REAL(KIND=r8), INTENT(out ) :: BT_CLD(nCols,kMAx) ! ! OUT A buoyancy parameter for cloudy air REAL(KIND=r8), INTENT(out ) :: BQ_GB(nCols,kMAx) ! OUT A grid-box mean buoyancy parameter REAL(KIND=r8), INTENT(out ) :: BT_GB(nCols,kMAx) ! OUT A grid-box mean buoyancy parameter REAL(KIND=r8), INTENT(out ) :: A_QS(nCols,kMAx) ! ! OUT Saturated lapse rate factor REAL(KIND=r8), INTENT(out ) :: A_DQSDT(nCols,kMAx) ! ! OUT Saturated lapse rate factor REAL(KIND=r8), INTENT(out ) :: DQSDT(nCols,kMAx) ! ! OUT Derivative of q_SAT w.r.t. T ! LOCAL VARIABLES. REAL(KIND=r8) :: qss(nCols) REAL(KIND=r8) & & QS(nCols) ! WORK Saturated mixing ratio. INTEGER & & I & &, K REAL(KIND=r8) & & BC ! EXTERNAL & ! & QSAT, TIMER !*L------------------COMDECK C_O_DG_C----------------------------------- ! ZERODEGC IS CONVERSION BETWEEN DEGREES CELSIUS AND KELVIN ! TFS IS TEMPERATURE AT WHICH SEA WATER FREEZES ! TM IS TEMPERATURE AT WHICH FRESH WATER FREEZES AND ICE MELTS REAL(KIND=r8), PARAMETER :: ZeroDegC = 273.15_r8 REAL(KIND=r8), PARAMETER :: TFS = 271.35_r8 REAL(KIND=r8), PARAMETER :: TM = 273.15_r8 !*---------------------------------------------------------------------- ! C_LHEAT start ! latent heat of condensation of water at 0degc REAL(KIND=r8),PARAMETER:: LC=2.501E6_r8 ! latent heat of fusion at 0degc REAL(KIND=r8),PARAMETER:: LF=0.334E6_r8 ! C_LHEAT end !*L------------------COMDECK C_G---------------------------------------- ! G IS MEAN ACCEL DUE TO GRAVITY AT EARTH'S SURFACE REAL(KIND=r8), PARAMETER :: G = 9.80665_r8 !*---------------------------------------------------------------------- !*L------------------COMDECK C_R_CP------------------------------------- ! History: ! Version Date Comment. ! 5.0 07/05/99 Add variable P_zero for consistency with ! conversion to C-P 'C' dynamics grid. R. Rawlins ! 5.1 07/03/00 Fixed/Free format conversion P. Selwood ! R IS GAS CONSTANT FOR DRY AIR ! CP IS SPECIFIC HEAT OF DRY AIR AT CONSTANT PRESSURE ! PREF IS REFERENCE SURFACE PRESSURE REAL(KIND=r8), PARAMETER :: R = 287.05_r8 REAL(KIND=r8), PARAMETER :: CP = 1005.0_r8 REAL(KIND=r8), PARAMETER :: Kappa = R/CP REAL(KIND=r8), PARAMETER :: Pref = 100000.0_r8 ! Reference surface pressure = PREF REAL(KIND=r8), PARAMETER :: P_zero = Pref REAL(KIND=r8), PARAMETER :: sclht = 6.8E+03_r8 ! Scale Height H !*---------------------------------------------------------------------- !*L------------------COMDECK C_EPSLON----------------------------------- ! EPSILON IS RATIO OF MOLECULAR WEIGHTS OF WATER AND DRY AIR REAL(KIND=r8), PARAMETER :: Epsilon = 0.62198_r8 REAL(KIND=r8), PARAMETER :: C_Virtual = 1.0_r8/Epsilon-1.0_r8 !*---------------------------------------------------------------------- ! C_VKMAN start REAL(KIND=r8),PARAMETER:: VKMAN=0.4_r8 ! Von Karman's constant ! C_VKMAN end ! C_SOILH start ! No. of soil layers (must = NSOIL). REAL(KIND=r8),PARAMETER:: PSOIL=4.0_r8 ! Tunable characteristic freq (rad/s) REAL(KIND=r8),PARAMETER:: OMEGA1=3.55088E-4_r8 ! Density of lying snow (kg per m**3) REAL(KIND=r8),PARAMETER:: RHO_SNOW=250.0_r8 ! Depth of `effective' snow surface layer (m) REAL(KIND=r8),PARAMETER:: DEFF_SNOW=0.1_r8 ! Thermal conductivity of lying snow (Watts per m per K). REAL(KIND=r8),PARAMETER:: SNOW_HCON=0.265_r8 ! Thermal capacity of lying snow (J/K/m3) REAL(KIND=r8),PARAMETER:: SNOW_HCAP=0.63E6_r8 ! C_SOILH end ! History: ! Version Date Comment ! 3.4 18/5/94 Add PP missing data indicator. J F Thomson ! 5.1 6/3/00 Convert to Free/Fixed format. P Selwood !*L------------------COMDECK C_MDI------------------------------------- ! PP missing data indicator (-1.0E+30) REAL(KIND=r8), PARAMETER :: RMDI_PP = -1.0E+30_r8 ! Old REAL(KIND=r8) missing data indicator (-32768.0) REAL(KIND=r8), PARAMETER :: RMDI_OLD = -32768.0_r8 ! New REAL(KIND=r8) missing data indicator (-2**30) REAL(KIND=r8), PARAMETER :: RMDI = -32768.0_r8*32768.0_r8 ! Integer missing data indicator INTEGER, PARAMETER :: IMDI = -32768 !*---------------------------------------------------------------------- REAL(KIND=r8) ETAR,GRCP,LCRCP,LFRCP,LS,LSRCP PARAMETER ( & & ETAR=1.0_r8/(1.0_r8-EPSILON) & ! Used in buoyancy parameter BETAC. &,GRCP=G/CP & ! Used in DZTL, FTL calculations. &,LCRCP=LC/CP & ! Latent heat of condensation / CP. &,LFRCP=LF/CP & ! Latent heat of fusion / CP. &,LS=LC+LF & ! Latent heat of sublimation. &,LSRCP=LS/CP & ! Latent heat of sublimation / CP. &) ! IF (LTIMER) THEN !! DEPENDS ON: timer ! CALL TIMER('BOUY_TQ ',3) ! ENDIF !----------------------------------------------------------------------- !! 1. Loop round levels. !----------------------------------------------------------------------- DO K=1,kMAx !----------------------------------------------------------------------- !! 1.1 Calculate saturated specific humidity at pressure and !! temperature of current level. !----------------------------------------------------------------------- ! DEPENDS ON: qsat_mix IF(QSAT_opt)THEN CALL QSAT_mix(QS,T(1,K),P(1,K),nCols,Lq_mix_bl,mskant) END IF ! DO I=1,nCols IF(mskant(i) == 1_i8)THEN IF(.NOT. QSAT_opt)THEN qss(i) = fpvs(T(i,k)) QS (i) = EPS * qss(i) / (P(i,k) + EPSM1 * qss(i)) END IF !ajm DO I=P1,P1+P_POINTS-1 !----------------------------------------------------------------------- !! 1.2 Calculate buoyancy parameters BT and BQ, required for the !! calculation of stability. !----------------------------------------------------------------------- BT(I,K) = 1.0_r8/T(I,K) BQ(I,K) = & & C_VIRTUAL/(1.0_r8+C_VIRTUAL*Q(I,K)-QCL(I,K)-QCF(I,K)) ! IF (T(I,K) > TM) THEN DQSDT(I,K) = (EPSILON * LC * QS(I)) & & / ( R * T(I,K) * T(I,K) ) ! ... (Clausius-Clapeyron) for T above freezing ! A_QS(I,K) = 1.0_r8 / (1.0_r8 + LCRCP*DQSDT(I,K)) ! A_DQSDT(I,K) = A_QS(I,K) * DQSDT(I,K) ! BC = LCRCP*BT(I,K) - ETAR*BQ(I,K) ! ELSE DQSDT(I,K) = (EPSILON * LS * QS(I)) & & / ( R * T(I,K) * T(I,K) ) ! ... (Clausius-Clapeyron) for T below freezing ! A_QS(I,K) = 1.0_r8 / (1.0_r8 + LSRCP*DQSDT(I,K)) ! A_DQSDT(I,K) = A_QS(I,K) * DQSDT(I,K) ! BC = LSRCP*BT(I,K) - ETAR*BQ(I,K) ! ENDIF ! !----------------------------------------------------------------------- !! 1.3 Calculate in-cloud buoyancy parameters. !----------------------------------------------------------------------- ! BT_CLD(I,K) = BT(I,K) - A_DQSDT(I,K) * BC BQ_CLD(I,K) = BQ(I,K) + A_QS(I,K) * BC !----------------------------------------------------------------------- !! 1.4 Calculate grid-box mean buoyancy parameters. !----------------------------------------------------------------------- ! BT_GB(I,K) = BT(I,K) + CF(I,K)*( BT_CLD(I,K) - BT(I,K) ) BQ_GB(I,K) = BQ(I,K) + CF(I,K)*( BQ_CLD(I,K) - BQ(I,K) ) ! END IF ENDDO ! p_points,j ENDDO ! kMAx ! IF (LTIMER) THEN !! DEPENDS ON: timer ! CALL TIMER('BOUY_TQ ',4) ! ENDIF RETURN END SUBROUTINE BOUY_TQ ! *****************************COPYRIGHT******************************* ! (C) Crown copyright Met Office. All rights reserved. ! For further details please refer to the file COPYRIGHT.txt ! which you should have received as part of this distribution. ! *****************************COPYRIGHT******************************* ! SUBROUTINE qsat_data() ! ! Local dynamic arrays------------------------------------------------- REAL(KIND=r8) ES_Aux(0:N+1) ! TABLE OF SATURATION WATER VAPOUR PRESSURE (PA) ! - SET BY DATA STATEMENT CALCULATED FROM THE ! GOFF-GRATCH FORMULAE AS TAKEN FROM LANDOLT- ! BORNSTEIN, 1987 NUMERICAL DATA AND FUNCTIONAL ! RELATIONSHIPS IN SCIENCE AND TECHNOLOGY. ! GROUP V/ VOL 4B METEOROLOGY. PHYSICAL AND ! CHEMICAL PROPERTIES OF AIR, P35 ! ! ---------------------------------------------------------------------- ! SATURATION WATER VAPOUR PRESSURE ! ! ABOVE 0 DEG C VALUES ARE OVER WATER ! ! BELOW 0 DEC C VALUES ARE OVER ICE ! ---------------------------------------------------------------------- ! Note: 0 element is a repeat of 1st element to cater for special case ! of low temperatures ( <= T_LOW) for which the array index is ! rounded down due to machine precision. DATA (ES_Aux(IES),IES= 0, 95) / 0.966483E-02_r8, & 0.966483E-02_r8,0.984279E-02_r8,0.100240E-01_r8,0.102082E-01_r8,0.103957E-01_r8, & 0.105865E-01_r8,0.107803E-01_r8,0.109777E-01_r8,0.111784E-01_r8,0.113825E-01_r8, & 0.115902E-01_r8,0.118016E-01_r8,0.120164E-01_r8,0.122348E-01_r8,0.124572E-01_r8, & 0.126831E-01_r8,0.129132E-01_r8,0.131470E-01_r8,0.133846E-01_r8,0.136264E-01_r8, & 0.138724E-01_r8,0.141225E-01_r8,0.143771E-01_r8,0.146356E-01_r8,0.148985E-01_r8, & 0.151661E-01_r8,0.154379E-01_r8,0.157145E-01_r8,0.159958E-01_r8,0.162817E-01_r8, & 0.165725E-01_r8,0.168680E-01_r8,0.171684E-01_r8,0.174742E-01_r8,0.177847E-01_r8, & 0.181008E-01_r8,0.184216E-01_r8,0.187481E-01_r8,0.190801E-01_r8,0.194175E-01_r8, & 0.197608E-01_r8,0.201094E-01_r8,0.204637E-01_r8,0.208242E-01_r8,0.211906E-01_r8, & 0.215631E-01_r8,0.219416E-01_r8,0.223263E-01_r8,0.227172E-01_r8,0.231146E-01_r8, & 0.235188E-01_r8,0.239296E-01_r8,0.243465E-01_r8,0.247708E-01_r8,0.252019E-01_r8, & 0.256405E-01_r8,0.260857E-01_r8,0.265385E-01_r8,0.269979E-01_r8,0.274656E-01_r8, & 0.279405E-01_r8,0.284232E-01_r8,0.289142E-01_r8,0.294124E-01_r8,0.299192E-01_r8, & 0.304341E-01_r8,0.309571E-01_r8,0.314886E-01_r8,0.320285E-01_r8,0.325769E-01_r8, & 0.331348E-01_r8,0.337014E-01_r8,0.342771E-01_r8,0.348618E-01_r8,0.354557E-01_r8, & 0.360598E-01_r8,0.366727E-01_r8,0.372958E-01_r8,0.379289E-01_r8,0.385717E-01_r8, & 0.392248E-01_r8,0.398889E-01_r8,0.405633E-01_r8,0.412474E-01_r8,0.419430E-01_r8, & 0.426505E-01_r8,0.433678E-01_r8,0.440974E-01_r8,0.448374E-01_r8,0.455896E-01_r8, & 0.463545E-01_r8,0.471303E-01_r8,0.479191E-01_r8,0.487190E-01_r8,0.495322E-01_r8/ DATA (ES_Aux(IES),IES= 96,190) / & 0.503591E-01_r8,0.511977E-01_r8,0.520490E-01_r8,0.529145E-01_r8,0.537931E-01_r8, & 0.546854E-01_r8,0.555924E-01_r8,0.565119E-01_r8,0.574467E-01_r8,0.583959E-01_r8, & 0.593592E-01_r8,0.603387E-01_r8,0.613316E-01_r8,0.623409E-01_r8,0.633655E-01_r8, & 0.644053E-01_r8,0.654624E-01_r8,0.665358E-01_r8,0.676233E-01_r8,0.687302E-01_r8, & 0.698524E-01_r8,0.709929E-01_r8,0.721490E-01_r8,0.733238E-01_r8,0.745180E-01_r8, & 0.757281E-01_r8,0.769578E-01_r8,0.782061E-01_r8,0.794728E-01_r8,0.807583E-01_r8, & 0.820647E-01_r8,0.833905E-01_r8,0.847358E-01_r8,0.861028E-01_r8,0.874882E-01_r8, & 0.888957E-01_r8,0.903243E-01_r8,0.917736E-01_r8,0.932464E-01_r8,0.947407E-01_r8, & 0.962571E-01_r8,0.977955E-01_r8,0.993584E-01_r8,0.100942E+00_r8,0.102551E+00_r8, & 0.104186E+00_r8,0.105842E+00_r8,0.107524E+00_r8,0.109231E+00_r8,0.110963E+00_r8, & 0.112722E+00_r8,0.114506E+00_r8,0.116317E+00_r8,0.118153E+00_r8,0.120019E+00_r8, & 0.121911E+00_r8,0.123831E+00_r8,0.125778E+00_r8,0.127755E+00_r8,0.129761E+00_r8, & 0.131796E+00_r8,0.133863E+00_r8,0.135956E+00_r8,0.138082E+00_r8,0.140241E+00_r8, & 0.142428E+00_r8,0.144649E+00_r8,0.146902E+00_r8,0.149190E+00_r8,0.151506E+00_r8, & 0.153859E+00_r8,0.156245E+00_r8,0.158669E+00_r8,0.161126E+00_r8,0.163618E+00_r8, & 0.166145E+00_r8,0.168711E+00_r8,0.171313E+00_r8,0.173951E+00_r8,0.176626E+00_r8, & 0.179342E+00_r8,0.182096E+00_r8,0.184893E+00_r8,0.187724E+00_r8,0.190600E+00_r8, & 0.193518E+00_r8,0.196473E+00_r8,0.199474E+00_r8,0.202516E+00_r8,0.205604E+00_r8, & 0.208730E+00_r8,0.211905E+00_r8,0.215127E+00_r8,0.218389E+00_r8,0.221701E+00_r8/ DATA (ES_Aux(IES),IES=191,285) / & 0.225063E+00_r8,0.228466E+00_r8,0.231920E+00_r8,0.235421E+00_r8,0.238976E+00_r8, & 0.242580E+00_r8,0.246232E+00_r8,0.249933E+00_r8,0.253691E+00_r8,0.257499E+00_r8, & 0.261359E+00_r8,0.265278E+00_r8,0.269249E+00_r8,0.273274E+00_r8,0.277358E+00_r8, & 0.281498E+00_r8,0.285694E+00_r8,0.289952E+00_r8,0.294268E+00_r8,0.298641E+00_r8, & 0.303078E+00_r8,0.307577E+00_r8,0.312135E+00_r8,0.316753E+00_r8,0.321440E+00_r8, & 0.326196E+00_r8,0.331009E+00_r8,0.335893E+00_r8,0.340842E+00_r8,0.345863E+00_r8, & 0.350951E+00_r8,0.356106E+00_r8,0.361337E+00_r8,0.366636E+00_r8,0.372006E+00_r8, & 0.377447E+00_r8,0.382966E+00_r8,0.388567E+00_r8,0.394233E+00_r8,0.399981E+00_r8, & 0.405806E+00_r8,0.411714E+00_r8,0.417699E+00_r8,0.423772E+00_r8,0.429914E+00_r8, & 0.436145E+00_r8,0.442468E+00_r8,0.448862E+00_r8,0.455359E+00_r8,0.461930E+00_r8, & 0.468596E+00_r8,0.475348E+00_r8,0.482186E+00_r8,0.489124E+00_r8,0.496160E+00_r8, & 0.503278E+00_r8,0.510497E+00_r8,0.517808E+00_r8,0.525224E+00_r8,0.532737E+00_r8, & 0.540355E+00_r8,0.548059E+00_r8,0.555886E+00_r8,0.563797E+00_r8,0.571825E+00_r8, & 0.579952E+00_r8,0.588198E+00_r8,0.596545E+00_r8,0.605000E+00_r8,0.613572E+00_r8, & 0.622255E+00_r8,0.631059E+00_r8,0.639962E+00_r8,0.649003E+00_r8,0.658144E+00_r8, & 0.667414E+00_r8,0.676815E+00_r8,0.686317E+00_r8,0.695956E+00_r8,0.705728E+00_r8, & 0.715622E+00_r8,0.725641E+00_r8,0.735799E+00_r8,0.746082E+00_r8,0.756495E+00_r8, & 0.767052E+00_r8,0.777741E+00_r8,0.788576E+00_r8,0.799549E+00_r8,0.810656E+00_r8, & 0.821914E+00_r8,0.833314E+00_r8,0.844854E+00_r8,0.856555E+00_r8,0.868415E+00_r8/ DATA (ES_Aux(IES),IES=286,380) / & 0.880404E+00_r8,0.892575E+00_r8,0.904877E+00_r8,0.917350E+00_r8,0.929974E+00_r8, & 0.942771E+00_r8,0.955724E+00_r8,0.968837E+00_r8,0.982127E+00_r8,0.995600E+00_r8, & 0.100921E+01_r8,0.102304E+01_r8,0.103700E+01_r8,0.105116E+01_r8,0.106549E+01_r8, & 0.108002E+01_r8,0.109471E+01_r8,0.110962E+01_r8,0.112469E+01_r8,0.113995E+01_r8, & 0.115542E+01_r8,0.117107E+01_r8,0.118693E+01_r8,0.120298E+01_r8,0.121923E+01_r8, & 0.123569E+01_r8,0.125234E+01_r8,0.126923E+01_r8,0.128631E+01_r8,0.130362E+01_r8, & 0.132114E+01_r8,0.133887E+01_r8,0.135683E+01_r8,0.137500E+01_r8,0.139342E+01_r8, & 0.141205E+01_r8,0.143091E+01_r8,0.145000E+01_r8,0.146933E+01_r8,0.148892E+01_r8, & 0.150874E+01_r8,0.152881E+01_r8,0.154912E+01_r8,0.156970E+01_r8,0.159049E+01_r8, & 0.161159E+01_r8,0.163293E+01_r8,0.165452E+01_r8,0.167640E+01_r8,0.169852E+01_r8, & 0.172091E+01_r8,0.174359E+01_r8,0.176653E+01_r8,0.178977E+01_r8,0.181332E+01_r8, & 0.183709E+01_r8,0.186119E+01_r8,0.188559E+01_r8,0.191028E+01_r8,0.193524E+01_r8, & 0.196054E+01_r8,0.198616E+01_r8,0.201208E+01_r8,0.203829E+01_r8,0.206485E+01_r8, & 0.209170E+01_r8,0.211885E+01_r8,0.214637E+01_r8,0.217424E+01_r8,0.220242E+01_r8, & 0.223092E+01_r8,0.225979E+01_r8,0.228899E+01_r8,0.231855E+01_r8,0.234845E+01_r8, & 0.237874E+01_r8,0.240937E+01_r8,0.244040E+01_r8,0.247176E+01_r8,0.250349E+01_r8, & 0.253560E+01_r8,0.256814E+01_r8,0.260099E+01_r8,0.263431E+01_r8,0.266800E+01_r8, & 0.270207E+01_r8,0.273656E+01_r8,0.277145E+01_r8,0.280671E+01_r8,0.284248E+01_r8, & 0.287859E+01_r8,0.291516E+01_r8,0.295219E+01_r8,0.298962E+01_r8,0.302746E+01_r8/ DATA (ES_Aux(IES),IES=381,475) / & 0.306579E+01_r8,0.310454E+01_r8,0.314377E+01_r8,0.318351E+01_r8,0.322360E+01_r8, & 0.326427E+01_r8,0.330538E+01_r8,0.334694E+01_r8,0.338894E+01_r8,0.343155E+01_r8, & 0.347456E+01_r8,0.351809E+01_r8,0.356216E+01_r8,0.360673E+01_r8,0.365184E+01_r8, & 0.369744E+01_r8,0.374352E+01_r8,0.379018E+01_r8,0.383743E+01_r8,0.388518E+01_r8, & 0.393344E+01_r8,0.398230E+01_r8,0.403177E+01_r8,0.408175E+01_r8,0.413229E+01_r8, & 0.418343E+01_r8,0.423514E+01_r8,0.428746E+01_r8,0.434034E+01_r8,0.439389E+01_r8, & 0.444808E+01_r8,0.450276E+01_r8,0.455820E+01_r8,0.461423E+01_r8,0.467084E+01_r8, & 0.472816E+01_r8,0.478607E+01_r8,0.484468E+01_r8,0.490393E+01_r8,0.496389E+01_r8, & 0.502446E+01_r8,0.508580E+01_r8,0.514776E+01_r8,0.521047E+01_r8,0.527385E+01_r8, & 0.533798E+01_r8,0.540279E+01_r8,0.546838E+01_r8,0.553466E+01_r8,0.560173E+01_r8, & 0.566949E+01_r8,0.573807E+01_r8,0.580750E+01_r8,0.587749E+01_r8,0.594846E+01_r8, & 0.602017E+01_r8,0.609260E+01_r8,0.616591E+01_r8,0.623995E+01_r8,0.631490E+01_r8, & 0.639061E+01_r8,0.646723E+01_r8,0.654477E+01_r8,0.662293E+01_r8,0.670220E+01_r8, & 0.678227E+01_r8,0.686313E+01_r8,0.694495E+01_r8,0.702777E+01_r8,0.711142E+01_r8, & 0.719592E+01_r8,0.728140E+01_r8,0.736790E+01_r8,0.745527E+01_r8,0.754352E+01_r8, & 0.763298E+01_r8,0.772316E+01_r8,0.781442E+01_r8,0.790676E+01_r8,0.800001E+01_r8, & 0.809435E+01_r8,0.818967E+01_r8,0.828606E+01_r8,0.838343E+01_r8,0.848194E+01_r8, & 0.858144E+01_r8,0.868207E+01_r8,0.878392E+01_r8,0.888673E+01_r8,0.899060E+01_r8, & 0.909567E+01_r8,0.920172E+01_r8,0.930909E+01_r8,0.941765E+01_r8,0.952730E+01_r8/ DATA (ES_Aux(IES),IES=476,570) / & 0.963821E+01_r8,0.975022E+01_r8,0.986352E+01_r8,0.997793E+01_r8,0.100937E+02_r8, & 0.102105E+02_r8,0.103287E+02_r8,0.104481E+02_r8,0.105688E+02_r8,0.106909E+02_r8, & 0.108143E+02_r8,0.109387E+02_r8,0.110647E+02_r8,0.111921E+02_r8,0.113207E+02_r8, & 0.114508E+02_r8,0.115821E+02_r8,0.117149E+02_r8,0.118490E+02_r8,0.119847E+02_r8, & 0.121216E+02_r8,0.122601E+02_r8,0.124002E+02_r8,0.125416E+02_r8,0.126846E+02_r8, & 0.128290E+02_r8,0.129747E+02_r8,0.131224E+02_r8,0.132712E+02_r8,0.134220E+02_r8, & 0.135742E+02_r8,0.137278E+02_r8,0.138831E+02_r8,0.140403E+02_r8,0.141989E+02_r8, & 0.143589E+02_r8,0.145211E+02_r8,0.146845E+02_r8,0.148501E+02_r8,0.150172E+02_r8, & 0.151858E+02_r8,0.153564E+02_r8,0.155288E+02_r8,0.157029E+02_r8,0.158786E+02_r8, & 0.160562E+02_r8,0.162358E+02_r8,0.164174E+02_r8,0.166004E+02_r8,0.167858E+02_r8, & 0.169728E+02_r8,0.171620E+02_r8,0.173528E+02_r8,0.175455E+02_r8,0.177406E+02_r8, & 0.179372E+02_r8,0.181363E+02_r8,0.183372E+02_r8,0.185400E+02_r8,0.187453E+02_r8, & 0.189523E+02_r8,0.191613E+02_r8,0.193728E+02_r8,0.195866E+02_r8,0.198024E+02_r8, & 0.200200E+02_r8,0.202401E+02_r8,0.204626E+02_r8,0.206871E+02_r8,0.209140E+02_r8, & 0.211430E+02_r8,0.213744E+02_r8,0.216085E+02_r8,0.218446E+02_r8,0.220828E+02_r8, & 0.223241E+02_r8,0.225671E+02_r8,0.228132E+02_r8,0.230615E+02_r8,0.233120E+02_r8, & 0.235651E+02_r8,0.238211E+02_r8,0.240794E+02_r8,0.243404E+02_r8,0.246042E+02_r8, & 0.248704E+02_r8,0.251390E+02_r8,0.254109E+02_r8,0.256847E+02_r8,0.259620E+02_r8, & 0.262418E+02_r8,0.265240E+02_r8,0.268092E+02_r8,0.270975E+02_r8,0.273883E+02_r8/ DATA (ES_Aux(IES),IES=571,665) / & 0.276822E+02_r8,0.279792E+02_r8,0.282789E+02_r8,0.285812E+02_r8,0.288867E+02_r8, & 0.291954E+02_r8,0.295075E+02_r8,0.298222E+02_r8,0.301398E+02_r8,0.304606E+02_r8, & 0.307848E+02_r8,0.311119E+02_r8,0.314424E+02_r8,0.317763E+02_r8,0.321133E+02_r8, & 0.324536E+02_r8,0.327971E+02_r8,0.331440E+02_r8,0.334940E+02_r8,0.338475E+02_r8, & 0.342050E+02_r8,0.345654E+02_r8,0.349295E+02_r8,0.352975E+02_r8,0.356687E+02_r8, & 0.360430E+02_r8,0.364221E+02_r8,0.368042E+02_r8,0.371896E+02_r8,0.375790E+02_r8, & 0.379725E+02_r8,0.383692E+02_r8,0.387702E+02_r8,0.391744E+02_r8,0.395839E+02_r8, & 0.399958E+02_r8,0.404118E+02_r8,0.408325E+02_r8,0.412574E+02_r8,0.416858E+02_r8, & 0.421188E+02_r8,0.425551E+02_r8,0.429962E+02_r8,0.434407E+02_r8,0.438910E+02_r8, & 0.443439E+02_r8,0.448024E+02_r8,0.452648E+02_r8,0.457308E+02_r8,0.462018E+02_r8, & 0.466775E+02_r8,0.471582E+02_r8,0.476428E+02_r8,0.481313E+02_r8,0.486249E+02_r8, & 0.491235E+02_r8,0.496272E+02_r8,0.501349E+02_r8,0.506479E+02_r8,0.511652E+02_r8, & 0.516876E+02_r8,0.522142E+02_r8,0.527474E+02_r8,0.532836E+02_r8,0.538266E+02_r8, & 0.543737E+02_r8,0.549254E+02_r8,0.554839E+02_r8,0.560456E+02_r8,0.566142E+02_r8, & 0.571872E+02_r8,0.577662E+02_r8,0.583498E+02_r8,0.589392E+02_r8,0.595347E+02_r8, & 0.601346E+02_r8,0.607410E+02_r8,0.613519E+02_r8,0.619689E+02_r8,0.625922E+02_r8, & 0.632204E+02_r8,0.638550E+02_r8,0.644959E+02_r8,0.651418E+02_r8,0.657942E+02_r8, & 0.664516E+02_r8,0.671158E+02_r8,0.677864E+02_r8,0.684624E+02_r8,0.691451E+02_r8, & 0.698345E+02_r8,0.705293E+02_r8,0.712312E+02_r8,0.719398E+02_r8,0.726542E+02_r8/ DATA (ES_Aux(IES),IES=666,760) / & 0.733754E+02_r8,0.741022E+02_r8,0.748363E+02_r8,0.755777E+02_r8,0.763247E+02_r8, & 0.770791E+02_r8,0.778394E+02_r8,0.786088E+02_r8,0.793824E+02_r8,0.801653E+02_r8, & 0.809542E+02_r8,0.817509E+02_r8,0.825536E+02_r8,0.833643E+02_r8,0.841828E+02_r8, & 0.850076E+02_r8,0.858405E+02_r8,0.866797E+02_r8,0.875289E+02_r8,0.883827E+02_r8, & 0.892467E+02_r8,0.901172E+02_r8,0.909962E+02_r8,0.918818E+02_r8,0.927760E+02_r8, & 0.936790E+02_r8,0.945887E+02_r8,0.955071E+02_r8,0.964346E+02_r8,0.973689E+02_r8, & 0.983123E+02_r8,0.992648E+02_r8,0.100224E+03_r8,0.101193E+03_r8,0.102169E+03_r8, & 0.103155E+03_r8,0.104150E+03_r8,0.105152E+03_r8,0.106164E+03_r8,0.107186E+03_r8, & 0.108217E+03_r8,0.109256E+03_r8,0.110303E+03_r8,0.111362E+03_r8,0.112429E+03_r8, & 0.113503E+03_r8,0.114588E+03_r8,0.115684E+03_r8,0.116789E+03_r8,0.117903E+03_r8, & 0.119028E+03_r8,0.120160E+03_r8,0.121306E+03_r8,0.122460E+03_r8,0.123623E+03_r8, & 0.124796E+03_r8,0.125981E+03_r8,0.127174E+03_r8,0.128381E+03_r8,0.129594E+03_r8, & 0.130822E+03_r8,0.132058E+03_r8,0.133306E+03_r8,0.134563E+03_r8,0.135828E+03_r8, & 0.137109E+03_r8,0.138402E+03_r8,0.139700E+03_r8,0.141017E+03_r8,0.142338E+03_r8, & 0.143676E+03_r8,0.145025E+03_r8,0.146382E+03_r8,0.147753E+03_r8,0.149133E+03_r8, & 0.150529E+03_r8,0.151935E+03_r8,0.153351E+03_r8,0.154783E+03_r8,0.156222E+03_r8, & 0.157678E+03_r8,0.159148E+03_r8,0.160624E+03_r8,0.162117E+03_r8,0.163621E+03_r8, & 0.165142E+03_r8,0.166674E+03_r8,0.168212E+03_r8,0.169772E+03_r8,0.171340E+03_r8, & 0.172921E+03_r8,0.174522E+03_r8,0.176129E+03_r8,0.177755E+03_r8,0.179388E+03_r8/ DATA (ES_Aux(IES),IES=761,855) / & 0.181040E+03_r8,0.182707E+03_r8,0.184382E+03_r8,0.186076E+03_r8,0.187782E+03_r8, & 0.189503E+03_r8,0.191240E+03_r8,0.192989E+03_r8,0.194758E+03_r8,0.196535E+03_r8, & 0.198332E+03_r8,0.200141E+03_r8,0.201963E+03_r8,0.203805E+03_r8,0.205656E+03_r8, & 0.207532E+03_r8,0.209416E+03_r8,0.211317E+03_r8,0.213236E+03_r8,0.215167E+03_r8, & 0.217121E+03_r8,0.219087E+03_r8,0.221067E+03_r8,0.223064E+03_r8,0.225080E+03_r8, & 0.227113E+03_r8,0.229160E+03_r8,0.231221E+03_r8,0.233305E+03_r8,0.235403E+03_r8, & 0.237520E+03_r8,0.239655E+03_r8,0.241805E+03_r8,0.243979E+03_r8,0.246163E+03_r8, & 0.248365E+03_r8,0.250593E+03_r8,0.252830E+03_r8,0.255093E+03_r8,0.257364E+03_r8, & 0.259667E+03_r8,0.261979E+03_r8,0.264312E+03_r8,0.266666E+03_r8,0.269034E+03_r8, & 0.271430E+03_r8,0.273841E+03_r8,0.276268E+03_r8,0.278722E+03_r8,0.281185E+03_r8, & 0.283677E+03_r8,0.286190E+03_r8,0.288714E+03_r8,0.291266E+03_r8,0.293834E+03_r8, & 0.296431E+03_r8,0.299045E+03_r8,0.301676E+03_r8,0.304329E+03_r8,0.307006E+03_r8, & 0.309706E+03_r8,0.312423E+03_r8,0.315165E+03_r8,0.317930E+03_r8,0.320705E+03_r8, & 0.323519E+03_r8,0.326350E+03_r8,0.329199E+03_r8,0.332073E+03_r8,0.334973E+03_r8, & 0.337897E+03_r8,0.340839E+03_r8,0.343800E+03_r8,0.346794E+03_r8,0.349806E+03_r8, & 0.352845E+03_r8,0.355918E+03_r8,0.358994E+03_r8,0.362112E+03_r8,0.365242E+03_r8, & 0.368407E+03_r8,0.371599E+03_r8,0.374802E+03_r8,0.378042E+03_r8,0.381293E+03_r8, & 0.384588E+03_r8,0.387904E+03_r8,0.391239E+03_r8,0.394604E+03_r8,0.397988E+03_r8, & 0.401411E+03_r8,0.404862E+03_r8,0.408326E+03_r8,0.411829E+03_r8,0.415352E+03_r8/ DATA (ES_Aux(IES),IES=856,950) / & 0.418906E+03_r8,0.422490E+03_r8,0.426095E+03_r8,0.429740E+03_r8,0.433398E+03_r8, & 0.437097E+03_r8,0.440827E+03_r8,0.444570E+03_r8,0.448354E+03_r8,0.452160E+03_r8, & 0.455999E+03_r8,0.459870E+03_r8,0.463765E+03_r8,0.467702E+03_r8,0.471652E+03_r8, & 0.475646E+03_r8,0.479674E+03_r8,0.483715E+03_r8,0.487811E+03_r8,0.491911E+03_r8, & 0.496065E+03_r8,0.500244E+03_r8,0.504448E+03_r8,0.508698E+03_r8,0.512961E+03_r8, & 0.517282E+03_r8,0.521617E+03_r8,0.525989E+03_r8,0.530397E+03_r8,0.534831E+03_r8, & 0.539313E+03_r8,0.543821E+03_r8,0.548355E+03_r8,0.552938E+03_r8,0.557549E+03_r8, & 0.562197E+03_r8,0.566884E+03_r8,0.571598E+03_r8,0.576351E+03_r8,0.581131E+03_r8, & 0.585963E+03_r8,0.590835E+03_r8,0.595722E+03_r8,0.600663E+03_r8,0.605631E+03_r8, & 0.610641E+03_r8,0.615151E+03_r8,0.619625E+03_r8,0.624140E+03_r8,0.628671E+03_r8, & 0.633243E+03_r8,0.637845E+03_r8,0.642465E+03_r8,0.647126E+03_r8,0.651806E+03_r8, & 0.656527E+03_r8,0.661279E+03_r8,0.666049E+03_r8,0.670861E+03_r8,0.675692E+03_r8, & 0.680566E+03_r8,0.685471E+03_r8,0.690396E+03_r8,0.695363E+03_r8,0.700350E+03_r8, & 0.705381E+03_r8,0.710444E+03_r8,0.715527E+03_r8,0.720654E+03_r8,0.725801E+03_r8, & 0.730994E+03_r8,0.736219E+03_r8,0.741465E+03_r8,0.746756E+03_r8,0.752068E+03_r8, & 0.757426E+03_r8,0.762819E+03_r8,0.768231E+03_r8,0.773692E+03_r8,0.779172E+03_r8, & 0.784701E+03_r8,0.790265E+03_r8,0.795849E+03_r8,0.801483E+03_r8,0.807137E+03_r8, & 0.812842E+03_r8,0.818582E+03_r8,0.824343E+03_r8,0.830153E+03_r8,0.835987E+03_r8, & 0.841871E+03_r8,0.847791E+03_r8,0.853733E+03_r8,0.859727E+03_r8,0.865743E+03_r8/ DATA (ES_Aux(IES),IES=951,1045) / & 0.871812E+03_r8,0.877918E+03_r8,0.884046E+03_r8,0.890228E+03_r8,0.896433E+03_r8, & 0.902690E+03_r8,0.908987E+03_r8,0.915307E+03_r8,0.921681E+03_r8,0.928078E+03_r8, & 0.934531E+03_r8,0.941023E+03_r8,0.947539E+03_r8,0.954112E+03_r8,0.960708E+03_r8, & 0.967361E+03_r8,0.974053E+03_r8,0.980771E+03_r8,0.987545E+03_r8,0.994345E+03_r8, & 0.100120E+04_r8,0.100810E+04_r8,0.101502E+04_r8,0.102201E+04_r8,0.102902E+04_r8, & 0.103608E+04_r8,0.104320E+04_r8,0.105033E+04_r8,0.105753E+04_r8,0.106475E+04_r8, & 0.107204E+04_r8,0.107936E+04_r8,0.108672E+04_r8,0.109414E+04_r8,0.110158E+04_r8, & 0.110908E+04_r8,0.111663E+04_r8,0.112421E+04_r8,0.113185E+04_r8,0.113952E+04_r8, & 0.114725E+04_r8,0.115503E+04_r8,0.116284E+04_r8,0.117071E+04_r8,0.117861E+04_r8, & 0.118658E+04_r8,0.119459E+04_r8,0.120264E+04_r8,0.121074E+04_r8,0.121888E+04_r8, & 0.122709E+04_r8,0.123534E+04_r8,0.124362E+04_r8,0.125198E+04_r8,0.126036E+04_r8, & 0.126881E+04_r8,0.127731E+04_r8,0.128584E+04_r8,0.129444E+04_r8,0.130307E+04_r8, & 0.131177E+04_r8,0.132053E+04_r8,0.132931E+04_r8,0.133817E+04_r8,0.134705E+04_r8, & 0.135602E+04_r8,0.136503E+04_r8,0.137407E+04_r8,0.138319E+04_r8,0.139234E+04_r8, & 0.140156E+04_r8,0.141084E+04_r8,0.142015E+04_r8,0.142954E+04_r8,0.143896E+04_r8, & 0.144845E+04_r8,0.145800E+04_r8,0.146759E+04_r8,0.147725E+04_r8,0.148694E+04_r8, & 0.149672E+04_r8,0.150655E+04_r8,0.151641E+04_r8,0.152635E+04_r8,0.153633E+04_r8, & 0.154639E+04_r8,0.155650E+04_r8,0.156665E+04_r8,0.157688E+04_r8,0.158715E+04_r8, & 0.159750E+04_r8,0.160791E+04_r8,0.161836E+04_r8,0.162888E+04_r8,0.163945E+04_r8/ DATA (ES_Aux(IES),IES=1046,1140) / & 0.165010E+04_r8,0.166081E+04_r8,0.167155E+04_r8,0.168238E+04_r8,0.169325E+04_r8, & 0.170420E+04_r8,0.171522E+04_r8,0.172627E+04_r8,0.173741E+04_r8,0.174859E+04_r8, & 0.175986E+04_r8,0.177119E+04_r8,0.178256E+04_r8,0.179402E+04_r8,0.180552E+04_r8, & 0.181711E+04_r8,0.182877E+04_r8,0.184046E+04_r8,0.185224E+04_r8,0.186407E+04_r8, & 0.187599E+04_r8,0.188797E+04_r8,0.190000E+04_r8,0.191212E+04_r8,0.192428E+04_r8, & 0.193653E+04_r8,0.194886E+04_r8,0.196122E+04_r8,0.197368E+04_r8,0.198618E+04_r8, & 0.199878E+04_r8,0.201145E+04_r8,0.202416E+04_r8,0.203698E+04_r8,0.204983E+04_r8, & 0.206278E+04_r8,0.207580E+04_r8,0.208887E+04_r8,0.210204E+04_r8,0.211525E+04_r8, & 0.212856E+04_r8,0.214195E+04_r8,0.215538E+04_r8,0.216892E+04_r8,0.218249E+04_r8, & 0.219618E+04_r8,0.220994E+04_r8,0.222375E+04_r8,0.223766E+04_r8,0.225161E+04_r8, & 0.226567E+04_r8,0.227981E+04_r8,0.229399E+04_r8,0.230829E+04_r8,0.232263E+04_r8, & 0.233708E+04_r8,0.235161E+04_r8,0.236618E+04_r8,0.238087E+04_r8,0.239560E+04_r8, & 0.241044E+04_r8,0.242538E+04_r8,0.244035E+04_r8,0.245544E+04_r8,0.247057E+04_r8, & 0.248583E+04_r8,0.250116E+04_r8,0.251654E+04_r8,0.253204E+04_r8,0.254759E+04_r8, & 0.256325E+04_r8,0.257901E+04_r8,0.259480E+04_r8,0.261073E+04_r8,0.262670E+04_r8, & 0.264279E+04_r8,0.265896E+04_r8,0.267519E+04_r8,0.269154E+04_r8,0.270794E+04_r8, & 0.272447E+04_r8,0.274108E+04_r8,0.275774E+04_r8,0.277453E+04_r8,0.279137E+04_r8, & 0.280834E+04_r8,0.282540E+04_r8,0.284251E+04_r8,0.285975E+04_r8,0.287704E+04_r8, & 0.289446E+04_r8,0.291198E+04_r8,0.292954E+04_r8,0.294725E+04_r8,0.296499E+04_r8/ DATA (ES_Aux(IES),IES=1141,1235) / & 0.298288E+04_r8,0.300087E+04_r8,0.301890E+04_r8,0.303707E+04_r8,0.305529E+04_r8, & 0.307365E+04_r8,0.309211E+04_r8,0.311062E+04_r8,0.312927E+04_r8,0.314798E+04_r8, & 0.316682E+04_r8,0.318577E+04_r8,0.320477E+04_r8,0.322391E+04_r8,0.324310E+04_r8, & 0.326245E+04_r8,0.328189E+04_r8,0.330138E+04_r8,0.332103E+04_r8,0.334073E+04_r8, & 0.336058E+04_r8,0.338053E+04_r8,0.340054E+04_r8,0.342069E+04_r8,0.344090E+04_r8, & 0.346127E+04_r8,0.348174E+04_r8,0.350227E+04_r8,0.352295E+04_r8,0.354369E+04_r8, & 0.356458E+04_r8,0.358559E+04_r8,0.360664E+04_r8,0.362787E+04_r8,0.364914E+04_r8, & 0.367058E+04_r8,0.369212E+04_r8,0.371373E+04_r8,0.373548E+04_r8,0.375731E+04_r8, & 0.377929E+04_r8,0.380139E+04_r8,0.382355E+04_r8,0.384588E+04_r8,0.386826E+04_r8, & 0.389081E+04_r8,0.391348E+04_r8,0.393620E+04_r8,0.395910E+04_r8,0.398205E+04_r8, & 0.400518E+04_r8,0.402843E+04_r8,0.405173E+04_r8,0.407520E+04_r8,0.409875E+04_r8, & 0.412246E+04_r8,0.414630E+04_r8,0.417019E+04_r8,0.419427E+04_r8,0.421840E+04_r8, & 0.424272E+04_r8,0.426715E+04_r8,0.429165E+04_r8,0.431634E+04_r8,0.434108E+04_r8, & 0.436602E+04_r8,0.439107E+04_r8,0.441618E+04_r8,0.444149E+04_r8,0.446685E+04_r8, & 0.449241E+04_r8,0.451810E+04_r8,0.454385E+04_r8,0.456977E+04_r8,0.459578E+04_r8, & 0.462197E+04_r8,0.464830E+04_r8,0.467468E+04_r8,0.470127E+04_r8,0.472792E+04_r8, & 0.475477E+04_r8,0.478175E+04_r8,0.480880E+04_r8,0.483605E+04_r8,0.486336E+04_r8, & 0.489087E+04_r8,0.491853E+04_r8,0.494623E+04_r8,0.497415E+04_r8,0.500215E+04_r8, & 0.503034E+04_r8,0.505867E+04_r8,0.508707E+04_r8,0.511568E+04_r8,0.514436E+04_r8/ DATA (ES_Aux(IES),IES=1236,1330) / & 0.517325E+04_r8,0.520227E+04_r8,0.523137E+04_r8,0.526068E+04_r8,0.529005E+04_r8, & 0.531965E+04_r8,0.534939E+04_r8,0.537921E+04_r8,0.540923E+04_r8,0.543932E+04_r8, & 0.546965E+04_r8,0.550011E+04_r8,0.553064E+04_r8,0.556139E+04_r8,0.559223E+04_r8, & 0.562329E+04_r8,0.565449E+04_r8,0.568577E+04_r8,0.571727E+04_r8,0.574884E+04_r8, & 0.578064E+04_r8,0.581261E+04_r8,0.584464E+04_r8,0.587692E+04_r8,0.590924E+04_r8, & 0.594182E+04_r8,0.597455E+04_r8,0.600736E+04_r8,0.604039E+04_r8,0.607350E+04_r8, & 0.610685E+04_r8,0.614036E+04_r8,0.617394E+04_r8,0.620777E+04_r8,0.624169E+04_r8, & 0.627584E+04_r8,0.631014E+04_r8,0.634454E+04_r8,0.637918E+04_r8,0.641390E+04_r8, & 0.644887E+04_r8,0.648400E+04_r8,0.651919E+04_r8,0.655467E+04_r8,0.659021E+04_r8, & 0.662599E+04_r8,0.666197E+04_r8,0.669800E+04_r8,0.673429E+04_r8,0.677069E+04_r8, & 0.680735E+04_r8,0.684415E+04_r8,0.688104E+04_r8,0.691819E+04_r8,0.695543E+04_r8, & 0.699292E+04_r8,0.703061E+04_r8,0.706837E+04_r8,0.710639E+04_r8,0.714451E+04_r8, & 0.718289E+04_r8,0.722143E+04_r8,0.726009E+04_r8,0.729903E+04_r8,0.733802E+04_r8, & 0.737729E+04_r8,0.741676E+04_r8,0.745631E+04_r8,0.749612E+04_r8,0.753602E+04_r8, & 0.757622E+04_r8,0.761659E+04_r8,0.765705E+04_r8,0.769780E+04_r8,0.773863E+04_r8, & 0.777975E+04_r8,0.782106E+04_r8,0.786246E+04_r8,0.790412E+04_r8,0.794593E+04_r8, & 0.798802E+04_r8,0.803028E+04_r8,0.807259E+04_r8,0.811525E+04_r8,0.815798E+04_r8, & 0.820102E+04_r8,0.824427E+04_r8,0.828757E+04_r8,0.833120E+04_r8,0.837493E+04_r8, & 0.841895E+04_r8,0.846313E+04_r8,0.850744E+04_r8,0.855208E+04_r8,0.859678E+04_r8/ DATA (ES_Aux(IES),IES=1331,1425) / & 0.864179E+04_r8,0.868705E+04_r8,0.873237E+04_r8,0.877800E+04_r8,0.882374E+04_r8, & 0.886979E+04_r8,0.891603E+04_r8,0.896237E+04_r8,0.900904E+04_r8,0.905579E+04_r8, & 0.910288E+04_r8,0.915018E+04_r8,0.919758E+04_r8,0.924529E+04_r8,0.929310E+04_r8, & 0.934122E+04_r8,0.938959E+04_r8,0.943804E+04_r8,0.948687E+04_r8,0.953575E+04_r8, & 0.958494E+04_r8,0.963442E+04_r8,0.968395E+04_r8,0.973384E+04_r8,0.978383E+04_r8, & 0.983412E+04_r8,0.988468E+04_r8,0.993534E+04_r8,0.998630E+04_r8,0.100374E+05_r8, & 0.100888E+05_r8,0.101406E+05_r8,0.101923E+05_r8,0.102444E+05_r8,0.102966E+05_r8, & 0.103492E+05_r8,0.104020E+05_r8,0.104550E+05_r8,0.105082E+05_r8,0.105616E+05_r8, & 0.106153E+05_r8,0.106693E+05_r8,0.107234E+05_r8,0.107779E+05_r8,0.108325E+05_r8, & 0.108874E+05_r8,0.109425E+05_r8,0.109978E+05_r8,0.110535E+05_r8,0.111092E+05_r8, & 0.111653E+05_r8,0.112217E+05_r8,0.112782E+05_r8,0.113350E+05_r8,0.113920E+05_r8, & 0.114493E+05_r8,0.115070E+05_r8,0.115646E+05_r8,0.116228E+05_r8,0.116809E+05_r8, & 0.117396E+05_r8,0.117984E+05_r8,0.118574E+05_r8,0.119167E+05_r8,0.119762E+05_r8, & 0.120360E+05_r8,0.120962E+05_r8,0.121564E+05_r8,0.122170E+05_r8,0.122778E+05_r8, & 0.123389E+05_r8,0.124004E+05_r8,0.124619E+05_r8,0.125238E+05_r8,0.125859E+05_r8, & 0.126484E+05_r8,0.127111E+05_r8,0.127739E+05_r8,0.128372E+05_r8,0.129006E+05_r8, & 0.129644E+05_r8,0.130285E+05_r8,0.130927E+05_r8,0.131573E+05_r8,0.132220E+05_r8, & 0.132872E+05_r8,0.133526E+05_r8,0.134182E+05_r8,0.134842E+05_r8,0.135503E+05_r8, & 0.136168E+05_r8,0.136836E+05_r8,0.137505E+05_r8,0.138180E+05_r8,0.138854E+05_r8/ DATA (ES_Aux(IES),IES=1426,1520) / & 0.139534E+05_r8,0.140216E+05_r8,0.140900E+05_r8,0.141588E+05_r8,0.142277E+05_r8, & 0.142971E+05_r8,0.143668E+05_r8,0.144366E+05_r8,0.145069E+05_r8,0.145773E+05_r8, & 0.146481E+05_r8,0.147192E+05_r8,0.147905E+05_r8,0.148622E+05_r8,0.149341E+05_r8, & 0.150064E+05_r8,0.150790E+05_r8,0.151517E+05_r8,0.152250E+05_r8,0.152983E+05_r8, & 0.153721E+05_r8,0.154462E+05_r8,0.155205E+05_r8,0.155952E+05_r8,0.156701E+05_r8, & 0.157454E+05_r8,0.158211E+05_r8,0.158969E+05_r8,0.159732E+05_r8,0.160496E+05_r8, & 0.161265E+05_r8,0.162037E+05_r8,0.162811E+05_r8,0.163589E+05_r8,0.164369E+05_r8, & 0.165154E+05_r8,0.165942E+05_r8,0.166732E+05_r8,0.167526E+05_r8,0.168322E+05_r8, & 0.169123E+05_r8,0.169927E+05_r8,0.170733E+05_r8,0.171543E+05_r8,0.172356E+05_r8, & 0.173173E+05_r8,0.173993E+05_r8,0.174815E+05_r8,0.175643E+05_r8,0.176471E+05_r8, & 0.177305E+05_r8,0.178143E+05_r8,0.178981E+05_r8,0.179826E+05_r8,0.180671E+05_r8, & 0.181522E+05_r8,0.182377E+05_r8,0.183232E+05_r8,0.184093E+05_r8,0.184955E+05_r8, & 0.185823E+05_r8,0.186695E+05_r8,0.187568E+05_r8,0.188447E+05_r8,0.189326E+05_r8, & 0.190212E+05_r8,0.191101E+05_r8,0.191991E+05_r8,0.192887E+05_r8,0.193785E+05_r8, & 0.194688E+05_r8,0.195595E+05_r8,0.196503E+05_r8,0.197417E+05_r8,0.198332E+05_r8, & 0.199253E+05_r8,0.200178E+05_r8,0.201105E+05_r8,0.202036E+05_r8,0.202971E+05_r8, & 0.203910E+05_r8,0.204853E+05_r8,0.205798E+05_r8,0.206749E+05_r8,0.207701E+05_r8, & 0.208659E+05_r8,0.209621E+05_r8,0.210584E+05_r8,0.211554E+05_r8,0.212524E+05_r8, & 0.213501E+05_r8,0.214482E+05_r8,0.215465E+05_r8,0.216452E+05_r8,0.217442E+05_r8/ DATA (ES_Aux(IES),IES=1521,1552) / & 0.218439E+05_r8,0.219439E+05_r8,0.220440E+05_r8,0.221449E+05_r8,0.222457E+05_r8, & 0.223473E+05_r8,0.224494E+05_r8,0.225514E+05_r8,0.226542E+05_r8,0.227571E+05_r8, & 0.228606E+05_r8,0.229646E+05_r8,0.230687E+05_r8,0.231734E+05_r8,0.232783E+05_r8, & 0.233839E+05_r8,0.234898E+05_r8,0.235960E+05_r8,0.237027E+05_r8,0.238097E+05_r8, & 0.239173E+05_r8,0.240254E+05_r8,0.241335E+05_r8,0.242424E+05_r8,0.243514E+05_r8, & 0.244611E+05_r8,0.245712E+05_r8,0.246814E+05_r8,0.247923E+05_r8,0.249034E+05_r8, & 0.250152E+05_r8,0.250152E+05_r8/ ! ES=ES_Aux END SUBROUTINE qsat_data ! ====================================================================== ! *****************************COPYRIGHT******************************* ! (C) Crown copyright Met Office. All rights reserved. ! For further details please refer to the file COPYRIGHT.txt ! which you should have received as part of this distribution. ! *****************************COPYRIGHT******************************* ! !+ Saturation Specific Humidity/mixing Scheme(Qsat):Vapour to Liquid/Ice ! ! Subroutine Interface: SUBROUTINE QSAT_mix ( & ! Output field & QmixS & ! Input fields &, T, P & ! Array dimensions &, NPNTS & ! logical control &, lq_mix & &, mskant & & ) IMPLICIT NONE ! ! Purpose: ! Returns a saturation specific humidity or mixing ratio given a ! temperature and pressure using the saturation vapour pressure ! calculated using the Goff-Gratch formulae, adopted by the WMO as ! taken from Landolt-Bornstein, 1987 Numerical Data and Functional ! Relationships in Science and Technolgy. Group V/vol 4B meteorology. ! Phyiscal and Chemical properties or air, P35. ! ! Values in the lookup table are over water above 0 deg C and over ice ! below this temperature. ! ! Method: ! Uses lookup tables to find eSAT, calculates qSAT directly from that. ! ! Current Owner of Code: OWNER OF LARGE_SCALE CLOUD CODE. ! ! Code Description: ! Language: FORTRAN 77 + common extensions also in Fortran90. ! This code is written to UMDP3 version 6 programming standards. ! ! Documentation: UMDP No.29 ! ! Declarations: ! ! Global Variables:---------------------------------------------------- !----------------------------------------------------------------------- ! Subroutine Arguments !----------------------------------------------------------------------- ! ! arguments with intent in. ie: input variables. ! INTEGER, INTENT(in) :: npnts ! Points (=horizontal dimensions) being processed by qSAT scheme. REAL(KIND=r8), INTENT(in) :: T(npnts) ! Temperature (K). REAL(KIND=r8), INTENT(in) :: P(npnts) ! Pressure (Pa). INTEGER(KIND=i8), INTENT(in) :: mskant(npnts) LOGICAL, INTENT(in) :: lq_mix ! .true. return qsat as a mixing ratio ! .false. return qsat as a specific humidity ! ! arguments with intent out ! REAL(KIND=r8), INTENT(out) :: QmixS(npnts) ! Output Saturation mixing ratio or saturation ! specific humidity at temperature T and pressure ! P (kg/kg). !----------------------------------------------------------------------- ! Local scalars !----------------------------------------------------------------------- INTEGER :: itable, itable_p1 ! Work variables REAL(KIND=r8) :: atable, atable_p1 ! Work variables REAL(KIND=r8) :: fsubw, fsubw_p1 ! FACTOR THAT CONVERTS FROM SAT VAPOUR PRESSURE IN A PURE ! WATER SYSTEM TO SAT VAPOUR PRESSURE IN AIR. REAL(KIND=r8) :: TT, TT_P1 REAL(KIND=r8), PARAMETER :: R_delta_T = 1.0_r8/delta_T INTEGER :: I, II ! !----------------------------------------------------------------------- !----------------------------------------------------------------------- ! ! loop over points ! DO i = 1, npnts-1,2 IF(mskant(i) == 1_i8)THEN ! Compute the factor that converts from sat vapour pressure in a ! pure water system to sat vapour pressure in air, fsubw. ! This formula is taken from equation A4.7 of Adrian Gill's book: ! Atmosphere-Ocean Dynamics. Note that his formula works in terms ! of pressure in mb and temperature in Celsius, so conversion of ! units leads to the slightly different equation used here. fsubw = 1.0_r8 + 1.0E-8_r8 * P(I) * ( 4.5_r8 + & 6.0E-4_r8*( T(I) - zerodegC ) * ( T(I) - zerodegC ) ) fsubw_p1 = 1.0_r8 + 1.0E-8_r8 * P(I+1) * ( 4.5_r8 + & 6.0E-4_r8*( T(I+1) - zerodegC ) * ( T(I+1) - zerodegC ) ) ! Use the lookup table to find saturated vapour pressure, and store ! it in qmixs. ! TT = MAX(T_low,T(I)) TT = MIN(T_high,TT) atable = (TT - T_low + delta_T) * R_delta_T itable = atable atable = atable - itable TT_p1 = MAX(T_LOW,T(I+1)) TT_p1 = MIN(T_HIGH,TT_p1) ATABLE_p1 = (TT_p1 - T_LOW + DELTA_T) * R_DELTA_T ITABLE_p1 = ATABLE_p1 ATABLE_p1 = ATABLE_p1 - ITABLE_p1 QmixS(I) = (1.0_r8 - atable) *ES(itable) + atable*ES(itable+1) QmixS(I+1) = (1.0_r8 - atable_p1)*ES(itable_p1) + & atable_p1*ES(itable_p1+1) ! Multiply by fsubw to convert to saturated vapour pressure in air ! (equation A4.6 of Adrian Gill's book). ! QmixS(I) = QmixS(I) * fsubw QmixS(I+1) = QmixS(I+1) * fsubw_p1 ! ! Now form the accurate expression for qmixs, which is a rearranged ! version of equation A4.3 of Gill's book. ! !----------------------------------------------------------------------- ! For mixing ratio, rsat = epsilon *e/(p-e) ! e - saturation vapour pressure ! Note applying the fix to qsat for specific humidity at low pressures ! is not possible, this implies mixing ratio qsat tends to infinity. ! If the pressure is very low then the mixing ratio value may become ! very large. !----------------------------------------------------------------------- IF (lq_mix) THEN QmixS(I) = ( epsilonk*QmixS(I) ) / & ( MAX(P(I), 1.1_r8*QmixS(I)) - QmixS(I) ) QmixS(I+1) = ( epsilonk*QmixS(I+1) ) / & ( MAX(P(I+1),1.1_r8*QmixS(I+1)) - QmixS(I+1) ) !----------------------------------------------------------------------- ! For specific humidity, qsat = epsilon*e/(p-(1-epsilon)e) ! ! Note that at very low pressure we apply a fix, to prevent a ! singularity (qsat tends to 1. kg/kg). !----------------------------------------------------------------------- ELSE QmixS(I) = ( epsilonk*QmixS(I) ) / & ( MAX(P(I), QmixS(I)) - one_minus_epsilon*QmixS(I) ) QmixS(I+1) = ( epsilonk*QmixS(I+1) ) / & ( MAX(P(I+1),QmixS(I+1)) - one_minus_epsilon*QmixS(I+1)) ENDIF ! test on lq_mix ! END IF END DO ! Npnts_do_1 II = I DO I = II, NPNTS IF(mskant(i) == 1_i8)THEN fsubw = 1.0_r8 + 1.0E-8_r8*P(I)*( 4.5_r8 + & & 6.0E-4_r8*( T(I) - zerodegC )*( T(I) - zerodegC ) ) ! TT = MAX(T_low,T(I)) TT = MIN(T_high,TT) atable = (TT - T_low + delta_T) * R_delta_T itable = atable atable = atable - itable QmixS(I) = (1.0_r8 - atable)*ES(itable) + atable*ES(itable+1) QmixS(I) = QmixS(I) * fsubw IF (lq_mix) THEN QmixS(I) = ( epsilonk*QmixS(I) ) / & & ( MAX(P(I),1.1_r8*QmixS(i)) - QmixS(I) ) ELSE QmixS(I) = ( epsilonk*QmixS(I) ) / & & ( MAX(P(I),QmixS(I)) - one_minus_epsilon*QmixS(I) ) ENDIF ! test on lq_mix END IF END DO ! Npnts_do_1 RETURN END SUBROUTINE QSAT_mix ! ====================================================================== END MODULE OceanModel !PROGRAM Main ! USE OceanModel ! IMPLICIT NONE ! !END PROGRAM Main