Back to home page

MITgcm
 
 

    

File indexing completed on 2018-03-02 18:37:47 UTC

view on githubraw file Latest commit add29e06 on 2018-01-31 20:35:05 UTC
b2ea1d2979 Jean*0001 ! ============================================================================
                0002 module surface_flux_mod
                0003 ! ============================================================================
                0004 
                0005 ! use             fms_mod, only: FATAL, close_file, mpp_pe, mpp_root_pe,        &
                0006 !                               write_version_number
                0007 ! use             fms_mod, only: file_exist, check_nml_error, open_namelist_file, stdlog
                0008 
                0009 use   monin_obukhov_mod, only: mo_drag, mo_profile
                0010 use  simple_sat_vapor_pres_mod, only: escomp, descomp
                0011 use      gcm_params_mod, only: gcm_LEN_MBUF, gcm_SQZ_R, gcm_stdMsgUnit
                0012 use       constants_mod, only: cp_air, hlv, stefan, rdgas, rvgas, grav
                0013 
                0014 implicit none
                0015 private
                0016 
                0017 ! ==== public interface ======================================================
                0018 public  surface_flux
                0019 ! ==== end of public interface ===============================================
                0020 
                0021 interface surface_flux
                0022 !    module procedure surface_flux_0d
                0023 !    module procedure surface_flux_1d
                0024     module procedure surface_flux_2d
                0025 end interface
                0026 
                0027 !-----------------------------------------------------------------------
                0028 
15388b052e Jean*0029 character(len=*), parameter :: version = '$Id: surface_flux_mod.F90,v 1.3 2017/08/11 20:48:51 jmc Exp $'
b2ea1d2979 Jean*0030 character(len=*), parameter :: tagname = '$Name:  $'
                0031 
                0032 logical :: do_init = .true.
                0033 
                0034 real, parameter :: d622   = rdgas/rvgas
                0035 real, parameter :: d378   = 1.-d622
                0036 real, parameter :: hlars  = hlv/rvgas
                0037 real, parameter :: gcp    = grav/cp_air
                0038 real, parameter :: kappa  = rdgas/cp_air
                0039 real            :: d608   = d378/d622
                0040       ! d608 set to zero at initialization if the use of
                0041       ! virtual temperatures is turned off in namelist
                0042 
                0043 ! ---- namelist with default values ------------------------------------------
                0044 logical :: no_neg_q         = .false.  ! for backwards compatibility
                0045 logical :: use_virtual_temp = .true.
                0046 logical :: alt_gustiness    = .false.
                0047 logical :: use_mixing_ratio = .false.
                0048 real    :: gust_const       =  1.0
                0049 
                0050 namelist /surface_flux_nml/ no_neg_q,         &
                0051                             use_virtual_temp, &
                0052                             alt_gustiness,    &
                0053                             gust_const,       &
                0054                             use_mixing_ratio
                0055 
                0056 contains
                0057 
                0058 ! ============================================================================
                0059 subroutine surface_flux_1d (                                           &
                0060      t_atm,     q_atm_in,   u_atm,     v_atm,     p_atm,     z_atm,    &
                0061      p_surf,    t_surf,     t_ca,      q_surf,                         &
                0062      u_surf,    v_surf,                                                &
                0063      rough_mom, rough_heat, rough_moist, gust,  &
                0064      flux_t, flux_q, flux_r, flux_u, flux_v,                &
                0065      cd_m,      cd_t,       cd_q,                                      &
                0066      w_atm,     u_star,     b_star,     q_star,                        &
                0067      dhdt_surf, dedt_surf,  dedq_surf,  drdt_surf,                     &
                0068      dhdt_atm,  dedq_atm,   dtaudv_atm,                                &
                0069 ! + slm Mar 28 2002 -- it is not necessary here since it is just cd_q*wind
                0070 !     drag_q,                                                           &
                0071 ! - slm Mar 28 2002
                0072      dt,        land,      avail, myThid  )
                0073 ! ============================================================================
                0074   ! ---- arguments -----------------------------------------------------------
                0075   logical, intent(in), dimension(:) :: land,  avail
                0076   real, intent(in),  dimension(:) :: &
                0077        t_atm,     q_atm_in,   u_atm,     v_atm,              &
                0078        p_atm,     z_atm,      t_ca,                          &
                0079        p_surf,    t_surf,     u_surf,    v_surf,  &
                0080        rough_mom, rough_heat, rough_moist,  gust
                0081   real, intent(out), dimension(:) :: &
                0082        flux_t,    flux_q,     flux_r,    flux_u,  flux_v,    &
                0083        dhdt_surf, dedt_surf,  dedq_surf, drdt_surf,          &
                0084        dhdt_atm,  dedq_atm,   dtaudv_atm,                    &
                0085        w_atm,     u_star,     b_star,    q_star,             &
                0086        cd_m,      cd_t,       cd_q
                0087   real, intent(inout), dimension(:) :: q_surf
                0088   real, intent(in) :: dt
                0089   integer, intent(in) :: myThid
                0090 
                0091   ! ---- local constants -----------------------------------------------------
                0092   ! temperature increment and its reciprocal value for comp. of derivatives
                0093   real, parameter:: del_temp=0.1, del_temp_inv=1.0/del_temp
                0094 
                0095   ! ---- local vars ----------------------------------------------------------
                0096   real, dimension(size(t_atm)) ::                          &
                0097        thv_atm,  th_atm,   tv_atm,    thv_surf,            &
                0098        e_sat,    e_sat1,   q_sat,     q_sat1,    p_ratio,  &
                0099        t_surf0,  t_surf1,  u_dif,     v_dif,               &
                0100        rho_drag, drag_t,    drag_m,   drag_q,    rho,      &
                0101        q_atm,    q_surf0
                0102 
                0103   integer :: i, nbad
                0104 
                0105   if (do_init) call surface_flux_init(myThid)
                0106 
                0107   !---- use local value of surf temp ----
                0108 
                0109   t_surf0 = 200.   !  avoids out-of-bounds in es lookup
                0110   where (avail)
                0111      where (land)
                0112         t_surf0 = t_ca
                0113      elsewhere
                0114         t_surf0 = t_surf
                0115      endwhere
                0116   endwhere
                0117 
                0118   t_surf1 = t_surf0 + del_temp
                0119 
                0120   call escomp ( t_surf0, e_sat  )  ! saturation vapor pressure
                0121   call escomp ( t_surf1, e_sat1 )  ! perturbed  vapor pressure
                0122 
                0123   if(use_mixing_ratio) then
                0124     ! surface mixing ratio at saturation
                0125     q_sat   = d622*e_sat /(p_surf-e_sat )
                0126     q_sat1  = d622*e_sat1/(p_surf-e_sat1)
                0127   else
                0128     ! surface specific humidity at saturation
                0129     q_sat   = d622*e_sat /(p_surf-d378*e_sat )
                0130     q_sat1  = d622*e_sat1/(p_surf-d378*e_sat1)
                0131   endif
                0132 
                0133   ! initilaize surface air humidity according to surface type
                0134   where (land)
                0135      q_surf0 = q_surf ! land calculates it
                0136   elsewhere
                0137      q_surf0 = q_sat  ! everything else assumes saturated sfc humidity
                0138   endwhere
                0139 
                0140   ! check for negative atmospheric humidities
                0141   where(avail) q_atm = q_atm_in
                0142   if(no_neg_q) then
                0143      where(avail .and. q_atm_in < 0.0) q_atm = 0.0
                0144   endif
                0145 
                0146   ! generate information needed by monin_obukhov
                0147   where (avail)
                0148      p_ratio = (p_surf/p_atm)**kappa
                0149 
                0150      tv_atm  = t_atm  * (1.0 + d608*q_atm)     ! virtual temperature
                0151      th_atm  = t_atm  * p_ratio  ! potential T, using p_surf as refernce
                0152      thv_atm = tv_atm * p_ratio  ! virt. potential T, using p_surf as reference
                0153      thv_surf= t_surf0 * (1.0 + d608*q_surf0 ) ! surface virtual (potential) T
                0154 !     thv_surf= t_surf0  ! surface virtual (potential) T -- just for testing tun off the q_surf
                0155 
                0156      u_dif = u_surf - u_atm      ! velocity components relative to surface
                0157      v_dif = v_surf - v_atm
                0158   endwhere
                0159 
                0160   if(alt_gustiness) then
                0161      where(avail) &
                0162           w_atm = max(sqrt(u_dif*u_dif + v_dif*v_dif), gust_const)
                0163   else
                0164      where(avail) &
                0165           w_atm = sqrt(u_dif*u_dif + v_dif*v_dif + gust*gust)
                0166   endif
                0167 
                0168   !  monin-obukhov similarity theory
                0169   call mo_drag (thv_atm, thv_surf, z_atm,                  &
                0170        rough_mom, rough_heat, rough_moist, w_atm,          &
                0171        cd_m, cd_t, cd_q, u_star, b_star, myThid, avail     )
                0172 
                0173   where (avail)
                0174      ! surface layer drag coefficients
                0175      drag_t = cd_t * w_atm
                0176      drag_q = cd_q * w_atm
                0177      drag_m = cd_m * w_atm
                0178 
                0179      ! density
                0180      rho = p_atm / (rdgas * tv_atm)
                0181 
                0182      ! sensible heat flux
                0183      rho_drag = cp_air * drag_t * rho
                0184      flux_t = rho_drag * (t_surf0 - th_atm)  ! flux of sensible heat (W/m**2)
                0185      dhdt_surf =  rho_drag                   ! d(sensible heat flux)/d(surface temperature)
                0186      dhdt_atm  = -rho_drag*p_ratio           ! d(sensible heat flux)/d(atmos temperature)
                0187 
                0188      ! evaporation
                0189      rho_drag  =  drag_q * rho
                0190      flux_q    =  rho_drag * (q_surf0 - q_atm) ! flux of water vapor  (Kg/(m**2 s))
                0191      where (land)
                0192         dedq_surf = rho_drag
                0193         dedt_surf = 0
                0194      elsewhere
                0195         dedq_surf = 0
                0196         dedt_surf =  rho_drag * (q_sat1 - q_sat) *del_temp_inv
                0197      endwhere
                0198 
                0199      dedq_atm  = -rho_drag   ! d(latent heat flux)/d(atmospheric mixing ratio)
                0200 
                0201      q_star = flux_q / (u_star * rho)             ! moisture scale
                0202      ! ask Chris and Steve K if we still want to keep this for diagnostics
                0203      q_surf = q_atm + flux_q / (rho*cd_q*w_atm)   ! surface specific humidity
                0204 
                0205      ! upward long wave radiation
                0206      flux_r    =   stefan*t_surf**4               ! (W/m**2)
                0207      drdt_surf = 4*stefan*t_surf**3               ! d(upward longwave)/d(surface temperature)
                0208 
                0209      ! stresses
                0210      rho_drag   = drag_m * rho
                0211      flux_u     = rho_drag * u_dif   ! zonal      component of stress (Nt/m**2)
                0212      flux_v     = rho_drag * v_dif   ! meridional component of stress
                0213      dtaudv_atm = -rho_drag          ! d(stress component)/d(atmos wind)
                0214 
                0215   elsewhere
                0216      ! zero-out un-available data in output only fields
                0217      flux_t     = 0.0
                0218      flux_q     = 0.0
                0219      flux_r     = 0.0
                0220      flux_u     = 0.0
                0221      flux_v     = 0.0
                0222      dhdt_surf  = 0.0
                0223      dedt_surf  = 0.0
                0224      dedq_surf  = 0.0
                0225      drdt_surf  = 0.0
                0226      dhdt_atm   = 0.0
                0227      dedq_atm   = 0.0
                0228      dtaudv_atm = 0.0
                0229      u_star     = 0.0
                0230      b_star     = 0.0
                0231      q_star     = 0.0
                0232      q_surf     = 0.0
                0233      w_atm      = 0.0
                0234   endwhere
                0235 
                0236 end subroutine surface_flux_1d
                0237 
                0238 !#######################################################################
                0239 
                0240 subroutine surface_flux_0d (                                                 &
                0241      t_atm_0,     q_atm_0,      u_atm_0,     v_atm_0,   p_atm_0, z_atm_0,    &
                0242      p_surf_0,    t_surf_0,     t_ca_0,      q_surf_0,                       &
                0243      u_surf_0,    v_surf_0,                                                  &
                0244      rough_mom_0, rough_heat_0, rough_moist_0, gust_0,                       &
                0245      flux_t_0,    flux_q_0,     flux_r_0,    flux_u_0,  flux_v_0,            &
                0246      cd_m_0,      cd_t_0,       cd_q_0,                                      &
                0247      w_atm_0,     u_star_0,     b_star_0,     q_star_0,                      &
                0248      dhdt_surf_0, dedt_surf_0,  dedq_surf_0,  drdt_surf_0,                   &
                0249      dhdt_atm_0,  dedq_atm_0,   dtaudv_atm_0,                                &
                0250      dt,          land_0,       avail_0, myThid  )
                0251 
                0252   ! ---- arguments -----------------------------------------------------------
                0253   logical, intent(in) :: land_0,  avail_0
                0254   real, intent(in) :: &
                0255        t_atm_0,     q_atm_0,      u_atm_0,     v_atm_0,              &
                0256        p_atm_0,     z_atm_0,      t_ca_0,                          &
                0257        p_surf_0,    t_surf_0,     u_surf_0,    v_surf_0,  &
                0258        rough_mom_0, rough_heat_0, rough_moist_0,  gust_0
                0259   real, intent(out) :: &
                0260        flux_t_0,    flux_q_0,     flux_r_0,    flux_u_0,  flux_v_0,    &
                0261        dhdt_surf_0, dedt_surf_0,  dedq_surf_0, drdt_surf_0,          &
                0262        dhdt_atm_0,  dedq_atm_0,   dtaudv_atm_0,                    &
                0263        w_atm_0,     u_star_0,     b_star_0,    q_star_0,             &
                0264        cd_m_0,      cd_t_0,       cd_q_0
                0265   real, intent(inout) :: q_surf_0
                0266   real, intent(in)    :: dt
                0267   integer, intent(in) :: myThid
                0268 
                0269   ! ---- local vars ----------------------------------------------------------
                0270   logical, dimension(1) :: land,  avail
                0271   real, dimension(1) :: &
                0272        t_atm,     q_atm,      u_atm,     v_atm,              &
                0273        p_atm,     z_atm,      t_ca,                          &
                0274        p_surf,    t_surf,     u_surf,    v_surf,  &
                0275        rough_mom, rough_heat, rough_moist,  gust
                0276   real, dimension(1) :: &
                0277        flux_t,    flux_q,     flux_r,    flux_u,  flux_v,    &
                0278        dhdt_surf, dedt_surf,  dedq_surf, drdt_surf,          &
                0279        dhdt_atm,  dedq_atm,   dtaudv_atm,                    &
                0280        w_atm,     u_star,     b_star,    q_star,             &
                0281        cd_m,      cd_t,       cd_q
                0282   real, dimension(1) :: q_surf
                0283 
                0284   avail = .true.
                0285 
                0286   t_atm(1)       = t_atm_0
                0287   q_atm(1)       = q_atm_0
                0288   u_atm(1)       = u_atm_0
                0289   v_atm(1)       = v_atm_0
                0290   p_atm(1)       = p_atm_0
                0291   z_atm(1)       = z_atm_0
                0292   p_surf(1)      = p_surf_0
                0293   t_surf(1)      = t_surf_0
68d8f0461a Jean*0294   t_ca(1)        = t_ca_0
                0295   q_surf(1)      = q_surf_0
b2ea1d2979 Jean*0296   u_surf(1)      = u_surf_0
                0297   v_surf(1)      = v_surf_0
                0298   rough_mom(1)   = rough_mom_0
                0299   rough_heat(1)  = rough_heat_0
                0300   rough_moist(1) = rough_moist_0
                0301   gust(1)        = gust_0
68d8f0461a Jean*0302   land(1)        = land_0
b2ea1d2979 Jean*0303 
                0304   call surface_flux_1d (                                           &
                0305        t_atm,     q_atm,      u_atm,     v_atm,     p_atm,     z_atm,    &
                0306        p_surf,    t_surf,     t_ca,      q_surf,                         &
                0307        u_surf,    v_surf,                                                &
                0308        rough_mom, rough_heat, rough_moist, gust,  &
                0309        flux_t, flux_q, flux_r, flux_u, flux_v,                &
                0310        cd_m,      cd_t,       cd_q,                                      &
                0311        w_atm,     u_star,     b_star,     q_star,                        &
                0312        dhdt_surf, dedt_surf,  dedq_surf,  drdt_surf,                     &
                0313        dhdt_atm,  dedq_atm,   dtaudv_atm,                                &
                0314        dt,        land,      avail, myThid  )
                0315 
                0316   flux_t_0     = flux_t(1)
                0317   flux_q_0     = flux_q(1)
                0318   flux_r_0     = flux_r(1)
                0319   flux_u_0     = flux_u(1)
                0320   flux_v_0     = flux_v(1)
68d8f0461a Jean*0321   cd_m_0       = cd_m(1)
                0322   cd_t_0       = cd_t(1)
                0323   cd_q_0       = cd_q(1)
                0324   w_atm_0      = w_atm(1)
                0325   u_star_0     = u_star(1)
                0326   b_star_0     = b_star(1)
                0327   q_star_0     = q_star(1)
                0328   q_surf_0     = q_surf(1)
b2ea1d2979 Jean*0329   dhdt_surf_0  = dhdt_surf(1)
                0330   dedt_surf_0  = dedt_surf(1)
                0331   dedq_surf_0  = dedq_surf(1)
68d8f0461a Jean*0332   drdt_surf_0  = drdt_surf(1)
b2ea1d2979 Jean*0333   dhdt_atm_0   = dhdt_atm(1)
                0334   dedq_atm_0   = dedq_atm(1)
                0335   dtaudv_atm_0 = dtaudv_atm(1)
                0336 
                0337 end subroutine surface_flux_0d
                0338 
                0339 subroutine surface_flux_2d (                                           &
                0340      t_atm,     q_atm_in,   u_atm,     v_atm,     p_atm,     z_atm,    &
                0341      p_surf,    t_surf,     t_ca,      q_surf,                         &
                0342      u_surf,    v_surf,                                                &
                0343      rough_mom, rough_heat, rough_moist, gust,                         &
                0344      flux_t,    flux_q,     flux_r,    flux_u,    flux_v,              &
                0345      cd_m,      cd_t,       cd_q,                                      &
                0346      w_atm,     u_star,     b_star,     q_star,                        &
                0347      dhdt_surf, dedt_surf,  dedq_surf,  drdt_surf,                     &
                0348      dhdt_atm,  dedq_atm,   dtaudv_atm,                                &
                0349      dt,        land,       avail, myThid  )
                0350 
                0351   ! ---- arguments -----------------------------------------------------------
                0352   logical, intent(in), dimension(:,:) :: land,  avail
                0353   real, intent(in),  dimension(:,:) :: &
                0354        t_atm,     q_atm_in,   u_atm,     v_atm,              &
                0355        p_atm,     z_atm,      t_ca,                          &
                0356        p_surf,    t_surf,     u_surf,    v_surf,  &
                0357        rough_mom, rough_heat, rough_moist,  gust
                0358   real, intent(out), dimension(:,:) :: &
                0359        flux_t,    flux_q,     flux_r,    flux_u,  flux_v,    &
                0360        dhdt_surf, dedt_surf,  dedq_surf, drdt_surf,          &
                0361        dhdt_atm,  dedq_atm,   dtaudv_atm,                    &
                0362        w_atm,     u_star,     b_star,    q_star,             &
                0363        cd_m,      cd_t,       cd_q
                0364   real, intent(inout), dimension(:,:) :: q_surf
                0365   real, intent(in) :: dt
                0366   integer, intent (in) :: myThid
                0367 
                0368   ! ---- local vars -----------------------------------------------------------
                0369   integer :: j
                0370   do j = 1, size(t_atm,2)
                0371      call surface_flux_1d (                                           &
                0372           t_atm(:,j),     q_atm_in(:,j),   u_atm(:,j),     v_atm(:,j),     p_atm(:,j),     z_atm(:,j),    &
                0373           p_surf(:,j),    t_surf(:,j),     t_ca(:,j),      q_surf(:,j),                         &
                0374           u_surf(:,j),    v_surf(:,j),                                                &
                0375           rough_mom(:,j), rough_heat(:,j), rough_moist(:,j), gust(:,j),                         &
                0376           flux_t(:,j),    flux_q(:,j),     flux_r(:,j),    flux_u(:,j),    flux_v(:,j),              &
                0377           cd_m(:,j),      cd_t(:,j),       cd_q(:,j),                                      &
                0378           w_atm(:,j),     u_star(:,j),     b_star(:,j),     q_star(:,j),                        &
                0379           dhdt_surf(:,j), dedt_surf(:,j),  dedq_surf(:,j),  drdt_surf(:,j),                     &
                0380           dhdt_atm(:,j),  dedq_atm(:,j),   dtaudv_atm(:,j),                                &
                0381           dt,             land(:,j),       avail(:,j), myThid  )
                0382   end do
                0383 end subroutine surface_flux_2d
                0384 
                0385 ! ============================================================================
                0386 subroutine surface_flux_init (myThid)
                0387 ! ============================================================================
                0388 ! initializes surface flux module
                0389   ! ---- local vars ----------------------------------------------------------
                0390  integer, intent(in) :: myThid
                0391  integer :: unit, ierr, io
                0392 !-- added:
                0393  integer         :: iUnit
                0394  CHARACTER*(gcm_LEN_MBUF) :: msgBuf
                0395 
                0396   ! read namelist
                0397 
                0398 !    _BARRIER
                0399 !    _BEGIN_MASTER(myThid)
                0400      CALL BARRIER(myThid)
                0401      IF ( myThid.EQ.1 ) THEN
                0402 
                0403      WRITE(msgBuf,'(A)') 'SURFACE_FLUX_INIT: opening data.atm_gray'
                0404      CALL PRINT_MESSAGE( msgBuf, gcm_stdMsgUnit, gcm_SQZ_R, myThid )
                0405      CALL OPEN_COPY_DATA_FILE(                                      &
                0406                            'data.atm_gray', 'SURFACE_FLUX_INIT',       &
                0407                            iUnit,                                   &
                0408                            myThid )
                0409 !    Read parameters from open data file
                0410      READ(UNIT=iUnit,NML=surface_flux_nml)
                0411      WRITE(msgBuf,'(A)')                                            &
                0412           'SURFACE_FLUX_INIT: finished reading data.atm_gray'
                0413      CALL PRINT_MESSAGE( msgBuf, gcm_stdMsgUnit, gcm_SQZ_R, myThid )
                0414 !    Close the open data file
15388b052e Jean*0415 #ifdef SINGLE_DISK_IO
b2ea1d2979 Jean*0416      CLOSE(iUnit)
15388b052e Jean*0417 #else
                0418      CLOSE(iUnit,STATUS='DELETE')
                0419 #endif /* SINGLE_DISK_IO */
b2ea1d2979 Jean*0420 
                0421 !  if ( file_exist('input.nml')) then
                0422 !     unit = open_namelist_file ()
                0423 !     ierr=1;
                0424 !     do while (ierr /= 0)
                0425 !        read  (unit, nml=surface_flux_nml, iostat=io, end=10)
                0426 !        ierr = check_nml_error(io,'surface_flux_nml')
                0427 !     enddo
                0428 !10   call close_file (unit)
                0429 !  endif
                0430 
                0431 !  ! write version number
                0432 !  call write_version_number(version, tagname)
                0433 
                0434 !  if ( mpp_pe() == mpp_root_pe() ) write (stdlog(), nml=surface_flux_nml)
                0435 
                0436 !  if(.not. use_virtual_temp) d608 = 0.0
                0437 
                0438     do_init = .false.
                0439 
                0440     ENDIF
                0441     CALL BARRIER (myThid)
                0442 
                0443 end subroutine surface_flux_init
                0444 
                0445 end module surface_flux_mod
[ 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!