Back to home page

MITgcm
 
 

    

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

view on githubraw file Latest commit add29e06 on 2018-01-31 20:35:05 UTC
ced0783fba Jean*0001 #include "CHEAPAML_OPTIONS.h"
                0002 #undef ALLOW_THSICE
                0003 
6f955c5285 Jean*0004 CBOP
                0005 C     !ROUTINE: CHEAPAML_LANL_FLUX
                0006 C     !INTERFACE:
                0007       SUBROUTINE CHEAPAML_LANL_FLUX(
                0008      I                    i,j,bi,bj,
                0009      O                    fsha, flha, evp, xolw, ssqt, q100 )
                0010 
                0011 C     !DESCRIPTION:
                0012 C     ==================================================================
                0013 C     SUBROUTINE cheapaml_LANL_flux
                0014 C     ==================================================================
                0015 C     o compute surface fluxes using LANL algorithms
                0016 
                0017 C     !USES:
                0018       IMPLICIT NONE
                0019 C     === Global variables ===
ced0783fba Jean*0020 #include "SIZE.h"
6f955c5285 Jean*0021 #include "EEPARAMS.h"
ced0783fba Jean*0022 #include "PARAMS.h"
                0023 #include "DYNVARS.h"
                0024 #include "GRID.h"
6f955c5285 Jean*0025 c#include "FFIELDS.h"
                0026 c#ifdef ALLOW_THSICE
                0027 c#include "THSICE_VARS.h"
                0028 c#endif
ced0783fba Jean*0029 #include "CHEAPAML.h"
                0030 
6f955c5285 Jean*0031 C     !INPUT PARAMETERS:
                0032       INTEGER i,j,bi,bj
                0033 C     !OUTPUT PARAMETERS:
                0034       _RL fsha, flha, evp, xolw, ssqt, q100
                0035 CEOP
                0036 
                0037 C       Output:
                0038 C       ustress, vstress - wind stress
                0039 C       fsha - sensible heat flux
                0040 C       flha - latent heat flux
                0041 C       xolw - oceanic upwelled long wave radiation
                0042 C       ssqt - sat. specific humidity at atm layer top
                0043 
                0044 C       Input:
                0045 C       uwind, vwind  - mean wind speed (m/s)
                0046 C       Tair  - mean air temperature (K)  at height ht (m)
                0047 C       theta(k=1) - sea surface temperature (C)
                0048 C       Qair - Specific humidity kg/kg
                0049 C       Solar - short wave net solar flux at surface (W/m^2)
                0050 C       Tr - relaxation profile for temperature on boundaries (C)
                0051 C       qr - relaxation profile for specific humidity (kg/kg)
                0052 C       i,j,bi,bj - indices of data
                0053 
                0054 C     !LOCAL VARIABLES:
                0055 C     iceornot :: variables to include seaice effect
                0056       INTEGER iceornot
                0057       _RL deltaTm
                0058       _RL uss,usm,uw,vw
                0059       _RL cheapaml_BulkCdn
                0060       _RL to
                0061       _RL t
                0062       _RL t0,QaR
                0063       _RL ssq, q
                0064       _RL deltap, delq, pt, psx100, z100ol
                0065       _RL rdn,ren,rhn,zice,zref
                0066       _RL rd,re,rh,tta,ttas,toa,ttt
                0067       _RL ustar,tstar,qstar,ht,hu,hq
                0068       _RL aln,cdalton,czol,psim_fac
                0069       _RL huol,stable,xsq,x,psimh,psixh
                0070       _RL clha, csha
                0071       INTEGER niter_bulk,iter
                0072 
                0073 C useful values
                0074 C hardwire atmospheric relative humidity at 80%
ced0783fba Jean*0075         QaR=0.8 _d 0
6f955c5285 Jean*0076 C factor to compute rainfall from specific humidity
                0077 C inverse of time step
ced0783fba Jean*0078         deltaTm=1. _d 0/deltaT
6f955c5285 Jean*0079 C reference values to compute turbulent flux
e900d93663 Jean*0080               ht=zt
                0081               hq=zq
                0082               hu=zu
ced0783fba Jean*0083               zref = zt
e900d93663 Jean*0084               zice=.0005 _d 0
ced0783fba Jean*0085               aln = log(ht/zref)
6f955c5285 Jean*0086 C for iterating on turbulence
ced0783fba Jean*0087               niter_bulk = 5
                0088               cdalton = 0.0346000 _d 0
                0089               czol = zref*xkar*gravity
                0090               psim_fac=5. _d 0
                0091 
6f955c5285 Jean*0092 C     determine wind stress
ced0783fba Jean*0093         IF(.NOT.useStressOption)THEN
                0094 
fac34f1c8c Nico*0095              if (maskC(i,j,1,bi,bj).ne.0. _d 0) then
ced0783fba Jean*0096 #ifdef ALLOW_THSICE
                0097                if (ICEMASK(i,j,bi,bj).gt.0. _d 0) then
                0098                 if (snowheight(i,j,bi,bj).gt.3. _d -1) then
                0099                    iceornot=2
                0100                  else
                0101                    iceornot=1
                0102                  endif
                0103                else
                0104                  iceornot=0
                0105                endif
                0106 #else
                0107                iceornot=0
                0108 #endif
                0109                        uw=uwind(i,j,bi,bj)
                0110                        vw=vwind(i,j,bi,bj)
                0111                        uss=sqrt(uw**2+vw**2)
                0112                        usm=max(uss,1. _d 0)
                0113                   cheapaml_BulkCdn = cdrag_1/usm + cdrag_2 + cdrag_3*usm
                0114                        ustress(i,j,bi,bj)= rhoa*cheapaml_BulkCdn*uss*uw
                0115                        vstress(i,j,bi,bj)= rhoa*cheapaml_BulkCdn*uss*vw
                0116              else
                0117                usm=0. _d 0
                0118                ustress(i,j,bi,bj) = 0. _d 0
                0119                vstress(i,j,bi,bj) = 0. _d 0
4fa4901be6 Nico*0120                 endif
6f955c5285 Jean*0121 C wind stress computed
ced0783fba Jean*0122                 ENDIF
6f955c5285 Jean*0123 C diabatic and freshwater flux forcing
ced0783fba Jean*0124         to=theta(i,j,1,bi,bj)
                0125         t=Tair(i,j,bi,bj)
4fa4901be6 Nico*0126         toa=to+Celsius2K
                0127         tta=t+Celsius2K
e900d93663 Jean*0128         ttas=tta+gamma_blk*zref
2616d73cb2 Nico*0129         ttt=tta-( cheaphgrid(i,j,bi,bj)- zref)*gamma_blk
                0130         pt=p0*(1-gamma_blk*cheaphgrid(i,j,bi,bj)/ttas)
6f955c5285 Jean*0131      &     **(gravity/gamma_blk/gasR)
ced0783fba Jean*0132 
6f955c5285 Jean*0133 C specific humidities
2616d73cb2 Nico*0134               ssq= ssq0*exp( lath*(ssq1-ssq2/toa) ) / p0
ced0783fba Jean*0135               ssqt = ssq0*exp( lath*(ssq1-ssq2/ttt) ) / pt
6f955c5285 Jean*0136 C     saturation no more at the top:
                0137               ssqt = 0.7 _d 0*ssq
ced0783fba Jean*0138 
6f955c5285 Jean*0139             if (useFreshWaterFlux) then
ced0783fba Jean*0140             q=qair(i,j,bi,bj)
                0141             else
                0142             q=QaR * ssq
                0143             endif
6f955c5285 Jean*0144 
                0145 C adjust temperature from reference height to formula height
ced0783fba Jean*0146             deltap = t  - to + gamma_blk*(zref-ht)
                0147             delq   = q - ssq
                0148             ttas   = tta+gamma_blk*(zref-ht)
                0149             t0     = ttas*(1. _d 0 + humid_fac*q)
                0150 
6f955c5285 Jean*0151 C initialize estimate exchange coefficients
ced0783fba Jean*0152               rdn=xkar/(log(zref/zice))
                0153               rhn=rdn
                0154               ren=rdn
6f955c5285 Jean*0155 C calculate turbulent scales
ced0783fba Jean*0156               ustar=rdn*usm
                0157               tstar=rhn*deltap
                0158               qstar=ren*delq
6f955c5285 Jean*0159 
                0160 C iteration with psi-functions to find transfer coefficients
ced0783fba Jean*0161               do iter=1,niter_bulk
                0162                  huol   = czol/ustar**2 *(tstar/t0 +
                0163      &                    qstar/(1. _d 0/humid_fac+q))
                0164                  huol   = sign( min(abs(huol),10. _d 0), huol)
                0165                  stable = 5. _d -1 + sign(5. _d -1 , huol)
                0166                  xsq = max(sqrt(abs(1. _d 0 - 16. _d 0*huol)),1. _d 0)
                0167                  x      = sqrt(xsq)
                0168                  psimh = -5. _d 0*huol*stable + (1. _d 0-stable)*
                0169      &                    (2. _d 0*log(5. _d -1*(1. _d 0+x)) +
                0170      &                     2. _d 0*log(5. _d -1*(1. _d 0+xsq)) -
                0171      &                     2. _d 0*atan(x) + pi*.5 _d 0)
                0172                  psixh  = -5. _d 0*huol*stable + (1. _d 0-stable)*
                0173      &                     (2. _d 0*log(5. _d -1*(1. _d 0+xsq)))
                0174 
6f955c5285 Jean*0175 C Update the transfer coefficients
ced0783fba Jean*0176 
                0177                  rd = rdn/(1. _d 0 + rdn*(aln-psimh)/xkar)
                0178                  rh = rhn/(1. _d 0 + rhn*(aln-psixh)/xkar)
                0179                  re = rh
6f955c5285 Jean*0180 C  Update ustar, tstar, qstar using updated, shifted coefficients.
ced0783fba Jean*0181                  ustar = rd*usm
                0182                  qstar = re*delq
                0183                  tstar = rh*deltap
                0184               enddo
6f955c5285 Jean*0185 
ced0783fba Jean*0186                         usm=max(uss,0.5 _d 0)
                0187                 csha   = rhoa*cpair*usm*rh*rd
                0188                 clha   = rhoa*lath*usm*re*rd
6f955c5285 Jean*0189 
ced0783fba Jean*0190                 fsha  = csha*deltap
                0191                 flha  = clha*delq
                0192                 evp   = -flha/lath
                0193 
6f955c5285 Jean*0194 C the sensible and latent heat fluxes, fsha and flha,
                0195 C are computed so that positive values are downward.
                0196 C the convention for cheapaml is upward fluxes are positive,
                0197 C so they must be multiplied by -1
ced0783fba Jean*0198         fsha=-fsha
                0199         flha=-flha
                0200 
6f955c5285 Jean*0201 C oceanic upwelled long wave
ced0783fba Jean*0202         xolw=stefan*(toa)**4
6f955c5285 Jean*0203 C compute specific humidity at 100m
ced0783fba Jean*0204                  huol   = czol/ustar**2 *(tstar/t0 +
                0205      &                    qstar/(1. _d 0/humid_fac+q))
                0206                  huol   = sign( min(abs(huol),10. _d 0), huol)
                0207                  stable = 5. _d -1 + sign(5. _d -1 , huol)
e900d93663 Jean*0208                  z100ol   = 100. _d 0 *xkar*gravity/ustar**2 *(tstar/t0
                0209      &                    + qstar/(1. _d 0/humid_fac+q))
ced0783fba Jean*0210                  xsq = max(sqrt(abs(1. _d 0 - 16. _d 0*z100ol)),1. _d 0)
                0211                  x      = sqrt(xsq)
                0212                  psx100  = -5. _d 0*z100ol*stable + (1. _d 0-stable)*
6f955c5285 Jean*0213      &                     (2. _d 0*log(5. _d -1*(1. _d 0+xsq)))
e900d93663 Jean*0214                  q100=ssq+qstar*(dlog(100. _d 0/zice)-psx100)
ced0783fba Jean*0215 
                0216       RETURN
                0217       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!