Back to home page

MITgcm
 
 

    

File indexing completed on 2018-03-02 18:38:13 UTC

view on githubraw file Latest commit add29e06 on 2018-01-31 20:35:05 UTC
6d54cf9ca1 Ed H*0001 #include "BULK_FORCE_OPTIONS.h"
6a1d3c464b Jean*0002 
69f66dfb04 Jean*0003 CBOP
                0004 C     !ROUTINE: BULKF_FORMULA_LANL
                0005 C     !INTERFACE:
                0006       SUBROUTINE BULKF_FORMULA_LANL(
7753507405 Curt*0007      I                           uw, vw, us, Ta, Qa, nc, tsf_in,
6a1d3c464b Jean*0008      O                           flwupa, flha, fsha, df0dT,
69f66dfb04 Jean*0009      O                           ust, vst, evp, ssq, dEvdT,
679d149d01 Jean*0010      I                           iceornot, myThid )
6a1d3c464b Jean*0011 
69f66dfb04 Jean*0012 C     !DESCRIPTION: \bv
                0013 C     *==========================================================*
                0014 C     | SUBROUTINE BULKF_FORMULA_LANL
                0015 c     | o Calculate bulk formula fluxes over open ocean or seaice
                0016 C     *==========================================================*
                0017 C     \ev
                0018 C swd -- bulkf formula used in bulkf and ice pkgs
                0019 C        taken from exf package
                0020 C
                0021 C     wind stress = (ust,vst) = rhoA * Cd * Ws * (del.u,del.v)
                0022 C     Sensib Heat flux = fsha = rhoA * Ch * Ws * del.T * CpAir
                0023 C     Latent Heat flux = flha = rhoA * Ce * Ws * del.Q * Lvap
                0024 C                      = -Evap * Lvap
                0025 C   with Ws = wind speed = sqrt(del.u^2 +del.v^2) ;
                0026 C        del.T = Tair - Tsurf ; del.Q = Qair - Qsurf
                0027 C        Cd,Ch,Ce = transfer coefficient for momentum, sensible
                0028 C                    & latent heat flux [no units]
                0029 C     *==========================================================*
7753507405 Curt*0030 
69f66dfb04 Jean*0031 C     !USES:
7753507405 Curt*0032        IMPLICIT NONE
69f66dfb04 Jean*0033 C     === Global variables ===
7753507405 Curt*0034 #include "EEPARAMS.h"
                0035 #include "SIZE.h"
                0036 #include "PARAMS.h"
6a1d3c464b Jean*0037 #include "BULKF_PARAMS.h"
7753507405 Curt*0038 
69f66dfb04 Jean*0039 C     !INPUT/OUTPUT PARAMETERS:
                0040 C     input:
                0041       _RL uw                 ! zonal wind speed (at grid center) [m/s]
                0042       _RL vw                 ! meridional wind speed (at grid center) [m/s]
                0043       _RL us                 ! wind speed        [m/s]   at height hu
                0044       _RL Ta                 ! air temperature   [K]     at height ht
                0045       _RL Qa                 ! specific humidity [kg/kg] at heigth ht
7753507405 Curt*0046       _RL nc                 ! fraction cloud cover
69f66dfb04 Jean*0047       _RL tsf_in             ! sea-ice or sea surface temperature [oC]
                0048       INTEGER iceornot       ! 0=open water, 1=sea-ice, 2=sea-ice with snow
                0049       INTEGER myThid         ! my Thread Id number
                0050 C     output:
                0051       _RL flwupa             ! upward long wave radiation (>0 upward) [W/m2]
                0052       _RL flha               ! latent heat flux         (>0 downward) [W/m2]
                0053       _RL fsha               ! sensible heat flux       (>0 downward) [W/m2]
                0054       _RL df0dT              ! derivative of heat flux with respect to Tsf [W/m2/K]
                0055       _RL ust                ! zonal wind stress (at grid center)     [N/m2]
                0056       _RL vst                ! meridional wind stress (at grid center)[N/m2]
548c63e38c Jean*0057       _RL evp                ! evaporation rate (over open water) [kg/m2/s]
69f66dfb04 Jean*0058       _RL ssq                ! surface specific humidity          [kg/kg]
548c63e38c Jean*0059       _RL dEvdT              ! derivative of evap. with respect to Tsf [kg/m2/s/K]
69f66dfb04 Jean*0060 CEOP
6a1d3c464b Jean*0061 
                0062 #ifdef ALLOW_BULK_FORCE
                0063 
69f66dfb04 Jean*0064 C     == Local variables ==
7753507405 Curt*0065       _RL dflhdT             ! derivative of latent heat with respect to T
                0066       _RL dfshdT             ! derivative of sensible heat with respect to T
                0067       _RL dflwupdT           ! derivative of long wave with respect to T
69f66dfb04 Jean*0068 
                0069       _RL tsf                ! surface temperature [K]
                0070       _RL ht                 ! height for air temperature [m]
                0071 c     _RL hq                 ! height for humidity [m]
                0072       _RL hu                 ! height for wind speed [m]
                0073 c     _RL zref               ! reference height [m]
                0074       _RL usm                ! wind speed limited [m/s]
                0075 c     _RL umin               ! minimum wind speed used for drag-coeff [m/s]
                0076       _RL lath               ! latent heat of vaporization or sublimation
                0077       _RL t0                 ! virtual temperature [K]
                0078       _RL deltap             ! potential temperature diff [K]
                0079       _RL delq               ! specific humidity difference [kg/kg]
                0080       _RL ustar              ! friction velocity [m/s]
                0081       _RL tstar              ! temperature scale [K]
                0082       _RL qstar              ! humidity scale  [kg/kg]
                0083       _RL rd                 ! = sqrt(Cd)          [-]
                0084       _RL re                 ! = Ce / sqrt(Cd)     [-]
                0085       _RL rh                 ! = Ch / sqrt(Cd)     [-]
                0086       _RL rdn, ren, rhn      ! initial (neutral) values of rd, re, rh
                0087       _RL stable             ! = 1 if stable ; = 0 if unstable
                0088       _RL huol               ! stability parameter [-]
                0089       _RL x                  ! stability function  [-]
                0090       _RL xsq                ! = x^2               [-]
                0091       _RL psimh              ! momentum stability function
                0092       _RL psixh              ! latent & sensib. stability function
                0093       _RL czol               ! = zref*Karman_cst*gravity
                0094       _RL aln                ! = log(ht/zref)
                0095 c     _RL cdalton            ! coeff to evaluate Dalton Number
a10fbde4a0 Jean*0096 c     _RL mixratio
                0097 c     _RL ea
69f66dfb04 Jean*0098 c     _RL psim_fac
                0099       _RL tau                ! surface stress  coef = ?
                0100       _RL csha               ! sensib.heat flx coef = rhoA * Ws * Ch * CpAir
                0101       _RL clha               ! latent heat flx coef = rhoA * Ws * Ce * Lvap
7753507405 Curt*0102       _RL zice
548c63e38c Jean*0103       _RL ssq0, ssq1, ssq2   ! constant used in surface specific humidity
69f66dfb04 Jean*0104       _RL p0                 ! reference sea-level atmospheric pressure [mb]
679d149d01 Jean*0105       _RL bulkf_Cdn          ! drag coefficient
69f66dfb04 Jean*0106       INTEGER niter_bulk, iter
7753507405 Curt*0107 
69f66dfb04 Jean*0108 C     == external Functions
679d149d01 Jean*0109 c     _RL       exf_BulkCdn
                0110 c     external  exf_BulkCdn
6a1d3c464b Jean*0111 c     _RL       exf_BulkqSat
                0112 c     external  exf_BulkqSat
                0113 c     _RL       exf_BulkRhn
                0114 c     external  exf_BulkRhn
7753507405 Curt*0115 
69f66dfb04 Jean*0116       DATA   ssq0,           ssq1,           ssq2
548c63e38c Jean*0117      &     / 3.797915 _d 0 , 7.93252 _d -6 , 2.166847 _d -3 /
69f66dfb04 Jean*0118       DATA   p0 / 1013. _d 0 /
548c63e38c Jean*0119 
69f66dfb04 Jean*0120 C--- Compute turbulent surface fluxes
6a1d3c464b Jean*0121               ht =  2. _d 0
69f66dfb04 Jean*0122 c             hq =  2. _d 0
6a1d3c464b Jean*0123               hu = 10. _d 0
69f66dfb04 Jean*0124 c             zref = 10. _d 0
6a1d3c464b Jean*0125               zice = 0.0005 _d 0
7753507405 Curt*0126               aln = log(ht/zref)
                0127               niter_bulk = 5
69f66dfb04 Jean*0128 c             cdalton = 0.0346000 _d 0
7753507405 Curt*0129               czol = zref*xkar*gravity
69f66dfb04 Jean*0130 c             psim_fac=5. _d 0
                0131 c             umin=1. _d 0
                0132 
7753507405 Curt*0133               lath=Lvap
                0134               if (iceornot.gt.0) lath=Lvap+Lfresh
                0135               Tsf=Tsf_in+Tf0kel
69f66dfb04 Jean*0136 C-   Wind speed
6a1d3c464b Jean*0137               if (us.eq.0. _d 0) then
7753507405 Curt*0138                 us = sqrt(uw*uw + vw*vw)
                0139               endif
                0140               usm = max(us,umin)
                0141 c
6a1d3c464b Jean*0142               t0     = Ta*(1. _d 0 + humid_fac*Qa)
                0143 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
69f66dfb04 Jean*0144 cQQ           ssq    = 0.622*6.11*exp(22.47*(1.d0-Tf0kel/tsf))/p0
548c63e38c Jean*0145 c             ssq = 3.797915 _d 0*exp(
                0146 c    &                lath*(7.93252 _d -6 - 2.166847 _d -3/Tsf)
                0147 c    &                               )/p0
                0148               ssq = ssq0*exp( lath*(ssq1-ssq2/Tsf) ) / p0
69f66dfb04 Jean*0149 
7753507405 Curt*0150               deltap = ta  - tsf + gamma_blk*ht
                0151               delq   = Qa - ssq
69f66dfb04 Jean*0152 
                0153 C--  initialize estimate exchange coefficients
7753507405 Curt*0154               rdn=xkar/(log(zref/zice))
                0155               rhn=rdn
                0156               ren=rdn
69f66dfb04 Jean*0157 C--  calculate turbulent scales
7753507405 Curt*0158               ustar=rdn*usm
                0159               tstar=rhn*deltap
                0160               qstar=ren*delq
69f66dfb04 Jean*0161 
                0162 C--  iteration with psi-functions to find transfer coefficients
7753507405 Curt*0163               do iter=1,niter_bulk
                0164                  huol   = czol/ustar**2 *(tstar/t0 +
6a1d3c464b Jean*0165      &                    qstar/(1. _d 0/humid_fac+Qa))
                0166                  huol   = sign( min(abs(huol),10. _d 0), huol)
7753507405 Curt*0167                  stable = 5. _d -1 + sign(5. _d -1 , huol)
6a1d3c464b Jean*0168                  xsq = max(sqrt(abs(1. _d 0 - 16. _d 0*huol)),1. _d 0)
7753507405 Curt*0169                  x      = sqrt(xsq)
6a1d3c464b Jean*0170                  psimh = -5. _d 0*huol*stable + (1. _d 0-stable)*
                0171      &                    (2. _d 0*log(5. _d -1*(1. _d 0+x)) +
                0172      &                     2. _d 0*log(5. _d -1*(1. _d 0+xsq)) -
                0173      &                     2. _d 0*atan(x) + pi*.5 _d 0)
                0174                  psixh  = -5. _d 0*huol*stable + (1. _d 0-stable)*
                0175      &                     (2. _d 0*log(5. _d -1*(1. _d 0+xsq)))
69f66dfb04 Jean*0176 
                0177 C--  Update the transfer coefficients
6a1d3c464b Jean*0178                  rd = rdn/(1. _d 0 + rdn*(aln-psimh)/xkar)
                0179                  rh = rhn/(1. _d 0 + rhn*(aln-psixh)/xkar)
7753507405 Curt*0180                  re = rh
69f66dfb04 Jean*0181 C--  Update ustar, tstar, qstar using updated, shifted coefficients.
7753507405 Curt*0182                  ustar = rd*usm
                0183                  qstar = re*delq
                0184                  tstar = rh*deltap
6a1d3c464b Jean*0185               enddo
69f66dfb04 Jean*0186 
                0187               tau   = rhoa*ustar**2
                0188               tau   = tau*us/usm
                0189               csha  = rhoa*cpair*us*rh*rd
                0190               clha  = rhoa*lath*us*re*rd
                0191 
                0192               fsha  = csha*deltap
                0193               flha  = clha*delq
                0194               evp   = -flha/lath
                0195 
                0196 C--  Upward long wave radiation
7753507405 Curt*0197 cQQ           mixratio=Qa/(1-Qa)
                0198 cQQ           ea=p0*mixratio/(0.62197+mixratio)
f4245d1665 Curt*0199 cQQ           flwupa=-0.985*stefan*tsf**4
7753507405 Curt*0200 cQQ  &                  *(0.39-0.05*sqrt(ea))
                0201 cQQ  &                  *(1-0.6*nc**2)
                0202               if (iceornot.eq.0) then
6a1d3c464b Jean*0203                flwupa=ocean_emissivity*stefan*tsf**4
                0204                dflwupdT=4. _d 0*ocean_emissivity*stefan*tsf**3
69f66dfb04 Jean*0205               elseif (iceornot.eq.2) then
6a1d3c464b Jean*0206                 flwupa=snow_emissivity*stefan*tsf**4
                0207                 dflwupdT=4. _d 0*snow_emissivity*stefan*tsf**3
69f66dfb04 Jean*0208               else
6a1d3c464b Jean*0209                 flwupa=ice_emissivity*stefan*tsf**4
                0210                 dflwupdT=4. _d 0*ice_emissivity*stefan*tsf**3
7753507405 Curt*0211               endif
                0212 cQQ           dflhdT = -clha*Tf0kel*ssq*22.47/(tsf**2)
                0213 c             dflhdT = -clha*Lath*ssq/(Rvap*tsf**2)
548c63e38c Jean*0214 c             dflhdT = -clha*ssq*Lath*2.166847 _d -3/(Tsf**2)
                0215               dEvdT  =  clha*ssq*ssq2/(Tsf*Tsf)
                0216               dflhdT = -lath*dEvdT
7753507405 Curt*0217               dfshdT = -csha
6a1d3c464b Jean*0218 cQQ           dflwupdT= 4.*0.985*stefan*tsf**3
7753507405 Curt*0219 cQQ  &                  *(0.39-0.05*sqrt(ea))
                0220 cQQ  &                  *(1-0.6*nc**2)
                0221 c total derivative with respect to surface temperature
6a1d3c464b Jean*0222               df0dT=-dflwupdT+dfshdT+dflhdT
69f66dfb04 Jean*0223 
                0224 C--  Wind stress at center points
                0225 C-   in-lining of function: exf_BulkCdn(umps) = cdrag_1/umps + cdrag_2 + cdrag_3*umps
679d149d01 Jean*0226               bulkf_Cdn = cdrag_1/usm + cdrag_2 + cdrag_3*usm
                0227               ust = rhoa*bulkf_Cdn*us*uw
                0228               vst = rhoa*bulkf_Cdn*us*vw
7753507405 Curt*0229 #endif /*ALLOW_BULK_FORCE*/
                0230 
6a1d3c464b Jean*0231       RETURN
                0232       END
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated from https://github.com/MITgcm/MITgcm by the 2.2.1-MITgcm-0.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!