MODULE GwddSchemeCPTEC PRIVATE ! 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(31)! Kind for 128-bits Real Numbers REAL(r8),PARAMETER :: SHR_CONST_CPDAIR = 1.00464e3_r8 ! specific heat of dry air ~ J/kg/K REAL(R8),PARAMETER :: SHR_CONST_BOLTZ = 1.38065e-23_R8 ! Boltzmann's constant ~ J/K/molecule REAL(R8),PARAMETER :: SHR_CONST_AVOGAD = 6.02214e26_R8 ! Avogadro's number ~ molecules/kmole REAL(R8),PARAMETER :: SHR_CONST_MWDAIR = 28.966_R8 ! molecular weight dry air ~ kg/kmole REAL(r8),PARAMETER :: SHR_CONST_MWWV = 18.016_r8 ! molecular weight water vapor REAL(R8),PARAMETER :: SHR_CONST_CPWV = 1.810e3_R8 ! specific heat of water vap ~ J/kg/K REAL(R8),PARAMETER :: SHR_CONST_G = 9.80616_R8 ! acceleration of gravity ~ m/s^2 REAL(R8),PARAMETER :: SHR_CONST_RGAS = SHR_CONST_AVOGAD*SHR_CONST_BOLTZ ! Universal gas constant ~ J/K/kmole REAL(R8),PARAMETER :: SHR_CONST_RDAIR = SHR_CONST_RGAS/SHR_CONST_MWDAIR ! Dry air gas constant ~ J/K/kg REAL(r8),PARAMETER :: SHR_CONST_RWV = SHR_CONST_RGAS/SHR_CONST_MWWV ! Water vapor gas constant ~ J/K/kg REAL (KIND=r8), PARAMETER :: cp = 1004.6_r8! specific heat of air (j/kg/k) REAL (KIND=r8), PARAMETER :: gasr = 287.05_r8! gas constant of dry air (j/kg/k) REAL(r8), PUBLIC, PARAMETER :: gravit = shr_const_g ! gravitational acceleration REAL(r8), PUBLIC, PARAMETER :: cpair = shr_const_cpdair ! specific heat of dry air (J/K/kg) REAL(r8), PUBLIC, PARAMETER :: rair = shr_const_rdair ! Gas constant for dry air (J/K/kg) REAL(r8), PUBLIC, PARAMETER :: cpwv = shr_const_cpwv REAL(r8), PUBLIC, PARAMETER :: zvir = SHR_CONST_RWV/rair - 1 ! rh2o/rair - 1 REAL(KIND=r8), PARAMETER :: Pi = 3.14159265358979323846_r8 ! gw_vert ---- gw_satn ! | ! gw_wake ! !*L------------------COMDECK C_A---------------------------------------- ! History: ! Version Date Comment. ! 5.0 07/05/99 Replace variable A by more meaningful name for ! conversion to C-P 'C' dynamics grid. R. Rawlins ! 5.1 07/03/00 Convert to Fixed/Free format. P. Selwood ! Mean radius of Earth in m. REAL(KIND=r8), PARAMETER :: Earth_Radius = 6371229.0_r8 REAL(KIND=r8), PARAMETER :: recip_a2 =1.0_r8/(earth_radius*earth_radius) ! 1/(radius of earth)^2 !*L------------------COMDECK C_G---------------------------------------- ! G IS MEAN ACCEL DUE TO GRAVITY AT EARTH'S SURFACE REAL(KIND=r8), PARAMETER :: G = 9.80665_r8 ! lambda=100000 m ! kwv =2*pi/lambda ! ! kwv = 6.28e-5_r8 ! 100 km wave length !GWD_FRC_max = 6.0 !GWD_FRC_min = 2.0 REAL(KIND=r8), PARAMETER :: Lambda_GWAVE_max = 2.0E+05_r8!highest wave length REAL(KIND=r8), PARAMETER :: Lambda_GWAVE_min = 2.5E+02_r8!lowet wave length REAL(KIND=r8) :: kay ! surface stress constant ( m**-1) REAL(KIND=r8) ,PARAMETER :: frc=5.5_r8! ! critical Froude number LOGICAL ,PARAMETER :: l_fix_gwsatn=.TRUE.! if true then invoke minor bug fixes in gwsatn LOGICAL ,PARAMETER :: l_gwd_40km=.TRUE.!.FALSE. ! if true then don't apply GWD above 40km LOGICAL ,PARAMETER :: l_taus_scale=.TRUE.! if true then surface stress is made to ! ! depend on the low level Froude number ! Number of standard deviations above the mean orography of top ! of sub-grid mountains REAL(KIND=r8),PARAMETER :: NSIGMA = 3.50_r8 !*---------------------------------------------------------------------- REAL(KIND=r8), ALLOCATABLE :: sigkiv(:) REAL(KIND=r8), ALLOCATABLE :: sl(:) REAL(KIND=r8), ALLOCATABLE :: delsig(:) REAL(KIND=r8), ALLOCATABLE :: si(:) !*---------------------------------------------------------------------- PUBLIC :: Init_Gwave PUBLIC :: g_wave CONTAINS SUBROUTINE Init_Gwave(kmax,sl_in,delsig_in,si_in ) INTEGER , INTENT(IN ) :: kmax REAL(KIND=r8), INTENT(IN ) :: sl_in (kmax) REAL(KIND=r8), INTENT(IN ) :: delsig_in(kmax) REAL(KIND=r8), INTENT(IN ) :: si_in (kmax+1) REAL(KIND=r8) :: akappa REAL(KIND=r8) :: sigk(kmax) REAL(KIND=r8) :: ku!is the hisghest wavenumber lying outside the range of trapped lee waves REAL(KIND=r8) :: kl!is the lowest unresolved wave number of model REAL(KIND=r8) :: pow1 REAL(KIND=r8) :: pow2 INTEGER :: k ALLOCATE(sigkiv(kmax)) ALLOCATE(sl(kmax)) ALLOCATE(delsig(kmax)) ALLOCATE(si(kmax+1)) sl=sl_in si=si_in delsig=delsig_in akappa=gasr/cp ! DO k = 1, kmax sigk (k)=sl(k)**akappa sigkiv(k)=1.0_r8/sigk(k) END DO pow1=1.0_r8/2.0_r8 pow2=3.0_r8/2.0_r8 ku=(2.0_r8*Pi)/Lambda_GWAVE_min !is the hisghest wavenumber lying outside the range of trapped lee waves kl=(2.0_r8*Pi)/Lambda_GWAVE_max !is the lowest unresolved wave number of model kay=(1.0_r8/3.0_r8)*(((ku**(pow2)) - (kl**(pow2)))/ & ((ku**(pow1)) - (kl**(pow1)))) !WRITE(0,*)'Init_Gwave',kay END SUBROUTINE Init_Gwave SUBROUTINE g_wave( & kMax , & !INTEGER, INTENT(IN ) :: kMax nCols , & !INTEGER, INTENT(IN ) :: nCols ! number of points per row gps , & !REAL(KIND=r8), INTENT(IN ) :: gps (nCols) !REAL(r8), INTENT(in) :: gps(nCols) Pa gt , & !REAL(KIND=r8), INTENT(IN ) :: gt (nCols,kMax) !REAL(r8), INTENT(in) :: gt (nCols,kMax) gq , & !REAL(KIND=r8), INTENT(IN ) :: gq (nCols,kMax) !REAL(r8), INTENT(in) :: gq (nCols,kMax) gu , & !REAL(KIND=r8), INTENT(IN ) :: gu (nCols,kMax) !REAL(r8), INTENT(in) :: gu (nCols,kMax) gv , & !REAL(KIND=r8), INTENT(IN ) :: gv (nCols,kMax) !REAL(r8), INTENT(in) :: gv (nCols,kMax) topo , & !REAL(KIND=r8), INTENT(IN ) :: topo (nCols) !colrad , & !REAL(KIND=r8), INTENT(IN ) :: colrad (nCols) sd_orog , & !REAL(KIND=r8) ,INTENT(in ) :: sd_orog (nCols)! standard deviation of orography (m) orog_grad_xx , & !REAL(KIND=r8), INTENT(in ) :: orog_grad_xx(nCols)! dh/dx squared gradient orography orog_grad_yy , & !REAL(KIND=r8), INTENT(in ) :: orog_grad_yy(nCols)! dh/dy squared gradient orography orog_grad_xy , & !REAL(KIND=r8), INTENT(in ) :: orog_grad_xy(nCols)! (dh/dx)(dh/dy) gradient orography imask , & !INTEGER , INTENT(in ) :: imask (nCols)! index for land points timestep , & !REAL(KIND=r8) ,INTENT(in ) :: timestep ! timestep (s) du_dt , & !REAL(KIND=r8) ,INTENT(out):: du_dt(nCols,kMax) ! total GWD du/dt on land/theta dv_dt , & !REAL(KIND=r8) ,INTENT(out):: dv_dt(nCols,kMax) ! total GWD dv/dt on land/theta iret ) !INTEGER ,INTENT(OUT ) :: iret ! return code : iret=0 normal exit IMPLICIT NONE ! ! Description: ! 1) Interpolate winds to theta points ! gather data for land points only ! 2) Call surface stress routine ! 3) Calculate stress profiles due to different components of the ! scheme. Calculate associated wind increments. ! 4) Interpolate acceleration to wind points and update winds ! 5) Gather diagnostics from land points and interpolate from p to ! u,v staggering on 'c' grid ! ! current code owner: S.Webster ! ! history: ! version date comment ! ------- ---- ------- ! 5.2 15/11/00 original deck. Stuart Webster ! 5.3 16/10/01 Remove code no longer required because of ! simplifications to scheme. Stuart Webster ! 5.3 27/07/01 Permit a LAM model to run with GWD scheme on. ! S. Cusack ! 5.4 28/08/02 Introduce numerical limiter for flow blocking ! scheme. S. Webster ! 6.2 17/05/06 Remove t_to_p reference. P.Selwood. ! 6.2 21/02/06 Introduce surface stress diagnostics. S.Webster ! 6.2 21/02/06 Pass l_taus_scale switch thru' to GWSURF4A ! and l_fix_gwsatn thru to GWSATN4A. S. Webster ! ! ! language: fortran 77 + common extensions. ! this code is written to umdp3 v6 programming standards. ! ! suitable for single column use, with calls to: uv_to_p removed ! p_to_uv removed ! suitable for rotated grids ! ! global variables (*called comdecks etc...): ! ! SUBROUTINE ARGUMENTS ! INTEGER, INTENT(IN ) :: kMax INTEGER, INTENT(IN ) :: nCols ! number of points per row REAL(KIND=r8), INTENT(IN ) :: gps (nCols) !REAL(r8), INTENT(in) :: gps(nCols) Pa REAL(KIND=r8), INTENT(IN ) :: gt (nCols,kMax) !REAL(r8), INTENT(in) :: gt (nCols,kMax) REAL(KIND=r8), INTENT(IN ) :: gq (nCols,kMax) !REAL(r8), INTENT(in) :: gq (nCols,kMax) REAL(KIND=r8), INTENT(IN ) :: gu (nCols,kMax) !REAL(r8), INTENT(in) :: gu (nCols,kMax) REAL(KIND=r8), INTENT(IN ) :: gv (nCols,kMax) !REAL(r8), INTENT(in) :: gv (nCols,kMax) REAL(KIND=r8), INTENT(IN ) :: topo (nCols) !REAL(KIND=r8), INTENT(IN ) :: colrad (nCols) INTEGER(KIND=i8), INTENT(in ) :: imask (nCols)! index for land points REAL(KIND=r8), INTENT(in ) :: orog_grad_xx(nCols)! dh/dx squared gradient orography REAL(KIND=r8), INTENT(in ) :: orog_grad_yy(nCols)! dh/dy squared gradient orography REAL(KIND=r8), INTENT(in ) :: orog_grad_xy(nCols)! (dh/dx)(dh/dy) gradient orography REAL(KIND=r8) ,INTENT(in ) :: sd_orog (nCols)! standard deviation of orography (m) REAL(KIND=r8) ,INTENT(in ) :: timestep ! timestep (s) INTEGER ,INTENT(OUT ) :: iret ! return code : iret=0 normal exit ! Start of full field diagnostic arrays. This space is allocated ! in GWD_CTL2 only if a diagnostic is called. ! REAL(KIND=r8) ,INTENT(out):: du_dt(nCols,kMax) ! total GWD du/dt on land/theta REAL(KIND=r8) ,INTENT(out):: dv_dt(nCols,kMax) ! total GWD dv/dt on land/theta REAL(KIND=r8) :: theta (nCols,kMax)! primary theta field (K) REAL(KIND=r8) :: rho (1:nCols,kMax)! density*(radius of earth)^2 REAL(KIND=r8) :: r_rho_levels (1:nCols, kMax)! height of rho level above earth's centre (m) REAL(KIND=r8) :: r_theta_levels (1:nCols,0:kMax)! height of theta level above earth's centre (m) ! ! ! END of full field diagnostic arrays ! Below are the stash flags for calculating diagnostics: ! INTEGER :: land_points ! number of land points LOGICAL :: stress_ud_on !u stress LOGICAL :: stress_ud_p_on !u on press stress LOGICAL :: stress_vd_on !v stress LOGICAL :: stress_ud_satn_on !u satn stress LOGICAL :: stress_vd_satn_on !v satn stress LOGICAL :: stress_ud_wake_on !u wake stress LOGICAL :: stress_vd_wake_on !v wake stress LOGICAL :: du_dt_satn_on !u accel (saturation) LOGICAL :: du_dt_satn_p_on !u accel (saturation) LOGICAL :: dv_dt_satn_on !v accel (saturation) LOGICAL :: du_dt_wake_on !u accel blocked flow LOGICAL :: dv_dt_wake_on !v accel blocked flow LOGICAL :: u_s_d_on !u_s_d diag switch LOGICAL :: v_s_d_on !v_s_d diag switch LOGICAL :: nsq_s_d_on !nsq_s_d diag switch LOGICAL :: fr_d_on !fr_d switch LOGICAL :: bld_d_on !bld_d switch LOGICAL :: bldt_d_on !bldt_d switch LOGICAL :: num_lim_d_on !num_lim_d switch LOGICAL :: num_fac_d_on !num_fac_d switch LOGICAL :: tausx_d_on !tausx_d switch LOGICAL :: tausy_d_on !tausy_d switch LOGICAL :: taus_scale_d_on !taus_scale_d switch ! ! The integers below are set to size nCols if the corresponding ! diagnostic is called or to 1 if it is not. These are set in GWD_CTL2 ! INTEGER :: points_stress_ud INTEGER :: points_stress_vd INTEGER :: points_stress_ud_satn INTEGER :: points_stress_vd_satn INTEGER :: points_stress_ud_wake INTEGER :: points_stress_vd_wake INTEGER :: points_du_dt_satn INTEGER :: points_dv_dt_satn INTEGER :: points_du_dt_wake INTEGER :: points_dv_dt_wake INTEGER :: points_u_s_d INTEGER :: points_v_s_d INTEGER :: points_nsq_s_d INTEGER :: points_fr_d INTEGER :: points_bld_d INTEGER :: points_bldt_d INTEGER :: points_num_lim_d INTEGER :: points_num_fac_d INTEGER :: points_tausx_d INTEGER :: points_tausy_d INTEGER :: points_taus_scale_d ! REAL(KIND=r8) :: stress_ud (nCols ,0:kMax) !u total stress REAL(KIND=r8) :: stress_vd (nCols ,0:kMax) !v total stress REAL(KIND=r8) :: stress_ud_satn(nCols ,0:kMax)!u satn stress REAL(KIND=r8) :: stress_vd_satn(nCols ,0:kMax)!v satn stress REAL(KIND=r8) :: stress_ud_wake(nCols ,0:kMax)!u wake stress REAL(KIND=r8) :: stress_vd_wake(nCols ,0:kMax)!v wake stress REAL(KIND=r8) :: du_dt_satn (nCols ,kMax)!u acceln (saturation) REAL(KIND=r8) :: dv_dt_satn (nCols ,kMax)!v acceln (saturation) REAL(KIND=r8) :: du_dt_land (nCols,kMax) ! total GWD du/dt on land/theta REAL(KIND=r8) :: dv_dt_land (nCols,kMax) ! total GWD dv/dt on land/theta REAL(KIND=r8) :: du_dt_wake (nCols ,kMax)!u acceln (blocked flow) REAL(KIND=r8) :: dv_dt_wake (nCols ,kMax)!v acceln (blocked flow) REAL(KIND=r8) :: u_s_d (nCols ) ! u_s diag at theta pts REAL(KIND=r8) :: v_s_d (nCols ) ! v_s diag at theta pts REAL(KIND=r8) :: nsq_s_d (nCols ) ! nsq_s diag at theta pts REAL(KIND=r8) :: fr_d (nCols ) ! Fr diag at theta pts REAL(KIND=r8) :: bld_d (nCols ) ! blocked layer depth at theta pts REAL(KIND=r8) :: bldt_d (nCols ) ! % of time blocked layer diagnosed REAL(KIND=r8) :: num_lim_d (nCols ) ! % of time numerical limiter invoked REAL(KIND=r8) :: num_fac_d (nCols ) ! % redn. of flow-blocking stress after numerical limiter invoked REAL(KIND=r8) :: tausx_d (nCols)!x-component of surface stress REAL(KIND=r8) :: tausy_d (nCols)!y-component of surface stress REAL(KIND=r8) :: taus_scale_d (nCols)! Factor surface stress scaled by ! if Froude no. dependence is on. ! LOGICAL ,intent(in ):: at_extremity(4) ! indicates if this processor is at north, ! ! south, east or west of the processor grid ! FLDTYPE definitions for the different field types recognised on the ! decomposition INTEGER,PARAMETER:: Nfld_max=7 ! maximum number of field types INTEGER,PARAMETER:: fld_type_p=1 ! grid on P points INTEGER,PARAMETER:: fld_type_u=2 ! grid on U points INTEGER,PARAMETER:: fld_type_v=3 ! grid on V points INTEGER,PARAMETER:: fld_type_comp_wave = 4 ! Compressed WAM Wave Field INTEGER,PARAMETER:: fld_type_full_wave = 5 ! Uncompressed WAM Wave Field INTEGER,PARAMETER:: fld_type_rim_wave = 6 ! Boundary data for WAM Wave Field INTEGER,PARAMETER:: fld_type_r=7 ! grid on river points INTEGER,PARAMETER:: fld_type_unknown=-1! non-standard grid ! FLDTYPE end ! DOMTYP contains different model domain types ! ! Author : P.Burton ! History: ! Version Date Comment. ! 5.0 15/04/99 New comdeck ! 5.2 15/11/00 add bi_cyclic_lam domain A. Malcolm INTEGER,PARAMETER:: mt_global = 1 INTEGER,PARAMETER:: mt_lam = 2 INTEGER,PARAMETER:: mt_cyclic_lam = 3 INTEGER,PARAMETER:: mt_bi_cyclic_lam = 4 INTEGER,PARAMETER:: mt_single_column = 5 ! DOMTYP end !-------------------------------------------------------------------- ! LOCAL DYNAMIC ARRAYS: !-------------------------------------------------------------------- INTEGER :: k_top (nCols) ! model level at mountain tops - ! ! exact definition given in gwsurf INTEGER :: k_top_max ! max(k_top) ! parameters for mpp INTEGER,PARAMETER :: pnorth= 1 ! north processor address in the neighbor array INTEGER,PARAMETER :: peast = 2 ! east processor address in the neighbor array INTEGER,PARAMETER :: psouth= 3 ! south processor address in the neighbor array INTEGER,PARAMETER :: pwest = 4! west processor address in the neighbor array INTEGER,PARAMETER :: nodomain = -1! value in neighbor array if the domain has ! no neighbor in this direction. otherwise ! the value will be the tid of the neighbor INTEGER :: i,k,points! loop counters in routine ! Work arrays REAL(KIND=r8) :: work_u (nCols ,kMax) REAL(KIND=r8) :: work_v (nCols ,kMax) ! REAL :: work_halo (1:nCols,0:kMax) ! REAL :: work_on_v_grid (nCols,kMax) REAL(KIND=r8) :: up_land (nCols,kMax)! interpolated u on theta grid REAL(KIND=r8) :: vp_land (nCols,kMax)! interpolated V on theta grid REAL(KIND=r8) :: theta_land (nCols,kMax)! land theta field on theta levels REAL(KIND=r8) :: r_rho_levels_land(nCols,kMax)! land field of heights of ! ! rho levels above z=0 REAL(KIND=r8) :: r_theta_levels_land(nCols,kMax)! land field of heights of ! ! theta levels above z=0 REAL(KIND=r8) :: rho_land (nCols,kMax) ! density at land points REAL(KIND=r8) :: s_x_lin_stress (nCols) ! 'surface' x_lin_stress land pnts REAL(KIND=r8) :: s_y_lin_stress (nCols) ! 'surface' y_lin_stress land pnts REAL(KIND=r8) :: s_x_wake_stress (nCols) ! 'surface' x_wake_stress land pts REAL(KIND=r8) :: s_y_wake_stress (nCols) ! 'surface' y_wake_stress land pts REAL(KIND=r8) :: s_x_orog (nCols) ! 'surface' x_orog on land points REAL(KIND=r8) :: s_y_orog (nCols) ! 'surface' y_orog on land points REAL(KIND=r8) :: lift (nCols) ! depth of blocked layer REAL(KIND=r8) :: fr (nCols) ! low level froude number REAL(KIND=r8) :: rho_s (nCols) ! low level density REAL(KIND=r8) :: orog_grad_xx_land(nCols)! dh/dx squared gradient orography REAL(KIND=r8) :: orog_grad_yy_land(nCols)! dh/dy squared gradient orography REAL(KIND=r8) :: orog_grad_xy_land(nCols)! (dh/dx)(dh/dy) gradient orography REAL(KIND=r8) :: sd_orog_land (nCols)! standard deviation of orography (m) ! ! Land points arrays below are for the total GWD stress and ! its 4 individual components. (x and y components for each) ! REAL(KIND=r8) :: stress_ud_land ( nCols , 0:kMax ) REAL(KIND=r8) :: stress_vd_land ( nCols , 0:kMax ) REAL(KIND=r8) :: stress_ud_satn_land ( nCols , 0:kMax ) REAL(KIND=r8) :: stress_vd_satn_land ( nCols , 0:kMax ) REAL(KIND=r8) :: stress_ud_wake_land ( nCols , 0:kMax ) REAL(KIND=r8) :: stress_vd_wake_land ( nCols , 0:kMax ) ! ! Land point arrays below are for the 4 individual components of ! the GWD wind increment. (x and y components for each) ! REAL(KIND=r8) :: du_dt_satn_land ( nCols , kMax ) REAL(KIND=r8) :: dv_dt_satn_land ( nCols , kMax ) REAL(KIND=r8) :: du_dt_wake_land ( nCols , kMax ) REAL(KIND=r8) :: dv_dt_wake_land ( nCols , kMax ) ! ! Land point arrays below are for the 9 GWD 'surface' diagnostics. ! REAL(KIND=r8) :: u_s_d_land ( nCols ) REAL(KIND=r8) :: v_s_d_land ( nCols ) REAL(KIND=r8) :: nsq_s_d_land ( nCols ) REAL(KIND=r8) :: fr_d_land ( nCols ) REAL(KIND=r8) :: bld_d_land ( nCols ) REAL(KIND=r8) :: bldt_d_land ( nCols ) REAL(KIND=r8) :: num_lim_d_land ( nCols ) REAL(KIND=r8) :: num_fac_d_land ( nCols ) REAL(KIND=r8) :: tausx_d_land ( nCols ) REAL(KIND=r8) :: tausy_d_land ( nCols ) REAL(KIND=r8) :: taus_scale_d_land ( nCols ) LOGICAL :: l_drag(nCols) ! whether point has a non-zero stress or not REAL(KIND=r8) :: pmid (nCols,kMax) REAL(KIND=r8) :: state_t (nCols,kMax) REAL(KIND=r8) :: state_q (nCols,kMax) REAL(KIND=r8) :: state_u (nCols,kMax) REAL(KIND=r8) :: state_v (nCols,kMax) REAL(KIND=r8) :: state_pmid (nCols,kMax) REAL(KIND=r8) :: state_pint (nCols,kMax+1) REAL(KIND=r8) :: state_lnpint(nCols,kMax+1) REAL(KIND=r8) :: state_pdel (nCols,kMax) REAL(KIND=r8) :: state_rpdel (nCols,kMax) REAL(KIND=r8) :: state_lnpmid(nCols,kMax) REAL(KIND=r8) :: state_zi (1:nCols,1:kMax+1) REAL(KIND=r8) :: state_zm (1:nCols,1:kMax) k_top=0 k_top_max=0 theta =0.0_r8 rho =0.0_r8 r_rho_levels =0.0_r8 r_theta_levels=0.0_r8 stress_ud (1:nCols ,0:kMax)=0.0_r8 stress_vd (1:nCols ,0:kMax)=0.0_r8 stress_ud_satn(1:nCols ,0:kMax)=0.0_r8 stress_vd_satn(1:nCols ,0:kMax)=0.0_r8 stress_ud_wake(1:nCols ,0:kMax)=0.0_r8 stress_vd_wake(1:nCols ,0:kMax)=0.0_r8 du_dt_satn =0.0_r8 dv_dt_satn =0.0_r8 du_dt_land =0.0_r8 dv_dt_land =0.0_r8 du_dt_wake =0.0_r8 dv_dt_wake =0.0_r8 u_s_d =0.0_r8 v_s_d =0.0_r8 nsq_s_d =0.0_r8 fr_d =0.0_r8 bld_d =0.0_r8 bldt_d =0.0_r8 num_lim_d =0.0_r8 num_fac_d =0.0_r8 tausx_d =0.0_r8 tausy_d =0.0_r8 taus_scale_d =0.0_r8 r_theta_levels_land=0.0_r8 s_x_lin_stress =0.0_r8 s_y_lin_stress =0.0_r8 s_x_wake_stress =0.0_r8 s_y_wake_stress =0.0_r8 s_x_orog=0.0_r8 s_y_orog=0.0_r8 lift=0.0_r8 fr=0.0_r8 rho_s=0.0_r8 orog_grad_xx_land=0.0_r8 orog_grad_yy_land=0.0_r8 orog_grad_xy_land=0.0_r8 sd_orog_land=0.0_r8 stress_ud_land (1: nCols , 0:kMax )= 0.0_r8 stress_vd_land (1: nCols , 0:kMax )= 0.0_r8 stress_ud_satn_land(1: nCols , 0:kMax )= 0.0_r8 stress_vd_satn_land(1: nCols , 0:kMax )= 0.0_r8 stress_ud_wake_land(1: nCols , 0:kMax )= 0.0_r8 stress_vd_wake_land(1: nCols , 0:kMax )= 0.0_r8 du_dt_satn_land=0.0_r8 dv_dt_satn_land=0.0_r8 du_dt_wake_land=0.0_r8 dv_dt_wake_land=0.0_r8 taus_scale_d_land=0.0_r8 land_points=0 DO i=1,ncols IF(imask(i).GE.1_i8) THEN land_points=land_points+1 END IF END DO ! ! END of full field diagnostic arrays ! Below are the stash flags for calculating diagnostics: ! stress_ud_on=.TRUE. !u stress points_stress_ud=land_points stress_ud (1:land_points ,0:kMax)=0.0_r8 !u total stress stress_ud_p_on=.TRUE. !u on press stress stress_vd_on=.TRUE. !v stress points_stress_vd=land_points stress_vd (1:land_points ,0:kMax)=0.0_r8 !v total stress stress_ud_satn_on=.TRUE. !u satn stress points_stress_ud_satn=land_points stress_ud_satn(1:land_points ,0:kMax)=0.0_r8!u satn stress stress_vd_satn_on=.TRUE. !v satn stress points_stress_vd_satn=land_points stress_vd_satn(1:land_points ,0:kMax)=0.0_r8!v satn stress stress_ud_wake_on=.TRUE. !u wake stress points_stress_ud_wake=land_points stress_ud_wake(1:land_points ,0:kMax)=0.0_r8!u wake stress stress_vd_wake_on=.TRUE. !v wake stress points_stress_vd_wake=land_points stress_vd_wake(1:land_points ,0:kMax)=0.0_r8!v wake stress du_dt_satn_on=.TRUE. !u accel (saturation) points_du_dt_satn=land_points du_dt_satn (1:land_points ,1:kMax)=0.0_r8!u acceln (saturation) du_dt_satn_p_on=.TRUE. !u accel (saturation) dv_dt_satn_on=.TRUE. !v accel (saturation) points_dv_dt_satn=land_points dv_dt_satn (1:land_points ,1:kMax)=0.0_r8!v acceln (saturation) du_dt_wake_on=.TRUE. !u accel blocked flow points_du_dt_wake=land_points du_dt_wake (1:land_points ,1:kMax)=0.0_r8!u acceln (blocked flow) dv_dt_wake_on=.TRUE. !v accel blocked flow points_dv_dt_wake=land_points dv_dt_wake (1:land_points ,1:kMax)=0.0_r8!v acceln (blocked flow) u_s_d_on=.TRUE. !u_s_d diag switch points_u_s_d=land_points u_s_d (1:land_points )=0.0_r8 ! u_s diag at theta pts v_s_d_on=.TRUE. !v_s_d diag switch points_v_s_d=land_points v_s_d (1:land_points ) =0.0_r8! v_s diag at theta pts nsq_s_d_on=.TRUE. !nsq_s_d diag switch points_nsq_s_d=land_points nsq_s_d (1:land_points ) =0.0_r8! nsq_s diag at theta pts fr_d_on=.TRUE. !fr_d switch points_fr_d=land_points fr_d (1:land_points ) =0.0_r8 ! Fr diag at theta pts bld_d_on=.TRUE. !bld_d switch points_bld_d=land_points bld_d (1:land_points )=0.0_r8 ! blocked layer depth at theta pts bldt_d_on=.TRUE. !bldt_d switch points_bldt_d=land_points bldt_d (1:land_points )=0.0_r8 ! % of time blocked layer diagnosed num_lim_d_on=.TRUE. !num_lim_d switch points_num_lim_d=land_points num_lim_d (1:land_points )=0.0_r8 ! % of time numerical limiter invoked num_fac_d_on=.TRUE. !num_fac_d switch points_num_fac_d=land_points num_fac_d (1:land_points )=0.0_r8 ! % redn. of flow-blocking stress ! after numerical limiter invoked tausx_d_on=.TRUE. !tausx_d switch points_tausx_d=land_points tausx_d (1:land_points)=0.0_r8!x-component of surface stress tausy_d_on=.TRUE. !tausy_d switch points_tausy_d=land_points tausy_d (1:land_points)=0.0_r8!y-component of surface stress taus_scale_d_on=.TRUE. !taus_scale_d switch points_taus_scale_d=land_points taus_scale_d(1:land_points)=0.0_r8! Factor surface stress scaled by ! if Froude no. dependence is on. !----------------------------------------------------------------------------- DO i=1,nCols state_pint (i,kMax+1) = gps(i)*si(1) END DO DO k=kMax,1,-1 DO i=1,nCols state_pint (i,k) = MAX(si(kMax+2-k)*gps(i) ,0.0001_r8) END DO END DO DO k=1,kMax DO i=1,nCols state_t (i,kMax+1-k) = gt (i,k) state_q (i,kMax+1-k) = gq (i,k) state_u (i,kMax+1-k) = gu (i,k) state_v (i,kMax+1-k) = gv (i,k) state_pmid(i,kMax+1-k) = sl(k)*gps (i) END DO END DO DO k=1,kMax DO i=1,nCols state_pdel (i,k) = MAX(state_pint(i,k+1) - state_pint(i,k),0.000000005_r8) state_rpdel (i,k) = 1.0_r8/MAX((state_pint(i,k+1) - state_pint(i,k)),0.00000000005_r8) state_lnpmid (i,k) = LOG(state_pmid(i,k)) END DO END DO DO k=1,kMax+1 DO i=1,nCols state_lnpint(i,k) = LOG(state_pint (i,k)) END DO END DO ! !..delsig k=2 ****gu,gv,gt,gq,gyu,gyv,gtd,gqd,sl*** } delsig(2) ! k=3/2----si,ric,rf,km,kh,b,l ----------- ! k=1 ****gu,gv,gt,gq,gyu,gyv,gtd,gqd,sl*** } delsig(1) ! k=1/2----si ---------------------------- ! Derive new temperature and geopotential fields CALL geopotential_t( & state_lnpint(1:nCols,1:kMax+1) , state_pint (1:nCols,1:kMax+1) , & state_pmid (1:nCols,1:kMax) , state_pdel (1:nCols,1:kMax) , state_rpdel(1:nCols,1:kMax) , & state_t (1:nCols,1:kMax) , state_q (1:nCols,1:kMax) , rair , gravit , zvir ,& state_zi (1:nCols,1:kMax+1) , state_zm (1:nCols,1:kMax) , nCols ,nCols, kMax) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! Parei aqui ! DO i=1,nCols R_THETA_LEVELS(i,0) = EARTH_RADIUS + MAX(topo(i),0.0_r8)+ state_zi(i,kMax+1) END DO DO k=1,kMax DO i=1,nCols work_u (i,k) = gu (i,k) work_v (i,k) = gv (i,k) R_THETA_LEVELS (i,k) = (EARTH_RADIUS + MAX(topo(i),0.0_r8)) + state_zi (i,kMax+1-k) r_rho_levels (i,k) = (EARTH_RADIUS + MAX(topo(i),0.0_r8)) + state_zm (i,kMax+1-k) THETA (i,k) = sigkiv(k)*gt(i,k) pmid (i,k) = gps(i)*sl(k) !PRINT*,R_THETA_LEVELS (i,k),state_zi (i,kMax+1-k) !j/kg/kelvin ! ! P = rho * R * T ! ! P ! rho = ------- ! R * T ! ! 1 R * T ! rrho = ----- = ------- ! rho P ! ! rho(i,k) = (pmid(i,k)/(gasr*gt(i,k)))*(earth_radius*earth_radius)! density*(radius of earth)^2 END DO END DO !------------------------------------------------------------------ !l 1.2 gather winds at land points !------------------------------------------------------------------ DO k=1,kMax land_points=0 DO i=1,ncols IF(imask(i).GE.1_i8) THEN land_points=land_points+1 !------------------------------------------------------------------ !l 1.2 gather winds at land points !------------------------------------------------------------------ up_land(land_points,k) =work_u(i,k) vp_land(land_points,k) =work_v(i,k) !------------------------------------------------------------------ !l 1.3 gather theta, rho and heights at land points !------------------------------------------------------------------ r_rho_levels_land (land_points,k)= r_rho_levels (i,k) - r_theta_levels(i,0) r_theta_levels_land(land_points,k)= r_theta_levels(i,k) - r_theta_levels(i,0) rho_land (land_points,k)= rho (i,k)*recip_a2 theta_land (land_points,k)= theta(i,k) END IF END DO END DO land_points=0 DO i=1,ncols IF(imask(i).GE.1_i8) THEN land_points=land_points+1 orog_grad_xx_land(land_points) = orog_grad_xx(i) orog_grad_yy_land(land_points) = orog_grad_yy(i) orog_grad_xy_land(land_points) = orog_grad_xy(i) sd_orog_land (land_points) = sd_orog (i) END IF END DO !------------------------------------------------------------------ !l 2. calculate anisotropic 'surface' stress,CALL gw_surf !------------------------------------------------------------------ ! PRINT*,'first gw_surf ' IF (land_points > 0) THEN points=land_points ! DEPENDS ON: gw_surf CALL gw_surf(& r_theta_levels_land(1:points,1:kMax), &!REAL,intent(in):: r_theta_levels (points,kMax)! heights on theta levels rho_land (1:points,1:kMax), &!REAL,intent(in):: rho (points,kMax)! density on rho levels theta_land (1:points,1:kMax), &!REAL,intent(in):: theta (points,kMax)! theta on theta levels up_land (1:points,1:kMax), &!REAL,intent(in):: u (points,kMax)! u on rho levels vp_land (1:points,1:kMax), &!REAL,intent(in):: v (points,kMax)! v on rho levels timestep , &!REAL,intent(in):: timestep ! timestep (s) sd_orog_land (1:points) , &!REAL,intent(in):: sd_orog (points)! standard deviation of the sub-grid orography orog_grad_xx_land (1:points) , &!REAL,intent(in):: sigma_xx (points)! (dh/dx)^2 grid box average of the orog_grad_xy_land (1:points) , &!REAL,intent(in):: sigma_xy (points)! (dh/dx)*(dh/dy) orog_grad_yy_land (1:points) , &!REAL,intent(in):: sigma_yy (points)! (dh/dy)^2 s_x_lin_stress (1:points) , &!REAL,intent(out):: s_x_lin_stress (points)! 'surface' linear stress in x-dirn s_y_lin_stress (1:points) , &!REAL,intent(out):: s_y_lin_stress (points)! 'surface' linear stress in y-dirn s_x_wake_stress (1:points) , &!REAL,intent(out):: s_x_wake_stress(points)! wake stress in x-dirn s_y_wake_stress (1:points) , &!REAL,intent(out):: s_y_wake_stress(points)! wake stress in y-dirn s_x_orog (1:points) , &!REAL,intent(out):: s_x_orog (points)! 'surface' stress/orog in x-dirn s_y_orog (1:points) , &!REAL,intent(out):: s_y_orog (points)! 'surface' stress/orog in y-dirn kMax , &!INTEGER,intent(in) :: kMax ! number of model levels land_points , &!INTEGER,intent(in) :: points ! number of gwd points on pe kay , &!REAL ,intent(in):: kay ! gwd surface stress constant (m-1) rho_s (1:points) , &!REAL,intent(out):: rho_s(points) ! density - av from z=0 to nsigma*sd_orog l_taus_scale , &!LOGICAL,intent(in):: l_taus_scale ! true allows variation of surface stress on the low level Froude number k_top (1:points) , &!INTEGER ,intent(out):: k_top(points) ! topmost model level including mountains. k_top_max , &!INTEGER ,intent(out):: k_top_max ! max(k_top) lift (1:points) , &! REAL,intent(out):: lift(points) ! blocked layer depth l_drag (1:points) , &!LOGICAL,intent(out):: l_drag(points) ! whether a point has a non-zero surface stress or not fr (1:points) , &!REAL,intent(out):: fr(points) ! low level froude number for linear frc , &!REAL ,intent(in):: frc ! critical froude number below which hydraulic jumps are triggered u_s_d_land (1:points) , &!REAL,intent(out):: u_s_d(points_u_s_d) ! 0-nsigma*sd_orog u_s diag u_s_d_on , &!LOGICAL,intent(in):: u_s_d_on points_u_s_d , &!INTEGER ,intent(in) ::points_u_s_d v_s_d_land (1:points) , &!REAL,intent(out):: v_s_d(points_v_s_d) ! 0-nsigma_sd_orog v_s diag v_s_d_on , &!LOGICAL,intent(in):: v_s_d_on points_v_s_d , &!INTEGER ,intent(in) ::points_v_s_d nsq_s_d_land (1:points) , &!REAL,intent(out):: nsq_s_d(points_nsq_s_d)! 0-nsigma*sd_orog nsq_s diagnostic nsq_s_d_on , &!LOGICAL,intent(in):: nsq_s_d_on points_nsq_s_d , &!INTEGER ,intent(in) ::points_nsq_s_d fr_d_land (1:points) , &!REAL,intent(out):: fr_d(points_fr_d) ! Froude no. diagnostic fr_d_on , &!LOGICAL,intent(in):: fr_d_on points_fr_d , &!INTEGER ,intent(in) ::points_fr_d bld_d_land (1:points) , &!REAL,intent(out):: bld_d(points_bld_d) ! blocked layer depth diagnostic bld_d_on , &!LOGICAL,intent(in):: bld_d_on points_bld_d , &!INTEGER ,intent(in) ::points_bld_d bldt_d_land (1:points) , &!REAL,intent(out):: bldt_d(points_bldt_d) ! % of time blocked flow param. invoked bldt_d_on , &!LOGICAL,intent(in):: bldt_d_on points_bldt_d , &!INTEGER ,intent(in) ::points_bldt_d num_lim_d_land (1:points) , &!REAL,intent(out):: num_lim_d(points_num_lim_d)! % of time numerical limiter invoked num_lim_d_on , &!LOGICAL,intent(in):: num_lim_d_on points_num_lim_d , &!INTEGER ,intent(in) ::points_num_lim_d num_fac_d_land (1:points) , &!REAL,intent(out):: num_fac_d(points_num_fac_d)! % reduction of flow blocking stressafter numerical limiter invoked num_fac_d_on , &!LOGICAL,intent(in):: num_fac_d_on points_num_fac_d , &!INTEGER ,intent(in) ::points_num_fac_d tausx_d_land (1:points) , &!REAL,intent(out):: tausx_d(points_tausx_d) ! x-component of total surface stress tausx_d_on , &!LOGICAL,intent(in):: tausx_d_on points_tausx_d , &!INTEGER ,intent(in) ::points_tausx_d tausy_d_land (1:points) , &!REAL,intent(out):: tausy_d(points_tausy_d) ! y-component of total surface stress tausy_d_on , &!LOGICAL,intent(in):: tausy_d_on points_tausy_d , &!INTEGER ,intent(in) ::points_tausy_d taus_scale_d_land (1:points) , &!REAL,intent(out):: taus_scale_d(points_taus_scale_d)! scaling factor for surface stress when Fr dependence of surface stress invoked. taus_scale_d_on , &!LOGICAL,intent(in):: taus_scale_d_on points_taus_scale_d )!INTEGER ,intent(in) ::points_taus_scale_d !PRINT*,'end gw_surf ' !------------------------------------------------------------------ !l 3. calculate stress profile and accelerations, !l CALL gw_vert !------------------------------------------------------------------ !PRINT*,'first gw_vert ' ! DEPENDS ON: gw_vert CALL gw_vert( & rho_land (1:points,1:kMax), &!REAL,intent(in):: rho (points,kMax)! density on rho levels r_rho_levels_land (1:points,1:kMax), &!REAL,intent(in):: r_rho_levels (points,kMax)! heights on rho levels r_theta_levels_land (1:points,1:kMax), &!REAL,intent(in):: r_theta_levels(points,kMax)! heights on theta levels theta_land (1:points,1:kMax), &!REAL,intent(in):: theta (points,kMax)! theta field up_land (1:points,1:kMax), &!REAL,intent(in):: u (points,kMax)! u field vp_land (1:points,1:kMax), &!REAL,intent(in):: v (points,kMax)! v field kMax , &!INTEGER ,intent(in) :: kMax ! number of model levels land_points , &!INTEGER ,intent(in) :: points ! number of points kay , &!REAL,intent(in):: kay ! GWD hydrostatic constant sd_orog_land (1:points) , &!REAL,intent(in):: sd_orog (points) ! standard deviation of the sub-grid orography s_x_lin_stress (1:points) , &!REAL,intent(in):: s_x_lin_stress (points) ! 'surface' lin x_stress s_y_lin_stress (1:points) , &!REAL,intent(in):: s_y_lin_stress (points)! 'surface' lin y_stress s_x_wake_stress (1:points) , &!REAL,intent(in):: s_x_wake_stress(points)! 'surface' x_wake_stress s_y_wake_stress (1:points) , &!REAL,intent(in):: s_y_wake_stress(points)! 'surface' y_wake_stress s_x_orog (1:points) , &!REAL,intent(in):: s_x_orog (points) ! 'surface' x_stress term s_y_orog (1:points) , &!REAL,intent(in):: s_y_orog (points) ! 'surface' y_stress term du_dt_land (1:points,1:kMax), &!REAL,intent(out):: du_dt (points,kMax) ! total GWD du/dt dv_dt_land (1:points,1:kMax), &!REAL,intent(out):: dv_dt (points,kMax) ! total GWD dv/dt k_top (1:points) , &!INTEGER,intent(in):: k_top (points) ! model level at mountain tops k_top_max , &!INTEGER,intent(in):: k_top_max! max(k_top) !lift (1:points) , &!REAL ,intent(in):: lift (points) ! blocked layer depth l_drag (1:points) , &!LOGICAL,intent(in):: l_drag (points)! true if a non-zero surface stress fr (1:points) , &!REAL ,intent(in):: fr (points)! low level froude number rho_s (1:points) , &!REAL ,intent(in):: rho_s (points)! Low level density calculated in gwsurf l_fix_gwsatn , &!LOGICAL,intent(in):: l_fix_gwsatn ! switch to include minor bug fixes l_gwd_40km , &!LOGICAL,intent(in):: l_gwd_40km ! switch to turn off GWD above 40km ! diagnostics stress_ud_land (1:points,0:kMax), &!REAL ,intent(out) :: stress_ud (points_stress_ud,0:kMax)! x total stress diag points_stress_ud , &!INTEGER,intent(in ) :: points_stress_ud stress_ud_on , &!LOGICAL,intent(in ) :: stress_ud_on stress_ud_p_on , &!LOGICAL,intent(in ) :: stress_ud_p_on stress_vd_land (1:points,0:kMax), &!REAL ,intent(out) :: stress_vd (points_stress_vd,0:kMax)! y total stress diag points_stress_vd , &!INTEGER,intent(in ) :: points_stress_vd stress_vd_on , &!LOGICAL,intent(in ) :: stress_vd_on stress_ud_satn_land (1:points,0:kMax), &!REAL ,intent(out) :: stress_ud_satn (points_stress_ud_satn,0:kMax)! x saturation stress diag points_stress_ud_satn , &!INTEGER,intent(in ) :: points_stress_ud_satn stress_ud_satn_on , &!LOGICAL,intent(in ) :: stress_ud_satn_on stress_vd_satn_land (1:points,0:kMax), &!REAL ,intent(out) :: stress_vd_satn (points_stress_vd_satn,0:kMax)! y saturation stress diag points_stress_vd_satn , &!INTEGER,intent(in ) :: points_stress_vd_satn stress_vd_satn_on , &!LOGICAL,intent(in ) :: stress_vd_satn_on stress_ud_wake_land (1:points,0:kMax), &!REAL ,intent(out) :: stress_ud_wake (points_stress_ud_wake,0:kMax)! x blocked flow stress diag points_stress_ud_wake , &!INTEGER,intent(in ) :: points_stress_ud_wake stress_ud_wake_on , &!LOGICAL,intent(in ) :: stress_ud_wake_on stress_vd_wake_land (1:points,0:kMax), &!REAL ,intent(out) :: stress_vd_wake (points_stress_vd_wake,0:kMax)! y blocked flow stress diag points_stress_vd_wake , &!INTEGER,intent(in ) :: points_stress_vd_wake stress_vd_wake_on , &!LOGICAL,intent(in ) :: stress_vd_wake_on du_dt_satn_land (1:points,1:kMax), &!REAL ,intent(out) :: du_dt_satn (points_du_dt_satn,kMax)! du/dt diagnostic (saturation) points_du_dt_satn , &!INTEGER,intent(in ) :: points_du_dt_satn du_dt_satn_on , &!LOGICAL,intent(in ) :: du_dt_satn_on du_dt_satn_p_on , &!LOGICAL,intent(in ) :: du_dt_satn_p_on dv_dt_satn_land (1:points,1:kMax), &!REAL ,intent(out) :: dv_dt_satn (points_dv_dt_satn,kMax)! dv/dt diagnostic (saturation) points_dv_dt_satn , &!INTEGER,intent(in ) :: points_dv_dt_satn dv_dt_satn_on , &!LOGICAL,intent(in ) :: dv_dt_satn_on du_dt_wake_land (1:points,1:kMax), &!REAL ,intent(out) :: du_dt_wake (points_du_dt_wake,kMax)! du/dt diagnostic (blocked flow) points_du_dt_wake , &!INTEGER,intent(in ) :: points_du_dt_wake du_dt_wake_on , &!LOGICAL,intent(in ) :: du_dt_wake_on dv_dt_wake_land (1:points,1:kMax), &!REAL ,intent(out) :: dv_dt_wake (points_dv_dt_wake,kMax)! dv/dt diagnostic (blocked flow) points_dv_dt_wake , &!INTEGER,intent(in ) :: points_dv_dt_wake dv_dt_wake_on )!LOGICAL,intent(in ) :: dv_dt_wake_on !PRINT*,'end gw_vert ' DO k=1,kMax land_points=0 DO i=1,ncols IF(imask(i).GE.1_i8) THEN land_points=land_points+1 !------------------------------------------------------------------ !l 1.2 gather winds at land points !------------------------------------------------------------------ du_dt(i,k) =du_dt_land(land_points,k) dv_dt(i,k) =dv_dt_land(land_points,k) !PRINT*,k,du_dt(i,k),dv_dt(i,k) END IF END DO END DO END IF ! on land_points > 0 iret=0 END SUBROUTINE g_wave ! *****************************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 gw_vert to calculate the vertical profile of gravity wave ! stress and associated wind increments ! SUBROUTINE gw_vert( & rho , &!REAL,intent(in):: rho(points,kMax)! density on rho levels r_rho_levels , &!REAL,intent(in):: r_rho_levels(points,kMax)! heights on rho levels r_theta_levels , &!REAL,intent(in):: r_theta_levels(points,kMax)! heights on theta levels theta , &!REAL,intent(in):: theta(points,kMax) ! theta field u , &!REAL,intent(in):: u(points,kMax)! u field v , &!REAL,intent(in):: v(points,kMax)! v field kMax , &!INTEGER ,intent(in) :: kMax ! number of model levels points , &!INTEGER ,intent(in) :: points ! number of points kay , &!REAL,intent(in):: kay ! GWD hydrostatic constant sd_orog , &!REAL,intent(in):: sd_orog(points) ! standard deviation of the sub-grid orography s_x_lin_stress , &!REAL,intent(in):: s_x_lin_stress(points) ! 'surface' lin x_stress s_y_lin_stress , &!REAL,intent(in):: s_y_lin_stress(points)! 'surface' lin y_stress s_x_wake_stress , &!REAL,intent(in):: s_x_wake_stress(points)! 'surface' x_wake_stress s_y_wake_stress , &!REAL,intent(in):: s_y_wake_stress(points)! 'surface' y_wake_stress s_x_orog , &!REAL,intent(in):: s_x_orog(points) ! 'surface' x_stress term s_y_orog , &!REAL,intent(in):: s_y_orog(points) ! 'surface' y_stress term du_dt , &!REAL,intent(out):: du_dt(points,kMax) ! total GWD du/dt dv_dt , &!REAL,intent(out):: dv_dt(points,kMax) ! total GWD dv/dt k_top , &!INTEGER ,intent(in) :: k_top(points) ! model level at mountain tops k_top_max , &!INTEGER ,intent(in) :: k_top_max ! max(k_top) !lift , &!REAL,intent(in):: lift (points) ! blocked layer depth l_drag , &!LOGICAL,intent(in):: l_drag(points) ! true if a non-zero surface stress fr , &!REAL,intent(in):: fr(points) ! low level froude number rho_s , &!REAL,intent(in):: rho_s(points) ! Low level density calculated in gwsurf l_fix_gwsatn , &!LOGICAL,intent(in):: l_fix_gwsatn ! switch to include minor bug fixes l_gwd_40km , &!LOGICAL,intent(in):: l_gwd_40km ! switch to turn off GWD above 40km ! diagnostics stress_ud , &!REAL,intent(out):: stress_ud (points_stress_ud,0:kMax)! x total stress diag points_stress_ud , &!INTEGER,intent(in):: points_stress_ud stress_ud_on , &!LOGICAL,intent(in) :: stress_ud_on stress_ud_p_on , &!LOGICAL,intent(in) :: stress_ud_p_on stress_vd , &!REAL,intent(out):: stress_vd (points_stress_vd,0:kMax)! y total stress diag points_stress_vd , &!INTEGER,intent(in):: points_stress_vd stress_vd_on , &!LOGICAL,intent(in) :: stress_vd_on stress_ud_satn , &!REAL,intent(out):: stress_ud_satn(points_stress_ud_satn,0:kMax)! x saturation stress diag points_stress_ud_satn , &!INTEGER,intent(in):: points_stress_ud_satn stress_ud_satn_on , &!LOGICAL,intent(in) :: stress_ud_satn_on stress_vd_satn , &!REAL,intent(out):: stress_vd_satn(points_stress_vd_satn,0:kMax)! y saturation stress diag points_stress_vd_satn , &!INTEGER,intent(in):: points_stress_vd_satn stress_vd_satn_on , &!LOGICAL,intent(in) :: stress_vd_satn_on stress_ud_wake , &!REAL,intent(out):: stress_ud_wake(points_stress_ud_wake,0:kMax)! x blocked flow stress diag points_stress_ud_wake , &!INTEGER,intent(in):: points_stress_ud_wake stress_ud_wake_on , &!LOGICAL,intent(in) :: stress_ud_wake_on stress_vd_wake , &!REAL,intent(out):: stress_vd_wake(points_stress_vd_wake,0:kMax)! y blocked flow stress diag points_stress_vd_wake , &!INTEGER,intent(in):: points_stress_vd_wake stress_vd_wake_on , &!LOGICAL,intent(in) :: stress_vd_wake_on du_dt_satn , &!REAL,intent(out):: du_dt_satn(points_du_dt_satn,kMax)! du/dt diagnostic (saturation) points_du_dt_satn , &!INTEGER,intent(in):: points_du_dt_satn du_dt_satn_on , &!LOGICAL,intent(in) :: du_dt_satn_on du_dt_satn_p_on , &!LOGICAL,intent(in) :: du_dt_satn_p_on dv_dt_satn , &!REAL,intent(out):: dv_dt_satn(points_dv_dt_satn,kMax)! dv/dt diagnostic (saturation) points_dv_dt_satn , &!INTEGER,intent(in):: points_dv_dt_satn dv_dt_satn_on , &!LOGICAL,intent(in) :: dv_dt_satn_on du_dt_wake , &!REAL,intent(out):: du_dt_wake(points_du_dt_wake,kMax)! du/dt diagnostic (blocked flow) points_du_dt_wake , &!INTEGER,intent(in):: points_du_dt_wake du_dt_wake_on , &!LOGICAL,intent(in) :: du_dt_wake_on dv_dt_wake , &!REAL,intent(out):: dv_dt_wake(points_dv_dt_wake,kMax)! dv/dt diagnostic (blocked flow) points_dv_dt_wake , &!INTEGER,intent(in):: points_dv_dt_wake dv_dt_wake_on )!LOGICAL,intent(in) :: dv_dt_wake_on IMPLICIT NONE ! Description: ! calculates the gwd stress profiles and wind increments. ! 1. calculate stress profile and wind increments for ! linear hydrostatic waves. ! 2. calculate stress profile and wind increments for ! the blocked flow. ! ! current code owner: S. Webster ! ! history: ! version date comment ! 5.2 15/11/00 original code. Stuart Webster. ! ! 5.3 16/10/01 Remove lee waves and hydraulic jumps as part ! of the simpler 4A scheme. Stuart Webster ! ! 6.2 21/02/06 Introduce minor bug fixes using l_fix_gwsatn. ! Stuart Webster ! ! code description: ! language: fortran 77 + common extensions ! this code is written to umdp3 v6 programming standards. ! suitable for single column use,rotated grids ! global variables ! local constants ! none ! ! SUBROUTINE ARGUMENTS: ! INTEGER ,INTENT(in) :: kMax ! number of model levels INTEGER ,INTENT(in) :: points ! number of points INTEGER ,INTENT(in) :: k_top(points) ! model level at mountain tops ! ! full definition in gwsurf INTEGER ,INTENT(in) :: k_top_max ! max(k_top) ! ! Integers below determine size of diagnostic arrays. They are ! set to points if called and to 1 if not. This is done in GWD_CTL2 ! INTEGER,INTENT(in):: points_stress_ud INTEGER,INTENT(in):: points_stress_vd INTEGER,INTENT(in):: points_stress_ud_satn INTEGER,INTENT(in):: points_stress_vd_satn INTEGER,INTENT(in):: points_stress_ud_wake INTEGER,INTENT(in):: points_stress_vd_wake INTEGER,INTENT(in):: points_du_dt_satn INTEGER,INTENT(in):: points_dv_dt_satn INTEGER,INTENT(in):: points_du_dt_wake INTEGER,INTENT(in):: points_dv_dt_wake REAL(KIND=r8),INTENT(in):: r_rho_levels(points,kMax)! heights on rho levels REAL(KIND=r8),INTENT(in):: r_theta_levels(points,kMax)! heights on theta levels REAL(KIND=r8),INTENT(in):: rho(points,kMax)! density on rho levels REAL(KIND=r8),INTENT(in):: theta(points,kMax) ! theta field REAL(KIND=r8),INTENT(in):: u(points,kMax)! u field REAL(KIND=r8),INTENT(in):: v(points,kMax)! v field REAL(KIND=r8),INTENT(in):: s_x_lin_stress(points) ! 'surface' lin x_stress REAL(KIND=r8),INTENT(in):: s_y_lin_stress(points)! 'surface' lin y_stress REAL(KIND=r8),INTENT(in):: s_x_wake_stress(points)! 'surface' x_wake_stress REAL(KIND=r8),INTENT(in):: s_y_wake_stress(points)! 'surface' y_wake_stress REAL(KIND=r8),INTENT(in):: s_x_orog(points) ! 'surface' x_stress term REAL(KIND=r8),INTENT(in):: s_y_orog(points) ! 'surface' y_stress term !REAL(KIND=r8),INTENT(in):: lift (points) ! blocked layer depth REAL(KIND=r8),INTENT(in):: sd_orog(points) ! standard deviation of the sub-grid orography REAL(KIND=r8),INTENT(in):: fr(points) ! low level froude number REAL(KIND=r8),INTENT(in):: kay ! GWD hydrostatic constant REAL(KIND=r8),INTENT(in):: rho_s(points) ! Low level density calculated in gwsurf REAL(KIND=r8),INTENT(out):: du_dt(points,kMax) ! total GWD du/dt REAL(KIND=r8),INTENT(out):: dv_dt(points,kMax) ! total GWD dv/dt ! ! Diagnostics ! REAL(KIND=r8),INTENT(out):: stress_ud (points_stress_ud,0:kMax)! x total stress diag REAL(KIND=r8),INTENT(out):: stress_vd (points_stress_vd,0:kMax)! y total stress diag REAL(KIND=r8),INTENT(out):: stress_ud_satn(points_stress_ud_satn,0:kMax)! x saturation stress diag REAL(KIND=r8),INTENT(out):: stress_vd_satn(points_stress_vd_satn,0:kMax)! y saturation stress diag REAL(KIND=r8),INTENT(out):: stress_ud_wake(points_stress_ud_wake,0:kMax)! x blocked flow stress diag REAL(KIND=r8),INTENT(out):: stress_vd_wake(points_stress_vd_wake,0:kMax)! y blocked flow stress diag REAL(KIND=r8),INTENT(out):: du_dt_satn(points_du_dt_satn,kMax)! du/dt diagnostic (saturation) REAL(KIND=r8),INTENT(out):: dv_dt_satn(points_dv_dt_satn,kMax)! dv/dt diagnostic (saturation) REAL(KIND=r8),INTENT(out):: du_dt_wake(points_du_dt_wake,kMax)! du/dt diagnostic (blocked flow) REAL(KIND=r8),INTENT(out):: dv_dt_wake(points_dv_dt_wake,kMax)! dv/dt diagnostic (blocked flow) LOGICAL,INTENT(in):: l_drag(points) ! true if a non-zero surface stress LOGICAL,INTENT(in):: l_fix_gwsatn ! switch to include minor bug fixes LOGICAL,INTENT(in):: l_gwd_40km ! switch to turn off GWD above 40km ! ! Diagnostic switches ! LOGICAL,INTENT(in) :: stress_ud_on LOGICAL,INTENT(in) :: stress_ud_p_on LOGICAL,INTENT(in) :: stress_vd_on LOGICAL,INTENT(in) :: stress_ud_satn_on LOGICAL,INTENT(in) :: stress_vd_satn_on LOGICAL,INTENT(in) :: stress_ud_wake_on LOGICAL,INTENT(in) :: stress_vd_wake_on LOGICAL,INTENT(in) :: du_dt_satn_on LOGICAL,INTENT(in) :: du_dt_satn_p_on LOGICAL,INTENT(in) :: dv_dt_satn_on LOGICAL,INTENT(in) :: du_dt_wake_on LOGICAL,INTENT(in) :: dv_dt_wake_on !-------------------------------------------------------------------- ! LOCAL ARRAYS AND SCALARS !-------------------------------------------------------------------- INTEGER i,k ! loop counter in routine ! REAL(KIND=r8) & ! & unit_x(points) & ! ! x_compnt of unit stress vector ! &,unit_y(points) ! y_compnt of unit stress vector ! function and subroutine calls: ! EXTERNAL gw_satn,gw_wake !------------------------------------------------------------------- ! 1.0 start preliminaries ! initialise increment and increment diagnostics !------------------------------------------------------------ DO k=1,kMax DO i=1,points du_dt(i,k)=0.0_r8 dv_dt(i,k)=0.0_r8 END DO IF( du_dt_satn_on .OR. du_dt_satn_p_on ) THEN DO i=1,points du_dt_satn(i,k)=0.0_r8 END DO END IF IF( dv_dt_satn_on ) THEN DO i=1,points dv_dt_satn(i,k)=0.0_r8 END DO END IF IF( du_dt_wake_on ) THEN DO i=1,points du_dt_wake(i,k)=0.0_r8 END DO END IF IF( dv_dt_wake_on ) THEN DO i=1,points dv_dt_wake(i,k)=0.0_r8 END DO END IF END DO ! kMax ! ! and now for stress diagnostics ! DO k=0,kMax IF (stress_ud_on .OR. stress_ud_p_on ) THEN DO i=1,points stress_ud(i,k) = 0.0_r8 END DO END IF IF (stress_vd_on ) THEN DO i=1,points stress_vd(i,k) = 0.0_r8 END DO END IF IF (stress_ud_satn_on ) THEN DO i=1,points stress_ud_satn(i,k) = 0.0_r8 END DO END IF IF (stress_vd_satn_on ) THEN DO i=1,points stress_vd_satn(i,k) = 0.0_r8 END DO END IF IF (stress_ud_wake_on ) THEN DO i=1,points stress_ud_wake(i,k) = 0.0_r8 END DO END IF IF (stress_vd_wake_on ) THEN DO i=1,points stress_vd_wake(i,k) = 0.0_r8 END DO END IF END DO !--------------------------------------------------------------------- ! 2. launch linear hydrostatic waves from level k_top and ! calculate stress profile due to wave saturation effects. !--------------------------------------------------------------------- ! DEPENDS ON: gw_satn CALL gw_satn(& rho (1:points,1:kMax) , & !REAL,intent(in):: rho(points,kMax)! density r_rho_levels (1:points,1:kMax) , & !REAL,intent(in):: r_rho_levels(points,kMax)! heights above z=0 on rho levels r_theta_levels (1:points,1:kMax) , & !REAL,intent(in):: r_theta_levels(points,kMax)! heights above z=0 on theta levels theta (1:points,1:kMax) , & !REAL,intent(in):: theta(points,kMax) ! theta field u (1:points,1:kMax) , & !REAL,intent(in):: u(points,kMax) ! u field v (1:points,1:kMax) , & !REAL,intent(in):: v(points,kMax) ! v field s_x_lin_stress (1:points) , & !REAL,intent(in):: s_x_stress(points)! linear surface x_stress s_y_lin_stress (1:points) , & !REAL,intent(in):: s_y_stress(points) ! linear surface y_stress kMax , & !INTEGER,intent(in):: kMax ! number of model levels points , & !INTEGER,intent(in):: points ! number of points kay , & !REAL,intent(in):: kay ! GWD Constant (m-1) sd_orog (1:points) , & !REAL,intent(in):: sd_orog(points) ! standard deviation of orography s_x_orog (1:points) , & !REAL,intent(in):: s_x_orog(points)! 'surface' x_orog s_y_orog (1:points) , & !REAL,intent(in):: s_y_orog(points)! 'surface' y_orog du_dt (1:points,1:kMax) , & !REAL,intent(out):: du_dt(points,kMax)! total GWD du/dt dv_dt (1:points,1:kMax) , & !REAL,intent(out):: dv_dt(points,kMax) ! total GWD dv/dt k_top (1:points) , & !INTEGER,intent(in):: k_top(points) ! model level at mountain tops k_top_max , & !INTEGER,intent(in):: k_top_max ! max(k_top) fr (1:points) , & !REAL,intent(in):: fr(points) !in low level Froude number l_drag (1:points) , & !LOGICAL,intent(in):: l_drag(points) ! whether a point has a rho_s (1:points) , & !REAL,intent(in):: rho_s(points) ! surface density (calculated in gwsurf) l_fix_gwsatn , & !LOGICAL,intent(in):: l_fix_gwsatn ! Switch at vn6.2 - if true then l_gwd_40km , & !LOGICAL,intent(in):: l_gwd_40km ! Switch at vn6.6 - if true then ! diagnostics stress_ud (1:points_stress_ud,0:kMax) , & !REAL,intent(out):: stress_ud(points_stress_ud,0:kMax)! u-stress diagnostic points_stress_ud , & !INTEGER,intent(in):: points_stress_ud stress_ud_on , & !LOGICAL,intent(in):: stress_ud_on stress_ud_p_on , & !LOGICAL,intent(in):: stress_ud_p_on stress_vd (1:points_stress_vd,0:kMax) , & !REAL,intent(out):: stress_vd(points_stress_vd,0:kMax)! v-stress diagnostic points_stress_vd , & !INTEGER,intent(in):: points_stress_vd stress_vd_on , & !LOGICAL,intent(in):: stress_vd_on stress_ud_satn (1:points_stress_ud_satn,0:kMax), & !REAL,intent(out):: stress_ud_satn(points_stress_ud_satn,0:kMax)! x saturation stress diag points_stress_ud_satn , & !INTEGER,intent(in):: points_stress_ud_satn stress_ud_satn_on , & !LOGICAL,intent(in):: stress_ud_satn_on stress_vd_satn (1:points_stress_vd_satn,0:kMax), & !REAL,intent(out):: stress_vd_satn(points_stress_vd_satn,0:kMax)! y saturation stress diag points_stress_vd_satn , & !INTEGER,intent(in):: points_stress_vd_satn stress_vd_satn_on , & !LOGICAL,intent(in):: stress_vd_satn_on du_dt_satn (1:points_du_dt_satn,1:kMax) , & !REAL,intent(out):: du_dt_satn(points_du_dt_satn,kMax)! Saturation du/dt points_du_dt_satn , & !INTEGER,intent(in):: points_du_dt_satn du_dt_satn_on , & !LOGICAL,intent(in):: du_dt_satn_on du_dt_satn_p_on , & !LOGICAL,intent(in):: du_dt_satn_p_on dv_dt_satn (1:points_dv_dt_satn,1:kMax) , & !REAL,intent(out):: dv_dt_satn(points_dv_dt_satn,kMax)! Saturation dv/dt points_dv_dt_satn , & !INTEGER,intent(in):: points_dv_dt_satn dv_dt_satn_on ) !LOGICAL,intent(in):: dv_dt_satn_on ! !------------------------------------------------------------------ ! 3 apply uniform stress reduction between surface and k_top ! for points where blocked flow was diagnosed. !------------------------------------------------------------------ ! DEPENDS ON: gw_wake CALL gw_wake( & s_x_wake_stress (1:points), &!REAL,intent(in):: s_x_stress(points) ! surface x_stress s_y_wake_stress (1:points), &!REAL,intent(in):: s_y_stress(points) ! surface y_stress kMax , &!INTEGER,intent(in):: kMax ! number of model levels rho (1:points,1:kMax), &!REAL,intent(in):: rho(points,kMax) ! density on rho levels ! r_rho_levels , &!REAL,intent(in):: r_rho_levels(points,kMax)! heights above z=0 on rho levels r_theta_levels (1:points,1:kMax), &!REAL,intent(in):: r_theta_levels(points,kMax)! heights above z=0 on theta levels points , &!INTEGER,intent(in):: points ! number of GWD points k_top (1:points), &!INTEGER,intent(in):: k_top(points) ! level of blocked layer top l_drag (1:points), &!LOGICAL,intent(in ) :: l_drag(points) ! true if a non-zero surface stress du_dt (1:points,1:kMax), &!REAL,intent(inout):: du_dt(points,kMax)! total GWD du/dt dv_dt (1:points,1:kMax), &!REAL,intent(inout):: dv_dt(points,kMax)! total GWD dv/dt ! diagnostics stress_ud (1:points_stress_ud,0:kMax), &!REAL ,INTENT(OUT ) :: stress_ud(points_stress_ud,0:kMax)! u stress diagnostic points_stress_ud , &!INTEGER,intent(in ) :: points_stress_ud stress_ud_on , &!LOGICAL,intent(in ) :: stress_ud_on stress_ud_p_on , &!LOGICAL,intent(in ) :: stress_ud_p_on stress_vd (1:points_stress_vd,0:kMax), &!REAL ,INTENT(OUT ) :: stress_vd(points_stress_vd,0:kMax)! v stress diagnostic points_stress_vd , &!INTEGER,intent(in ) :: points_stress_vd stress_vd_on , &!LOGICAL,intent(in ) :: stress_vd_on stress_ud_wake (1:points_stress_ud_wake,0:kMax), &!REAL ,INTENT(OUT ) :: stress_ud_wake(points_stress_ud_wake,0:kMax)! u blocked flow stress diag points_stress_ud_wake , &!INTEGER,intent(in ) :: points_stress_ud_wake stress_ud_wake_on , &!LOGICAL,intent(in ) :: stress_ud_wake_on stress_vd_wake (1:points_stress_vd_wake,0:kMax), &!REAL ,INTENT(OUT ) :: stress_vd_wake(points_stress_vd_wake,0:kMax)! v blocked flow stress diag points_stress_vd_wake , &!INTEGER,intent(in ) :: points_stress_vd_wake stress_vd_wake_on , &!LOGICAL,intent(in ) :: stress_vd_wake_on du_dt_wake (1:points_du_dt_wake,1:kMax), &!REAL ,INTENT(OUT ) :: du_dt_wake(points_du_dt_wake,kMax)! u-acceln diagnostic points_du_dt_wake , &!INTEGER,intent(in ) :: points_du_dt_wake du_dt_wake_on , &!LOGICAL,intent(in ) :: du_dt_wake_on dv_dt_wake (1:points_dv_dt_wake,1:kMax), &!REAL ,INTENT(OUT ) :: dv_dt_wake(points_dv_dt_wake,kMax)! v-acceln diagnostic points_dv_dt_wake , &!INTEGER,intent(in ) :: points_dv_dt_wake dv_dt_wake_on )!LOGICAL,intent(in ) :: dv_dt_wake_on RETURN END SUBROUTINE gw_vert ! *****************************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 gw_wake to deposit blocked flow stress uniformly ! from ground up to sub-grid mountain tops ! SUBROUTINE gw_wake( & s_x_stress , &!REAL,intent(in):: s_x_stress(points) ! surface x_stress s_y_stress , &!REAL,intent(in):: s_y_stress(points) ! surface y_stress kMax , &!INTEGER,intent(in):: kMax ! number of model levels rho , &!REAL,intent(in):: rho(points,kMax) ! density on rho levels ! r_rho_levels , &!REAL,intent(in):: r_rho_levels(points,kMax)! heights above z=0 on rho levels r_theta_levels , &!REAL,intent(in):: r_theta_levels(points,kMax)! heights above z=0 on theta levels points , &!INTEGER,intent(in):: points ! number of GWD points k_top , &!INTEGER,intent(in):: k_top(points) ! level of blocked layer top l_drag , &!LOGICAL,intent(in ) :: l_drag(points) ! true if a non-zero surface stress du_dt , &!REAL,intent(inout):: du_dt(points,kMax)! total GWD du/dt dv_dt , &!REAL,intent(inout):: dv_dt(points,kMax)! total GWD dv/dt ! diagnostics stress_ud , &!REAL ,INTENT(OUT ) :: stress_ud(points_stress_ud,0:kMax)! u stress diagnostic points_stress_ud , &!INTEGER,intent(in ) :: points_stress_ud stress_ud_on , &!LOGICAL,intent(in ) :: stress_ud_on stress_ud_p_on , &!LOGICAL,intent(in ) :: stress_ud_p_on stress_vd , &!REAL ,INTENT(OUT ) :: stress_vd(points_stress_vd,0:kMax)! v stress diagnostic points_stress_vd , &!INTEGER,intent(in ) :: points_stress_vd stress_vd_on , &!LOGICAL,intent(in ) :: stress_vd_on stress_ud_wake , &!REAL ,INTENT(OUT ) :: stress_ud_wake(points_stress_ud_wake,0:kMax)! u blocked flow stress diag points_stress_ud_wake, &!INTEGER,intent(in ) :: points_stress_ud_wake stress_ud_wake_on , &!LOGICAL,intent(in ) :: stress_ud_wake_on stress_vd_wake , &!REAL ,INTENT(OUT ) :: stress_vd_wake(points_stress_vd_wake,0:kMax)! v blocked flow stress diag points_stress_vd_wake, &!INTEGER,intent(in ) :: points_stress_vd_wake stress_vd_wake_on , &!LOGICAL,intent(in ) :: stress_vd_wake_on du_dt_wake , &!REAL ,INTENT(OUT ) :: du_dt_wake(points_du_dt_wake,kMax)! u-acceln diagnostic points_du_dt_wake , &!INTEGER,intent(in ) :: points_du_dt_wake du_dt_wake_on , &!LOGICAL,intent(in ) :: du_dt_wake_on dv_dt_wake , &!REAL ,INTENT(OUT ) :: dv_dt_wake(points_dv_dt_wake,kMax)! v-acceln diagnostic points_dv_dt_wake , &!INTEGER,intent(in ) :: points_dv_dt_wake dv_dt_wake_on )!LOGICAL,intent(in ) :: dv_dt_wake_on IMPLICIT NONE ! Description: deposit blocked flow stress uniformly between ! the surface and the blocked layer top (k_top) ! ! current code owner: S. Webster ! ! history: ! version date comment ! 5.2 15/11/00 Original code. Stuart Webster ! 5.3 16/10/01 Change k_bot to k_top. Stuart Webster ! 5.4 28/08/02 Correct header comments. Stuart Webster ! ! ! code description: ! language: fortran 77 + common extensions ! this code is written to umdp3 v6 programming standards. ! suitable for single column use,rotated grids ! local constants ! none ! ! SUBROUTINE ARGUMENTS ! INTEGER,INTENT(in):: kMax ! number of model levels INTEGER,INTENT(in):: points ! number of GWD points ! ! integers below set diagnostic array sizes to points if diagnostic ! is called or to 1 if not. These are set in GWD_CTL2 ! INTEGER,INTENT(in):: points_stress_ud INTEGER,INTENT(in):: points_stress_vd INTEGER,INTENT(in):: points_stress_ud_wake INTEGER,INTENT(in):: points_stress_vd_wake INTEGER,INTENT(in):: points_du_dt_wake INTEGER,INTENT(in):: points_dv_dt_wake INTEGER,INTENT(in):: k_top(points) ! level of blocked layer top REAL(KIND=r8),INTENT(in):: s_x_stress(points) ! surface x_stress REAL(KIND=r8),INTENT(in):: s_y_stress(points) ! surface y_stress ! REAL(KIND=r8),intent(in):: r_rho_levels(points,kMax)! heights above z=0 on rho levels REAL(KIND=r8),INTENT(in):: r_theta_levels(points,kMax)! heights above z=0 on theta levels REAL(KIND=r8),INTENT(in):: rho(points,kMax) ! density on rho levels REAL(KIND=r8),INTENT(inout):: du_dt(points,kMax)! total GWD du/dt REAL(KIND=r8),INTENT(inout):: dv_dt(points,kMax)! total GWD dv/dt ! ! Diagnostic arrays ! REAL(KIND=r8), INTENT(OUT ) :: du_dt_wake(points_du_dt_wake,kMax)! u-acceln diagnostic REAL(KIND=r8), INTENT(OUT ) :: dv_dt_wake(points_dv_dt_wake,kMax)! v-acceln diagnostic REAL(KIND=r8), INTENT(INOUT) :: stress_ud(points_stress_ud,0:kMax)! u stress diagnostic REAL(KIND=r8), INTENT(INOUT) :: stress_vd(points_stress_vd,0:kMax)! v stress diagnostic REAL(KIND=r8), INTENT(OUT ) :: stress_ud_wake(points_stress_ud_wake,0:kMax)! u blocked flow stress diag REAL(KIND=r8), INTENT(OUT ) :: stress_vd_wake(points_stress_vd_wake,0:kMax)! v blocked flow stress diag LOGICAL,INTENT(in ) :: l_drag(points) ! true if a non-zero surface stress ! ! switches for diagnostics. ! LOGICAL,INTENT(in) :: stress_ud_on LOGICAL,INTENT(in) :: stress_ud_p_on LOGICAL,INTENT(in) :: stress_vd_on LOGICAL,INTENT(in) :: stress_ud_wake_on LOGICAL,INTENT(in) :: stress_vd_wake_on LOGICAL,INTENT(in) :: du_dt_wake_on LOGICAL,INTENT(in) :: dv_dt_wake_on !---------------------------------------------------------------- ! LOCAL ARRAYS AND SCALARS !---------------------------------------------------------------- REAL(KIND=r8) :: x_stress (points,2)! x_stresses (layer boundaries) REAL(KIND=r8) :: y_stress (points,2)! y_stresses (layer boundaries) REAL(KIND=r8) :: dz_x_stress(points)! x component of stress gradient REAL(KIND=r8) :: dz_y_stress(points)! y component of stress gradient REAL(KIND=r8) :: delta_z ! difference in height across layer(s) INTEGER :: i,k ! loop counters INTEGER :: kk,kl,ku ! level counters in routine x_stress =0.0_r8 y_stress =0.0_r8 dz_x_stress=0.0_r8 dz_y_stress=0.0_r8 delta_z =0.0_r8 delta_z =0.0_r8 ! function and subroutine calls ! none !------------------------------------------------------------------- ! 1. start level preliminaries !------------------------------------------------------------------- kl=1 ku=2 DO i=1,points IF( l_drag(i) ) THEN delta_z = r_theta_levels(i,k_top(i)) dz_x_stress(i) = s_x_stress(i) / delta_z dz_y_stress(i) = s_y_stress(i) / delta_z END IF ! if l_drag END DO ! points IF( stress_ud_on .OR. stress_ud_p_on .OR. stress_ud_wake_on) THEN DO i=1,points IF( l_drag(i) ) THEN x_stress(i,kl) = s_x_stress(i) END IF END DO IF( stress_ud_on .OR. stress_ud_p_on) THEN DO i=1,points IF( l_drag(i) ) THEN stress_ud(i,0) = stress_ud(i,0) + x_stress(i,kl) END IF END DO END IF IF( stress_ud_wake_on ) THEN DO i=1,points IF( l_drag(i) ) THEN stress_ud_wake(i,0) = x_stress(i,kl) END IF END DO END IF END IF ! stress_ud_on .or. stress_ud_p_on .or. stress_ud_wake_on IF( stress_vd_on .OR. stress_vd_wake_on) THEN DO i=1,points IF( l_drag(i) ) THEN y_stress(i,kl) = s_y_stress(i) END IF END DO IF( stress_vd_on ) THEN DO i=1,points IF( l_drag(i) ) THEN stress_vd(i,0) = stress_vd(i,0) + y_stress(i,kl) END IF END DO END IF IF( stress_vd_wake_on ) THEN DO i=1,points IF( l_drag(i) ) THEN stress_vd_wake(i,0) = y_stress(i,kl) END IF END DO END IF END IF ! stress_vd_on .or. stress_vd_wake_on !------------------------------------------------------------------ ! 2 loop levels ! k is the rho level counter !------------------------------------------------------------------ DO k=1,kMax DO i=1,points IF( l_drag(i) .AND. k <= k_top(i) ) THEN ! ! note that a constant stress drop across a level in height coords => ! a non-uniform drag because of rho variation. ! du_dt(i,k) = du_dt(i,k) - dz_x_stress(i)/rho(i,k) dv_dt(i,k) = dv_dt(i,k) - dz_y_stress(i)/rho(i,k) !PRINT*,k, du_dt(i,k),dv_dt(i,k),rho(i,k) END IF ! if l_drag(i) .and. k<=k_top(i) END DO ! points ! diagnostics IF( du_dt_wake_on ) THEN DO i=1,points IF( l_drag(i) .AND. k <= k_top(i) ) THEN du_dt_wake(i,k) = - dz_x_stress(i)/rho(i,k) END IF END DO END IF IF( dv_dt_wake_on ) THEN DO i=1,points IF( l_drag(i) .AND. k <= k_top(i) ) THEN dv_dt_wake(i,k) = - dz_y_stress(i)/rho(i,k) END IF END DO END IF IF( stress_ud_on .OR. stress_ud_p_on .OR. stress_ud_wake_on) THEN DO i=1,points IF( l_drag(i) .AND. k <= k_top(i) ) THEN IF ( k == 1 ) THEN delta_z = r_theta_levels(i,k) ELSE delta_z= r_theta_levels(i,k) - r_theta_levels(i,k-1) END IF x_stress(i,ku) = x_stress(i,kl) x_stress(i,ku) = x_stress(i,kl)-dz_x_stress(i)*delta_z END IF ! l_drag & k < k_top END DO IF( stress_ud_on .OR. stress_ud_p_on) THEN DO i=1,points IF( l_drag(i) .AND. k <= k_top(i) ) THEN stress_ud(i,k) = stress_ud(i,k) + x_stress(i,ku) END IF ! l_drag & k < k_top END DO END IF ! stress_ud on IF( stress_ud_wake_on ) THEN DO i=1,points IF( l_drag(i) .AND. k <= k_top(i) ) THEN stress_ud_wake(i,k) = x_stress(i,ku) END IF ! l_drag & k < k_top END DO END IF ! stress_ud_wake on END IF ! stress_ud_on .or. stress_ud_p_on .or. stress_wake_ud_on IF( stress_vd_on .OR. stress_vd_wake_on) THEN DO i=1,points IF( l_drag(i) .AND. k <= k_top(i) ) THEN IF ( k == 1 ) THEN delta_z = r_theta_levels(i,k) ELSE delta_z= r_theta_levels(i,k) - r_theta_levels(i,k-1) END IF y_stress(i,ku) = y_stress(i,kl) y_stress(i,ku) = y_stress(i,kl)-dz_y_stress(i)*delta_z END IF ! l_drag & k < k_top END DO IF( stress_vd_on ) THEN DO i=1,points IF( l_drag(i) .AND. k <= k_top(i) ) THEN stress_vd(i,k) = stress_vd(i,k) + y_stress(i,ku) END IF ! l_drag & k < k_top END DO END IF ! stress_vd on IF( stress_vd_wake_on ) THEN DO i=1,points IF( l_drag(i) .AND. k <= k_top(i) ) THEN stress_vd_wake(i,k) = y_stress(i,ku) END IF ! l_drag & k < k_top END DO END IF ! stress_vd_wake on END IF ! stress_vd_on .or. stress_wake_vd_on ! swap storage for lower and upper layers kk=kl kl=ku ku=kk END DO ! end loop kMax RETURN END SUBROUTINE gw_wake ! *****************************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 gw_surf to calculate surface stress vector for gwd. ! also, calculate blocked layer surface stress. ! SUBROUTINE gw_surf(& r_theta_levels , &!REAL ,intent(in):: r_theta_levels(points,kMax)! heights on theta levels rho , &!REAL ,intent(in):: rho(points,kMax) ! density on rho levels theta , &!REAL ,intent(in):: theta(points,kMax) ! theta on theta levels u , &!REAL ,intent(in):: u(points,kMax) ! u on rho levels v , &!REAL ,intent(in):: v(points,kMax) ! v on rho levels timestep , &!REAL ,intent(in):: timestep ! timestep (s) sd_orog , &!REAL,intent(in):: sd_orog(points) ! standard deviation of the sub-grid orography sigma_xx , &!REAL,intent(in):: sigma_xx(points) ! (dh/dx)^2 grid box average of the sigma_xy , &!REAL,intent(in):: sigma_xy(points) ! (dh/dx)*(dh/dy) sigma_yy , &!REAL,intent(in):: sigma_yy(points) ! (dh/dy)^2 s_x_lin_stress , &!REAL,intent(out):: s_x_lin_stress(points)! 'surface' linear stress in x-dirn s_y_lin_stress , &!REAL,intent(out):: s_y_lin_stress(points)! 'surface' linear stress in y-dirn s_x_wake_stress , &!REAL,intent(out):: s_x_wake_stress(points)! wake stress in x-dirn s_y_wake_stress , &!REAL,intent(out):: s_y_wake_stress(points)! wake stress in y-dirn s_x_orog , &!REAL,intent(out):: s_x_orog(points) ! 'surface' stress/orog in x-dirn s_y_orog , &!REAL ,intent(out):: s_y_orog(points) ! 'surface' stress/orog in y-dirn kMax , &!INTEGER ,intent(in) :: kMax ! number of model levels points , &!INTEGER ,intent(in) :: points ! number of gwd points on pe kay , &!REAL ,intent(in):: kay ! gwd surface stress constant (m-1) rho_s , &!REAL,intent(out):: rho_s(points) ! density - av from z=0 to nsigma*sd_orog l_taus_scale , &!LOGICAL,intent(in):: l_taus_scale ! true allows variation of surface stress on the low level Froude number k_top , &!INTEGER ,intent(out):: k_top(points) ! topmost model level including mountains. k_top_max , &!INTEGER ,intent(out):: k_top_max ! max(k_top) lift , &! REAL,intent(out):: lift(points) ! blocked layer depth l_drag , &!LOGICAL,intent(out):: l_drag(points) ! whether a point has a non-zero surface stress or not fr , &!REAL,intent(out):: fr(points) ! low level froude number for linear frc , &!REAL ,intent(in):: frc ! critical froude number below which hydraulic jumps are triggered u_s_d , &!REAL,intent(out):: u_s_d(points_u_s_d) ! 0-nsigma*sd_orog u_s diag u_s_d_on , &!LOGICAL,intent(in):: u_s_d_on points_u_s_d , &!INTEGER ,intent(in) ::points_u_s_d v_s_d , &!REAL,intent(out):: v_s_d(points_v_s_d) ! 0-nsigma_sd_orog v_s diag v_s_d_on , &!LOGICAL,intent(in):: v_s_d_on points_v_s_d , &!INTEGER ,intent(in) ::points_v_s_d nsq_s_d , &!REAL,intent(out):: nsq_s_d(points_nsq_s_d)! 0-nsigma*sd_orog nsq_s diagnostic nsq_s_d_on , &!LOGICAL,intent(in):: nsq_s_d_on points_nsq_s_d , &!INTEGER ,intent(in) ::points_nsq_s_d fr_d , &!REAL,intent(out):: fr_d(points_fr_d) ! Froude no. diagnostic fr_d_on , &!LOGICAL,intent(in):: fr_d_on points_fr_d , &!INTEGER ,intent(in) ::points_fr_d bld_d , &!REAL,intent(out):: bld_d(points_bld_d) ! blocked layer depth diagnostic bld_d_on , &!LOGICAL,intent(in):: bld_d_on points_bld_d , &!INTEGER ,intent(in) ::points_bld_d bldt_d , &!REAL,intent(out):: bldt_d(points_bldt_d) ! % of time blocked flow param. invoked bldt_d_on , &!LOGICAL,intent(in):: bldt_d_on points_bldt_d , &!INTEGER ,intent(in) ::points_bldt_d num_lim_d , &!REAL,intent(out):: num_lim_d(points_num_lim_d)! % of time numerical limiter invoked num_lim_d_on , &!LOGICAL,intent(in):: num_lim_d_on points_num_lim_d , &!INTEGER ,intent(in) ::points_num_lim_d num_fac_d , &!REAL,intent(out):: num_fac_d(points_num_fac_d)! % reduction of flow blocking stressafter numerical limiter invoked num_fac_d_on , &!LOGICAL,intent(in):: num_fac_d_on points_num_fac_d , &!INTEGER ,intent(in) ::points_num_fac_d tausx_d , &!REAL,intent(out):: tausx_d(points_tausx_d) ! x-component of total surface stress tausx_d_on , &!LOGICAL,intent(in):: tausx_d_on points_tausx_d , &!INTEGER ,intent(in) ::points_tausx_d tausy_d , &!REAL,intent(out):: tausy_d(points_tausy_d) ! y-component of total surface stress tausy_d_on , &!LOGICAL,intent(in):: tausy_d_on points_tausy_d , &!INTEGER ,intent(in) ::points_tausy_d taus_scale_d , &!REAL,intent(out):: taus_scale_d(points_taus_scale_d)! scaling factor for surface stress when Fr dependence of surface stress invoked. taus_scale_d_on , &!LOGICAL,intent(in):: taus_scale_d_on points_taus_scale_d)!INTEGER ,intent(in) ::points_taus_scale_d IMPLICIT NONE ! Description: ! calculates the total surface stress using the ! linear hydrostatic GWD surface stress equation ! for a 2-d hill in the absence of rotation and ! friction. The blocked layer depth is then used ! to partition this stress into a linear hydrostatic ! GWD component and a blocked flow component. ! nsigma*sd_orog is taken to be the top of the sub-grid ! mountains and the low level U and N are calculated ! as averages over this layer. This U and N are used ! in the calculation of the surface stress and in the ! determination of the blocked layer and low level ! Froude number. ! From vn6.2 it is possible to allow the surface drag to ! be dependent on the low level Froude Number. This is ! motivated by the results of ! Wells et al.2005 (QJRMS, 131, 1321-1338). ! ! Current code owner: S. Webster ! ! history: ! version date comment ! 5.2 15/11/00 original code. loosely based on gwsurf3a at vn5.1. ! major overhaul though. see above for more details. ! Stuart Webster ! ! 5.3 16/10/01 Another major overhaul. Total stress now calculated ! using the usual linear hydrostatic GWD surface ! stress equation and this is then split into a ! blocked flow stress and a linear hydrostatic ! GWD stress. ! Stuart Webster ! 5.4 28/08/02 Introduce numerical limiter and associated ! diagnostics for flow blocking scheme. ! Stuart Webster ! 6.2 15/08/05 Free format fixes. P.Selwood ! 6.2 21/02/06 Introduce surface stress diagnostics. S. Webster ! 6.2 21/02/06 Introduce dependence of the total surface drag on ! the low level Froude number. Stuart Webster ! ! code description: ! language: fortran 77 + common extensions ! this code is written to umdp3 v6 programming standards. ! suitable for single column use,rotated grids ! global variables !*L------------------COMDECK C_G---------------------------------------- ! G IS MEAN ACCEL DUE TO GRAVITY AT EARTH'S SURFACE ! REAL(KIND=r8), Parameter :: G = 9.80665 !*---------------------------------------------------------------------- !*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 !*---------------------------------------------------------------------- ! local constants ! ! Description: This comdeck defines the constants for the 3B and 4A ! versions of the Gravity Wave Drag Code. These are ! tuneable parameters but are unlikely to be changed. ! ! History: ! Version Date Comment ! ------- ---- ------- ! 3.4 18/10/94 Original Version J.R. Mitchell ! 4.3 7/03/97 Remove KAY_LEE (now set in RUNCNST) S.Webster ! 4.5 03/08/98 Add GAMMA_SATN (Used in 06_3B). D. Robinson ! 5.0 14/07/99 Remove redundant switches/variables: only keep ! s-version 06_3B parameters. R Rawlins ! 5.2 15/11/00 Set parameters for the 4A scheme. ! Stuart Webster. ! 5.3 16/10/01 Partition 3B and 4A scheme parameters. ! Stuart Webster. ! 5.3 21/09/01 Set parameters for the spectral (middle ! atmosphere) gravity wave forcing scheme. ! Warner and McIntyre, J. Atm. Sci., 2001, ! Scaife et al., Geophys. Res. Lett., 2000 ! Scaife et al, J. Atm. Sci., 2002 give ! further details. ! Adam Scaife ! 5.4 Alters c_gwave.h include file to change from launch level ! to launch eta to make less model level dependent. ! Adam Scaife ! 5.5 25/02/03 Remove 3B GWD parameter settings. Stuart Webster ! ! 6.2 16/06/05 Move CCL parameter to gw_ussp. Andrew Bushell ! ! ! Number of standard deviations above the mean orography of top ! of sub-grid mountains !REAL(KIND=r8),PARAMETER :: NSIGMA = 2.5_r8 ! SPECTRAL GRAVITY WAVE SCHEME PARAMETERS: ! LMINL = 1/max wavelength at launch ! LSTAR = 1/characteristic wavelength ! ETALAUNCH = eta for model launch REAL(KIND=r8),PARAMETER:: LMINL = 1.0_r8/20000.0_r8 REAL(KIND=r8),PARAMETER:: LSTAR = 1.0_r8/4300.0_r8 REAL(KIND=r8),PARAMETER:: ETALAUNCH = 0.045_r8 ! ! SUBROUTINE ARGUMENTS ! INTEGER ,INTENT(in) :: kMax ! number of model levels INTEGER ,INTENT(in) :: points ! number of gwd points on pe INTEGER ,INTENT(out):: k_top(points) ! topmost model level including mountains. ! ! defined to be the kth rho level in ! ! which the nsigma*sd_orog height lies, ! ! i.e.nsigma*sd_orog lies between ! ! r_theta_levels(k-1) and ! ! r_theta_levels(k). INTEGER ,INTENT(out):: k_top_max ! max(k_top) ! ! Dimensions of diagnostic arrays. ! dimension = points if diag called, ! = 1 if diag not called. ! These are set in GWD_CTL2 ! INTEGER ,INTENT(in) ::points_nsq_s_d INTEGER ,INTENT(in) ::points_u_s_d INTEGER ,INTENT(in) ::points_v_s_d INTEGER ,INTENT(in) ::points_fr_d INTEGER ,INTENT(in) ::points_bld_d INTEGER ,INTENT(in) ::points_bldt_d INTEGER ,INTENT(in) ::points_num_lim_d INTEGER ,INTENT(in) ::points_num_fac_d INTEGER ,INTENT(in) ::points_tausx_d INTEGER ,INTENT(in) ::points_tausy_d INTEGER ,INTENT(in) ::points_taus_scale_d REAL(KIND=r8) ,INTENT(in):: r_theta_levels(points,kMax)! heights on theta levels REAL(KIND=r8) ,INTENT(in):: rho(points,kMax) ! density on rho levels REAL(KIND=r8) ,INTENT(in):: theta(points,kMax) ! theta on theta levels REAL(KIND=r8) ,INTENT(in):: u(points,kMax) ! u on rho levels REAL(KIND=r8) ,INTENT(in):: v(points,kMax) ! v on rho levels REAL(KIND=r8) ,INTENT(in):: timestep ! timestep (s) REAL(KIND=r8) ,INTENT(in):: kay ! gwd surface stress constant (m-1) REAL(KIND=r8) ,INTENT(in):: frc ! critical froude number below which ! ! hydraulic jumps are triggered REAL(KIND=r8),INTENT(in):: sd_orog(points) ! standard deviation of the sub-grid ! ! orography REAL(KIND=r8),INTENT(in):: sigma_xx(points) ! (dh/dx)^2 grid box average of the ! ! the hi-res source dataset. here, ! ! dh is the height change of the hi-res ! ! data and dx is the distance between ! ! adjacent data points. REAL(KIND=r8),INTENT(in):: sigma_xy(points) ! (dh/dx)*(dh/dy) REAL(KIND=r8),INTENT(in):: sigma_yy(points) ! (dh/dy)^2 REAL(KIND=r8),INTENT(out):: s_x_lin_stress (points) ! 'surface' linear stress in x-dirn REAL(KIND=r8),INTENT(out):: s_y_lin_stress (points) ! 'surface' linear stress in y-dirn REAL(KIND=r8),INTENT(out):: s_x_wake_stress (points)! wake stress in x-dirn REAL(KIND=r8),INTENT(out):: s_y_wake_stress (points)! wake stress in y-dirn REAL(KIND=r8),INTENT(out):: s_x_orog (points) ! 'surface' stress/orog in x-dirn REAL(KIND=r8),INTENT(out):: s_y_orog (points) ! 'surface' stress/orog in y-dirn REAL(KIND=r8),INTENT(out):: rho_s (points) ! density - av from z=0 to nsigma*sd_orog REAL(KIND=r8),INTENT(out):: lift (points) ! blocked layer depth REAL(KIND=r8),INTENT(out):: fr (points) ! low level froude number for linear ! ! hydrostatic waves - so nsigma*sd_orog ! ! is replaced by nsigma*sd_orog-lift REAL(KIND=r8),INTENT(out):: u_s_d (points_u_s_d) ! 0-nsigma*sd_orog u_s diag REAL(KIND=r8),INTENT(out):: v_s_d (points_v_s_d) ! 0-nsigma_sd_orog v_s diag REAL(KIND=r8),INTENT(out):: nsq_s_d (points_nsq_s_d)! 0-nsigma*sd_orog nsq_s diagnostic REAL(KIND=r8),INTENT(out):: fr_d ( points_fr_d) ! Froude no. diagnostic REAL(KIND=r8),INTENT(out):: bld_d (points_bld_d) ! blocked layer depth diagnostic REAL(KIND=r8),INTENT(out):: bldt_d (points_bldt_d) ! % of time blocked flow param. invoked REAL(KIND=r8),INTENT(out):: num_lim_d (points_num_lim_d)! % of time numerical limiter invoked REAL(KIND=r8),INTENT(out):: num_fac_d (points_num_fac_d)! % reduction of flow blocking stress ! after numerical limiter invoked REAL(KIND=r8),INTENT(out):: tausx_d (points_tausx_d) ! x-component of total surface stress REAL(KIND=r8),INTENT(out):: tausy_d (points_tausy_d) ! y-component of total surface stress REAL(KIND=r8),INTENT(out):: taus_scale_d (points_taus_scale_d)! scaling factor for surface stress when ! ! Fr dependence of surface stress invoked. LOGICAL,INTENT(in):: l_taus_scale ! true allows variation of surface stress ! ! on the low level Froude number LOGICAL,INTENT(out):: l_drag(points) ! whether a point has a non-zero surface ! ! stress or not ! ! Diagnostic switches ! LOGICAL,INTENT(in):: u_s_d_on LOGICAL,INTENT(in):: v_s_d_on LOGICAL,INTENT(in):: nsq_s_d_on LOGICAL,INTENT(in):: fr_d_on LOGICAL,INTENT(in):: bld_d_on LOGICAL,INTENT(in):: bldt_d_on LOGICAL,INTENT(in):: num_lim_d_on LOGICAL,INTENT(in):: num_fac_d_on LOGICAL,INTENT(in):: tausx_d_on LOGICAL,INTENT(in):: tausy_d_on LOGICAL,INTENT(in):: taus_scale_d_on !------------------------------------------------------------------ ! LOCAL ARRAYS AND SCALARS !------------------------------------------------------------------ REAL(KIND=r8):: u_s(points) ! u-winds - av from z=0 to nsigma*sd_orog REAL(KIND=r8):: v_s(points) ! v-winds - av from z=0 to nsigma*sd_orog REAL(KIND=r8):: nsq_s(points) ! n squared av from z=0 to nsigma*sd_orog REAL(KIND=r8):: fr_for_diag(points) ! low level froude number REAL(KIND=r8):: num_lim(points) ! 0 if limiter not invoked, 100 if it is. REAL(KIND=r8):: num_fac(points) ! percentage reduction of the ! ! flow-blocking stress after the limiter ! ! was invoked REAL(KIND=r8):: taus_scale(points) ! factor by which surf stress is scaled ! ! when Froude no. dependency is invoked. INTEGER :: i,k ! loop counter in routine REAL(KIND=r8) :: speed ! wind speed / wind speed in dirn stress REAL(KIND=r8) :: speedcalc ! numerator of calcuation for speed REAL(KIND=r8) :: s_stress_sq ! denominater of calculation for speed REAL(KIND=r8) :: s_wake_stress! amplitude of surface wake stress REAL(KIND=r8) :: s_wake_limit ! numerical limit for s_wake_stress REAL(KIND=r8) :: n ! brunt_vaisala frequency REAL(KIND=r8) :: r_frc ! reciprocal of the critical froude number REAL(KIND=r8) :: calc ! calculation for surface stress magnitude REAL(KIND=r8) :: calc1 ! as per calc REAL(KIND=r8) :: n_squared ! n^2 on rho levels REAL(KIND=r8) :: dzb ! relevant depth to vertical average ! ! of current theta level REAL(KIND=r8) :: dzt ! height from surface to kth theta level ! ! or nsigma*sd_orog LOGICAL l_cont(points) ! level continue ! function and subroutine calls: none speed = 0.0_r8 speedcalc = 0.0_r8 s_stress_sq = 0.0_r8 s_wake_stress = 0.0_r8 s_wake_limit = 0.0_r8 n = 0.0_r8 r_frc = 0.0_r8 calc = 0.0_r8 calc1 = 0.0_r8 n_squared = 0.0_r8 dzb = 0.0_r8 dzt = 0.0_r8 s_x_lin_stress = 0.0_r8 s_y_lin_stress = 0.0_r8 s_x_wake_stress = 0.0_r8 s_y_wake_stress = 0.0_r8 s_x_orog = 0.0_r8 s_y_orog = 0.0_r8 rho_s = 0.0_r8 lift = 0.0_r8 fr = 0.0_r8 u_s_d = 0.0_r8 v_s_d = 0.0_r8 nsq_s_d = 0.0_r8 fr_d = 0.0_r8 bld_d = 0.0_r8 bldt_d = 0.0_r8 num_lim_d = 0.0_r8 num_fac_d = 0.0_r8 tausx_d = 0.0_r8 tausy_d = 0.0_r8 taus_scale_d = 0.0_r8 !--------------------------------------------------------------------- ! 1.0 initialisation !--------------------------------------------------------------------- k_top_max = 1 r_frc = 1.0_r8/frc DO i=1,points nsq_s(i) = 0.0_r8 END DO DO i=1,points rho_s(i) = 0.0_r8 END DO DO i=1,points u_s(i) = 0.0_r8 END DO DO i=1,points v_s(i) = 0.0_r8 END DO DO i=1,points l_cont(i) =.TRUE. END DO DO i=1,points s_x_lin_stress(i) = 0.0_r8 END DO DO i=1,points s_y_lin_stress(i) = 0.0_r8 END DO DO i=1,points s_x_wake_stress(i) = 0.0_r8 END DO DO i=1,points s_y_wake_stress(i) = 0.0_r8 END DO DO i=1,points l_drag(i) = .TRUE. END DO DO i=1,points num_lim(i) = 0.0_r8 END DO DO i=1,points num_fac(i) = 0.0_r8 END DO DO i=1,points taus_scale(i) = 1.0_r8 END DO !------------------------------------------------------------------ ! 2. ! Calculation of the average surface quantities and ! the nsigma*sd_orog point values. The surface ! is the average from 0 up to nsigma*sd_orog. ! All interpolations now done in terms of height. ! k_top is defined to be the rho level containing nsigma*sd_orog ! dzb=dzt when k=2 because include level 1 depth in z, but ! don't calculate u and n at rho level 1. ! !--------------------------------------------------------------------- DO i=1,points IF (sd_orog(i) <= 0.0_r8 ) THEN l_drag(i) = .FALSE. l_cont(i) = .FALSE. k_top(i) = 2 END IF END DO DO k=2,kMax-1 DO i=1,points IF ( l_cont(i) ) THEN IF ( r_theta_levels(i,k) < nsigma*sd_orog(i) ) THEN dzt = r_theta_levels(i,k) IF ( k == 2 ) THEN dzb = dzt ELSE dzb = r_theta_levels(i,k) - r_theta_levels(i,k-1) END IF ELSE dzt = nsigma * sd_orog(i) IF ( k == 2 ) THEN dzb = dzt ELSE dzb = nsigma*sd_orog(i) - r_theta_levels(i,k-1) END IF k_top(i) = k IF( k_top_max < k ) k_top_max=k l_cont(i) = .FALSE. END IF ! r_theta_levels(i,k) > nsigma*sd_orog(i) n_squared = 2.0_r8*g*(theta(i,k)-theta(i,k-1) ) & & /( (theta(i,k)+theta(i,k-1)) * & & (r_theta_levels(i,k)-r_theta_levels(i,k-1)) ) ! ! Form z=0 to current level averages which ! when k=k_top become 0-nsigma*sd_orog averages ! u_s(i) = ( u_s(i) * r_theta_levels(i,k-1) + u(i,k) * dzb ) / dzt v_s(i) = ( v_s(i) * r_theta_levels(i,k-1) + v(i,k) * dzb ) / dzt nsq_s(i) = ( nsq_s(i) * r_theta_levels(i,k-1) + n_squared* dzb ) / dzt rho_s(i) = ( rho_s(i) * r_theta_levels(i,k-1) + rho(i,k) * dzb ) / dzt END IF ! l_cont END DO END DO !------------------------------------------------------------------ ! 3.1 calculation of total surface stress and blocked ! layer depth, and the split of the total stress into a linear ! hydrostatic gravity wave stress and a blocked flow stress !----------------------------------------------------------------- DO i=1,points ! ! First calculate the speed of the wind for the linear ! hydrostatic code ! speed = u_s(i)*u_s(i) + v_s(i)*v_s(i) IF ( speed <= 0.0_r8 ) THEN s_x_orog(i) = 0.0_r8 s_y_orog(i) = 0.0_r8 ELSE speed = SQRT(speed) s_x_orog(i)= (u_s(i)*sigma_xx(i) + v_s(i)*sigma_xy(i)) /speed s_y_orog(i)= (u_s(i)*sigma_xy(i) + v_s(i)*sigma_yy(i)) /speed s_stress_sq= s_x_orog(i)*s_x_orog(i) + s_y_orog(i)*s_y_orog(i) speedcalc = u_s(i)*s_x_orog(i) + v_s(i)*s_y_orog(i) IF ( s_stress_sq <= 0.0_r8 ) THEN speed = 0.0_r8 ELSE ! speed is the component of the wind perpendicular to the major ! axis of the orography. speed = speedcalc / SQRT( s_stress_sq ) END IF END IF IF ( nsq_s(i) > 0.0_r8 .AND. speed > 0.0_r8 .AND. l_drag(i) ) THEN n = SQRT( nsq_s(i) ) fr_for_diag(i) = speed / (n*nsigma*sd_orog(i)) lift(i) = nsigma * sd_orog(i) - r_frc * speed/n IF ( lift(i) < 0.0_r8 ) THEN lift(i) = 0.0_r8 END IF ! ! fr here is for the linear hydrostatic surface stress and ! then also for the critical stress calculation in gw_satn ! fr(i) = speed / ( n*(nsigma*sd_orog(i)-lift(i))) ! ! Calculate factor by which to scale surface stresses when dependence ! of surface stress on Froude number is invoked. Empirical expression ! is a simple fit to the experimental results plotted in Fig. 7 of ! Wells et al. (2005). ! ! This function was modified at vn6.6 to keep ratio of drag:linear drag ! equal to 1 for Froude numbers>2 ! IF ( l_taus_scale ) THEN IF ( fr_for_diag(i) <= 1.0_r8 ) THEN taus_scale(i) = 1.6_r8*fr_for_diag(i) ELSE IF ( fr_for_diag(i) <= 2.0_r8 ) THEN taus_scale(i) = 2.2_r8 - 0.6_r8*fr_for_diag(i) END IF END IF calc = kay * rho_s(i) * speed**3 * taus_scale(i) / & & (n*nsigma*nsigma*sd_orog(i)*sd_orog(i)*fr(i)*fr(i)) s_x_lin_stress(i) = s_x_orog(i) * calc s_y_lin_stress(i) = s_y_orog(i) * calc ! ! s_x_orog(i)*calc1 = linear hydrostatic prediction of the stress ! calc1 = kay*rho_s(i)* speed * n * taus_scale(i) s_x_wake_stress(i) = s_x_orog(i) * calc1 - s_x_lin_stress(i) s_y_wake_stress(i) = s_y_orog(i) * calc1 - s_y_lin_stress(i) ! ! Limit wake_stress so that numerical instability is not permitted ! s_wake_stress = SQRT( s_x_wake_stress(i)*s_x_wake_stress(i) & & +s_y_wake_stress(i)*s_y_wake_stress(i)) s_wake_limit = 2*nsigma*sd_orog(i)*speed*rho_s(i)/timestep IF (s_wake_stress > s_wake_limit) THEN s_x_wake_stress(i) = s_x_wake_stress(i) * s_wake_limit & & /s_wake_stress s_y_wake_stress(i) = s_y_wake_stress(i) * s_wake_limit & & /s_wake_stress num_lim(i) = 100.0_r8 num_fac(i) = 100.0_r8 * (s_wake_stress - s_wake_limit ) & & / s_wake_stress END IF ELSE l_drag(i) = .FALSE. lift(i) = 0.0_r8 fr_for_diag(i) = 0.0_r8 fr(i) = 0.0_r8 END IF ! speed or n or orog > 0.0_r8 END DO ! i=points !----------------------------------------------------------------- ! 4 diagnostics !----------------------------------------------------------------- IF ( u_s_d_on ) THEN DO i=1,points u_s_d(i) = u_s(i) END DO END IF IF ( v_s_d_on ) THEN DO i=1,points v_s_d(i) = v_s(i) END DO END IF IF ( nsq_s_d_on ) THEN DO i=1,points IF ( nsq_s(i) < 0.0_r8 ) nsq_s(i) = 0.0_r8 nsq_s_d(i) = SQRT( nsq_s(i) ) END DO END IF IF ( bld_d_on ) THEN DO i=1,points bld_d(i) = lift(i) END DO END IF IF ( bldt_d_on ) THEN DO i=1,points IF( s_x_wake_stress(i) /= 0.0_r8 .OR. & & s_y_wake_stress(i) /= 0.0_r8 ) THEN bldt_d(i) = 100.0_r8 ELSE bldt_d(i) = 0.0_r8 END IF END DO END IF IF ( fr_d_on ) THEN DO i=1,points fr_d(i) = fr_for_diag(i) ! limit Fr to sensible-ish values as previously had problems ! when creating time means which include very large values. IF (fr_d(i) > 1000.0_r8) fr_d(i) = 1000.0_r8 END DO END IF IF ( num_lim_d_on ) THEN DO i=1,points num_lim_d(i) = num_lim(i) END DO END IF IF ( num_fac_d_on ) THEN DO i=1,points num_fac_d(i) = num_fac(i) END DO END IF IF ( tausx_d_on ) THEN DO i=1,points tausx_d(i) = s_x_lin_stress(i) + s_x_wake_stress(i) END DO END IF IF ( tausy_d_on ) THEN DO i=1,points tausy_d(i) = s_y_lin_stress(i) + s_y_wake_stress(i) END DO END IF IF ( taus_scale_d_on ) THEN DO i=1,points taus_scale_d(i) = taus_scale(i) END DO END IF RETURN END SUBROUTINE gw_surf ! *****************************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 GW_SATN: Calculates stress profile for linear hydrostatic ! Waves ! SUBROUTINE gw_satn(& rho , & !REAL,intent(in):: rho(points,kMax)! density r_rho_levels , & !REAL,intent(in):: r_rho_levels(points,kMax)! heights above z=0 on rho levels r_theta_levels , & !REAL,intent(in):: r_theta_levels(points,kMax)! heights above z=0 on theta levels theta , & !REAL,intent(in):: theta(points,kMax) ! theta field u , & !REAL,intent(in):: u(points,kMax) ! u field v , & !REAL,intent(in):: v(points,kMax) ! v field s_x_stress , & !REAL,intent(in):: s_x_stress(points)! linear surface x_stress s_y_stress , & !REAL,intent(in):: s_y_stress(points) ! linear surface y_stress kMax , & !INTEGER,intent(in):: kMax ! number of model levels points , & !INTEGER,intent(in):: points ! number of points kay , & !REAL,intent(in):: kay ! GWD Constant (m-1) sd_orog , & !REAL,intent(in):: sd_orog(points) ! standard deviation of orography s_x_orog , & !REAL,intent(in):: s_x_orog(points)! 'surface' x_orog s_y_orog , & !REAL,intent(in):: s_y_orog(points)! 'surface' y_orog du_dt , & !REAL,intent(out):: du_dt(points,kMax)! total GWD du/dt dv_dt , & !REAL,intent(out):: dv_dt(points,kMax) ! total GWD dv/dt k_top , & !INTEGER,intent(in):: k_top(points) ! model level at mountain tops k_top_max , & !INTEGER,intent(in):: k_top_max ! max(k_top) fr , & !REAL,intent(in):: fr(points) !in low level Froude number l_drag , & !LOGICAL,intent(in):: l_drag(points) ! whether a point has a rho_s , & !REAL,intent(in):: rho_s(points) ! surface density (calculated in gwsurf) l_fix_gwsatn , & !LOGICAL,intent(in):: l_fix_gwsatn ! Switch at vn6.2 - if true then l_gwd_40km , & !LOGICAL,intent(in):: l_gwd_40km ! Switch at vn6.6 - if true then ! diagnostics stress_ud , & !REAL,intent(out):: stress_ud(points_stress_ud,0:kMax)! u-stress diagnostic points_stress_ud , & !INTEGER,intent(in):: points_stress_ud stress_ud_on , & !LOGICAL,intent(in):: stress_ud_on stress_ud_p_on , & !LOGICAL,intent(in):: stress_ud_p_on stress_vd , & !REAL,intent(out):: stress_vd(points_stress_vd,0:kMax)! v-stress diagnostic points_stress_vd , & !INTEGER,intent(in):: points_stress_vd stress_vd_on , & !LOGICAL,intent(in):: stress_vd_on stress_ud_satn , & !REAL,intent(out):: stress_ud_satn(points_stress_ud_satn,0:kMax)! x saturation stress diag points_stress_ud_satn , & !INTEGER,intent(in):: points_stress_ud_satn stress_ud_satn_on , & !LOGICAL,intent(in):: stress_ud_satn_on stress_vd_satn , & !REAL,intent(out):: stress_vd_satn(points_stress_vd_satn,0:kMax)! y saturation stress diag points_stress_vd_satn , & !INTEGER,intent(in):: points_stress_vd_satn stress_vd_satn_on , & !LOGICAL,intent(in):: stress_vd_satn_on du_dt_satn , & !REAL,intent(out):: du_dt_satn(points_du_dt_satn,kMax)! Saturation du/dt points_du_dt_satn , & !INTEGER,intent(in):: points_du_dt_satn du_dt_satn_on , & !LOGICAL,intent(in):: du_dt_satn_on du_dt_satn_p_on , & !LOGICAL,intent(in):: du_dt_satn_p_on dv_dt_satn , & !REAL,intent(out):: dv_dt_satn(points_dv_dt_satn,kMax)! Saturation dv/dt points_dv_dt_satn , & !INTEGER,intent(in):: points_dv_dt_satn dv_dt_satn_on ) !LOGICAL,intent(in):: dv_dt_satn_on IMPLICIT NONE ! Description: ! ! Calculate stress profiles and hence drags due to linear ! hydrostatic gravity waves. in contrast to the 3a/3b schemes, ! linear hydrostatic waves may occur for any Froude number now. ! When the Froude number drops below Frc (as specified in c_gwave.h), ! the wave amplitude is reduced to be that of the air that is not ! blocked, rather than be the full depth of the sub-grid mountains. ! Thus, as the Froude number reduces to zero, so the linear ! hydrostatic wave amplitude reduces to zero. ! ! ! current code owner: S. Webster ! ! history: ! version date comment ! 5.2 15/11/00 Original code. Based on 5.1 GWSATN3B ! code now includes the option to use the ! turning of the wind algorithm. However, ! by default, l_gwd_tow is set to false ! in c_gwave.h. The saturation part of the code ! is also now on a switch in the same header. ! L_gwd_sat is true by default. ! Stuart Webster. ! 5.3 16/10/01 Change to nsigma*sd_orog for the mountain top ! heights. Stuart Webster. ! 6.2 15/08/05 Free format fixes. P.Selwood ! 6.2 21/02/06 Two minor bug fixes - invoked using l_fix_gwsatn ! Stuart Webster ! ! ! Language: fortran 77 + common extensions ! This code is written to umdp3 v6 programming standards. ! Suitable for single column use,rotated grids ! Global variables !*L------------------COMDECK C_G---------------------------------------- ! G IS MEAN ACCEL DUE TO GRAVITY AT EARTH'S SURFACE REAL(KIND=r8), PARAMETER :: G = 9.80665_r8 !*---------------------------------------------------------------------- ! Local constants ! ! Description: This comdeck defines the constants for the 3B and 4A ! versions of the Gravity Wave Drag Code. These are ! tuneable parameters but are unlikely to be changed. ! ! History: ! Version Date Comment ! ------- ---- ------- ! 3.4 18/10/94 Original Version J.R. Mitchell ! 4.3 7/03/97 Remove KAY_LEE (now set in RUNCNST) S.Webster ! 4.5 03/08/98 Add GAMMA_SATN (Used in 06_3B). D. Robinson ! 5.0 14/07/99 Remove redundant switches/variables: only keep ! s-version 06_3B parameters. R Rawlins ! 5.2 15/11/00 Set parameters for the 4A scheme. ! Stuart Webster. ! 5.3 16/10/01 Partition 3B and 4A scheme parameters. ! Stuart Webster. ! 5.3 21/09/01 Set parameters for the spectral (middle ! atmosphere) gravity wave forcing scheme. ! Warner and McIntyre, J. Atm. Sci., 2001, ! Scaife et al., Geophys. Res. Lett., 2000 ! Scaife et al, J. Atm. Sci., 2002 give ! further details. ! Adam Scaife ! 5.4 Alters c_gwave.h include file to change from launch level ! to launch eta to make less model level dependent. ! Adam Scaife ! 5.5 25/02/03 Remove 3B GWD parameter settings. Stuart Webster ! ! 6.2 16/06/05 Move CCL parameter to gw_ussp. Andrew Bushell ! ! ! Number of standard deviations above the mean orography of top ! of sub-grid mountains !REAL(KIND=r8),PARAMETER :: NSIGMA = 2.5_r8 ! SPECTRAL GRAVITY WAVE SCHEME PARAMETERS: ! LMINL = 1/max wavelength at launch ! LSTAR = 1/characteristic wavelength ! ETALAUNCH = eta for model launch REAL(KIND=r8),PARAMETER:: LMINL = 1.0_r8/20000.0_r8 REAL(KIND=r8),PARAMETER:: LSTAR = 1.0_r8/4300.0_r8 REAL(KIND=r8),PARAMETER:: ETALAUNCH = 0.045_r8 ! subroutine arguements; INTEGER,INTENT(in):: kMax ! number of model levels INTEGER,INTENT(in):: points ! number of points INTEGER,INTENT(in):: k_top(points) ! model level at mountain tops INTEGER,INTENT(in):: k_top_max ! max(k_top) ! ! The integers below are set in GWD_CTL2 to points if the diagnostics ! are called and to 1 if not. ! INTEGER,INTENT(in):: points_stress_ud INTEGER,INTENT(in):: points_stress_vd INTEGER,INTENT(in):: points_stress_ud_satn INTEGER,INTENT(in):: points_stress_vd_satn INTEGER,INTENT(in):: points_du_dt_satn INTEGER,INTENT(in):: points_dv_dt_satn REAL(KIND=r8),INTENT(in):: r_rho_levels(points,kMax)! heights above z=0 on rho levels REAL(KIND=r8),INTENT(in):: r_theta_levels(points,kMax)! heights above z=0 on theta levels REAL(KIND=r8),INTENT(in):: rho(points,kMax)! density REAL(KIND=r8),INTENT(in):: rho_s(points) ! surface density (calculated in gwsurf) REAL(KIND=r8),INTENT(in):: theta(points,kMax) ! theta field REAL(KIND=r8),INTENT(in):: u(points,kMax) ! u field REAL(KIND=r8),INTENT(in):: v(points,kMax) ! v field REAL(KIND=r8),INTENT(in):: s_x_stress(points)! linear surface x_stress REAL(KIND=r8),INTENT(in):: s_y_stress(points) ! linear surface y_stress REAL(KIND=r8),INTENT(in):: s_x_orog(points)! 'surface' x_orog REAL(KIND=r8),INTENT(in):: s_y_orog(points)! 'surface' y_orog REAL(KIND=r8),INTENT(in):: sd_orog(points) ! standard deviation of orography REAL(KIND=r8),INTENT(in):: kay ! GWD Constant (m-1) REAL(KIND=r8),INTENT(in):: fr(points) !in low level Froude number REAL(KIND=r8),INTENT(out):: du_dt (points,kMax)! total GWD du/dt REAL(KIND=r8),INTENT(out):: dv_dt (points,kMax) ! total GWD dv/dt ! Diagnostics REAL(KIND=r8),INTENT(out):: du_dt_satn (points_du_dt_satn,kMax) ! Saturation du/dt REAL(KIND=r8),INTENT(out):: dv_dt_satn (points_dv_dt_satn,kMax) ! Saturation dv/dt REAL(KIND=r8),INTENT(out):: stress_ud (points_stress_ud,0:kMax) ! u-stress diagnostic REAL(KIND=r8),INTENT(out):: stress_vd (points_stress_vd,0:kMax) ! v-stress diagnostic REAL(KIND=r8),INTENT(out):: stress_ud_satn (points_stress_ud_satn,0:kMax)! x saturation stress diag REAL(KIND=r8),INTENT(out):: stress_vd_satn (points_stress_vd_satn,0:kMax)! y saturation stress diag LOGICAL,INTENT(in):: l_drag(points) ! whether a point has a ! !non-zero surface stress LOGICAL,INTENT(in):: l_fix_gwsatn ! Switch at vn6.2 - if true then ! ! minor big fixes invoked LOGICAL,INTENT(in):: l_gwd_40km ! Switch at vn6.6 - if true then ! ! GWD not applied above z=40km ! ! These logical switches determine whether a diagnostic is to be ! calculated. ! LOGICAL,INTENT(in):: stress_ud_on LOGICAL,INTENT(in):: stress_ud_p_on LOGICAL,INTENT(in):: stress_vd_on LOGICAL,INTENT(in):: stress_ud_satn_on LOGICAL,INTENT(in):: stress_vd_satn_on LOGICAL,INTENT(in):: du_dt_satn_on LOGICAL,INTENT(in):: du_dt_satn_p_on LOGICAL,INTENT(in):: dv_dt_satn_on !----------------------------------------------------------------- ! LOCAL ARRAYS AND SCALARS !----------------------------------------------------------------- INTEGER :: i,k ! loop counters in routine INTEGER :: kk,kl,ku ! level counters in routine REAL(KIND=r8) :: x_stress(points,2)! x_stresses (layer boundaries) REAL(KIND=r8) :: y_stress(points,2) ! y_stresses (layer boundaries) REAL(KIND=r8) :: x_s_const(points) ! level independent constants for REAL(KIND=r8) :: y_s_const(points) ! calculation of critical stresses. REAL(KIND=r8) :: dzkrho(points) ! depth of level k *rho ! ! (for diagnostic calculations) REAL(KIND=r8) :: dzb ! height difference across rho layer REAL(KIND=r8) :: delta_z ! height difference across theta layer REAL(KIND=r8) :: dzt ! height difference across 2 theta layers REAL(KIND=r8) :: dzu,dzl ! height differences in half layers REAL(KIND=r8) :: ub ! u-wind at u layer boundary REAL(KIND=r8) :: vb ! v-wind at u layer boundary REAL(KIND=r8) :: n ! brunt_vaisala frequency REAL(KIND=r8) :: n_sq ! square of brunt_vaisala frequency REAL(KIND=r8) :: c_x_stress ! critical x_stress REAL(KIND=r8) :: c_y_stress ! critical y_stress REAL(KIND=r8) :: c_stress ! critical stress amplitude REAL(KIND=r8) :: stress ! current stress amplitude REAL(KIND=r8) :: r_stress ! inverse of current stress amplitude REAL(KIND=r8) :: s_stress_sq ! current stress amplitude squared REAL(KIND=r8) :: stress_factor ! 1. or 1.000001 when using more ! ! numerically robust option REAL(KIND=r8) :: speedcalc ! dot product calculation for speed/stress REAL(KIND=r8) :: maxgwdlev ! set to 40km if l_gwd_40km is true ! ! set to 1e10 if l_gwd_40km is false ! ! (so that if test is never true) x_stress=0.0_r8 y_stress=0.0_r8 x_s_const=0.0_r8 y_s_const=0.0_r8 dzkrho=0.0_r8 dzb =0.0_r8 delta_z =0.0_r8 dzt =0.0_r8 dzu =0.0_r8 dzl =0.0_r8 ub =0.0_r8 vb =0.0_r8 n =0.0_r8 n_sq =0.0_r8 c_x_stress =0.0_r8 c_y_stress =0.0_r8 c_stress =0.0_r8 stress =0.0_r8 r_stress =0.0_r8 s_stress_sq =0.0_r8 stress_factor =0.0_r8 speedcalc =0.0_r8 maxgwdlev =0.0_r8 ! function and subroutine calls ! none du_dt=0.0_r8 dv_dt=0.0_r8 du_dt_satn =0.0_r8 dv_dt_satn =0.0_r8 stress_ud=0.0_r8 stress_vd=0.0_r8 stress_ud_satn =0.0_r8 stress_vd_satn =0.0_r8 !------------------------------------------------------------------- ! 1. preliminaries !------------------------------------------------------------------- ! ! Set critical stress_factor according to switch l_fix_gwsatn ! IF ( L_fix_gwsatn ) THEN stress_factor = 1.000001_r8 ELSE stress_factor = 1.0_r8 END IF ! ! Set maxgwdlev according to l_gwd_40km ! IF ( L_gwd_40km ) THEN maxgwdlev = 40000.0_r8 ELSE maxgwdlev = 1.e10_r8 END IF kl=1 ku=2 DO i=1,points x_stress(i,kl) = s_x_stress(i) y_stress(i,kl) = s_y_stress(i) END DO ! ! Calculate x_s_const. ! This includes all terms in the calculation of the critical ! stress that are independent of model level ! DO i=1,points IF ( l_drag(i) .AND. kay /= 0.0_r8 ) THEN s_stress_sq = SQRT( s_x_stress(i) * s_x_stress(i) & & + s_y_stress(i) * s_y_stress(i) ) ! ! 2 is from rho average in following loop over levels ! fr is passed in from gw_satn and scales the critical stress ! so that it is equal to the surface stress if rho, U and N ! were the same as in the surface stress calculation ! r_stress = 1.0_r8 / ( 2.0_r8 * s_stress_sq**3 * nsigma * sd_orog(i) & & * nsigma * sd_orog(i) * fr(i) * fr(i) ) x_s_const(i) = s_x_orog(i) * r_stress y_s_const(i) = s_y_orog(i) * r_stress ELSE x_s_const(i) = 0.0_r8 y_s_const(i) = 0.0_r8 END IF END DO IF( stress_ud_on .OR. stress_ud_p_on ) THEN DO i=1,points stress_ud(i,0) = s_x_stress(i) END DO END IF IF( stress_vd_on ) THEN DO i=1,points stress_vd(i,0) = s_y_stress(i) END DO END IF IF( stress_ud_satn_on ) THEN DO i=1,points stress_ud_satn(i,0) = s_x_stress(i) END DO END IF IF( stress_vd_satn_on ) THEN DO i=1,points stress_vd_satn(i,0) = s_y_stress(i) END DO END IF !------------------------------------------------------------------ ! 2 Loop over levels calculating critical stresses at each theta ! level and reducing the actual stress to that level if it ! exceeds the critical stress. The change in stress ! across the model level is then used to determine the ! deceleration of the wind. ! ! note: k is looping over rho levels !------------------------------------------------------------------ DO k=2,kMax-1 ! min value of k_top is 2 so first check must ! ! be at the top of that level DO i=1,points x_stress(i,ku) = x_stress(i,kl) y_stress(i,ku) = y_stress(i,kl) ! ! First level for wave breaking test is the first ! one above the mountain tops. ! Top level for wavebreaking now either 40km or level below model lid. ! IF ( k>=k_top(i) .AND. r_theta_levels(i,k)<=maxgwdlev ) THEN IF ( ( x_stress(i,kl) /= 0.0_r8) & & .OR.( y_stress(i,kl) /= 0.0_r8) ) THEN dzl = r_theta_levels(i,k ) - r_rho_levels(i,k ) dzu = r_rho_levels (i,k+1) - r_theta_levels(i,k ) ! ub here actually is ub * dzb. Only divide speedcalc by dzb if ! needed in if test below ub = dzu*u(i,k) + dzl*u(i,k+1) vb = dzu*v(i,k) + dzl*v(i,k+1) speedcalc = ub*s_x_stress(i) + vb*s_y_stress(i) ! ! All wave stress deposited if dth/dz < 0 or wind more than ! pi/2 to direction of surface stress ! IF ( theta(i,k+1) <= theta(i,k-1) .OR. & & speedcalc <= 0.0_r8 ) THEN x_stress(i,ku) = 0.0_r8 y_stress(i,ku) = 0.0_r8 ELSE dzt = r_theta_levels(i,k+1) - & & r_theta_levels(i,k-1) n_sq = g * ( theta(i,k+1) - theta(i,k-1) ) & & / ( theta(i,k) * dzt ) n = SQRT( n_sq ) dzb = r_rho_levels(i,k+1) - r_rho_levels(i,k) speedcalc = speedcalc / dzb ! ! rho here is just the average of the full level values. Should ! interpolate in z if this error wasn't so relatively small! ! c_y_stress = kay*(rho(i,k)+rho(i,k+1))* & & stress_factor * speedcalc**3/n c_x_stress = x_s_const(i) * c_y_stress c_y_stress = y_s_const(i) * c_y_stress stress = x_stress(i,kl) * x_stress(i,kl) + & & y_stress(i,kl) * y_stress(i,kl) c_stress = c_x_stress * c_x_stress + & & c_y_stress * c_y_stress IF ( stress > c_stress ) THEN x_stress(i,ku) = x_stress(i,kl) & & *SQRT(c_stress/stress) y_stress(i,ku) = y_stress(i,kl) & & *SQRT(c_stress/stress) ELSE x_stress(i,ku) = x_stress(i,kl) y_stress(i,ku) = y_stress(i,kl) END IF! stress < c_stress END IF ! n_sq > 0 and speedcalc > 0 !------------------------------------------------------------------ ! Calculate drag from vertical stress convergence ! Note that any stress drop at the first level above the mountains ! is applied uniformly from that level down to the ground. !------------------------------------------------------------------ IF( k == k_top(i) ) THEN delta_z = r_theta_levels(i,k) IF ( L_fix_gwsatn ) THEN delta_z = delta_z * rho_s(i)/ rho(i,k) END IF ELSE delta_z = r_theta_levels(i,k) - r_theta_levels(i,k-1) END IF du_dt(i,k) = (x_stress(i,ku) - x_stress(i,kl)) & & /(delta_z*rho(i,k)) dv_dt(i,k) = (y_stress(i,ku) - y_stress(i,kl)) & & /(delta_z*rho(i,k)) !PRINT*, k,du_dt(i,k),dv_dt(i,k),delta_z END IF ! stress x or y ne 0 END IF ! k ge k_top END DO ! diagnostics IF( stress_ud_on ) THEN DO i=1,points IF ( k >= k_top(i) ) THEN stress_ud(i,k) = x_stress(i,ku) END IF END DO END IF IF( stress_vd_on ) THEN DO i=1,points IF ( k >= k_top(i) ) THEN stress_vd(i,k) = y_stress(i,ku) END IF END DO END IF IF( stress_ud_satn_on ) THEN DO i=1,points IF ( k >= k_top(i) ) THEN stress_ud_satn(i,k) = x_stress(i,ku) END IF END DO END IF IF( stress_vd_satn_on ) THEN DO i=1,points IF ( k >= k_top(i) ) THEN stress_vd_satn(i,k) = y_stress(i,ku) END IF END DO END IF IF( du_dt_satn_on ) THEN DO i=1,points du_dt_satn(i,k) = du_dt(i,k) END DO END IF IF( dv_dt_satn_on ) THEN DO i=1,points dv_dt_satn(i,k) = dv_dt(i,k) END DO END IF ! swap storage for lower and upper layers kk=kl kl=ku ku=kk END DO ! end loop kMax ! ! Top level diagnostics - assume now that no drag is applied at ! the top level ! IF( stress_ud_on .OR. stress_ud_p_on ) THEN DO i=1,points stress_ud(i,kMax) = x_stress(i,kl) END DO END IF IF( stress_vd_on ) THEN DO i=1,points stress_vd(i,kMax) = y_stress(i,kl) END DO END IF IF( stress_ud_satn_on ) THEN DO i=1,points stress_ud_satn(i,kMax) = x_stress(i,kl) END DO END IF IF( stress_vd_satn_on ) THEN DO i=1,points stress_vd_satn(i,kMax) = y_stress(i,kl) END DO END IF !-------------------------------------------------------------------- ! 3.1 set drags and stresses below k_top ! can only be Done now because need drag at k_top first !-------------------------------------------------------------------- DO k=1,k_top_max-1 DO i=1,points IF ( k < k_top(i) ) THEN du_dt(i,k) = du_dt(i,k_top(i)) dv_dt(i,k) = dv_dt(i,k_top(i)) END IF END DO END DO ! ! saturation drag diagnostic ! IF( du_dt_satn_on .OR. du_dt_satn_p_on ) THEN DO k=1,k_top_max-1 DO i=1,points IF ( k < k_top(i) ) THEN du_dt_satn(i,k) = du_dt_satn(i,k_top(i)) END IF END DO END DO END IF IF( dv_dt_satn_on ) THEN DO k=1,k_top_max-1 DO i=1,points IF ( k < k_top(i) ) THEN dv_dt_satn(i,k) = dv_dt_satn(i,k_top(i)) END IF END DO END DO END IF ! ! stress diagnostics ! IF( stress_ud_on .OR. stress_ud_p_on .OR. stress_ud_satn_on .OR. & & stress_vd_on .OR. stress_vd_satn_on) THEN DO k=1,k_top_max-1 DO i=1,points IF ( k < k_top(i) ) THEN IF ( k == 1 ) THEN dzkrho(i) = r_theta_levels(i,1)* rho(i,k) ELSE dzkrho(i) = (r_theta_levels(i,k) & & -r_theta_levels(i,k-1))* rho(i,k) END IF END IF END DO IF( stress_ud_on ) THEN DO i=1,points IF ( k < k_top(i) ) THEN stress_ud(i,k) = stress_ud(i,k-1)+du_dt(i,k)*dzkrho(i) END IF END DO END IF IF( stress_ud_satn_on ) THEN DO i=1,points IF ( k < k_top(i) ) THEN stress_ud_satn(i,k)= stress_ud_satn(i,k-1) & & + du_dt(i,k) * dzkrho(i) END IF END DO END IF IF( stress_vd_on ) THEN DO i=1,points IF ( k < k_top(i) ) THEN stress_vd(i,k)= stress_vd(i,k-1) + dv_dt(i,k) *dzkrho(i) END IF END DO END IF IF( stress_vd_satn_on ) THEN DO i=1,points IF ( k < k_top(i) ) THEN stress_vd_satn(i,k)= stress_vd_satn(i,k-1) & & + dv_dt(i,k) * dzkrho(i) END IF END DO END IF END DO ! k=1,k_top_max-1 END IF ! stress_ud_on .or. stress_ud_p_on .or. stess_ud_satn_on .or. ! ! stress_vd_on .or. stess_vd_satn_on RETURN END SUBROUTINE gw_satn !=============================================================================== SUBROUTINE geopotential_t( & piln , pint , pmid , pdel , rpdel , & t , q , rair , gravit , zvir , & zi , zm , ncol ,nCols, kMax) !----------------------------------------------------------------------- ! ! Purpose: ! Compute the geopotential height (above the surface) at the midpoints and ! interfaces using the input temperatures and pressures. ! !----------------------------------------------------------------------- IMPLICIT NONE !------------------------------Arguments-------------------------------- ! ! Input arguments ! INTEGER, INTENT(in) :: ncol ! Number of longitudes INTEGER, INTENT(in) :: nCols INTEGER, INTENT(in) :: kMax REAL(r8), INTENT(in) :: piln (nCols,kMax+1) ! Log interface pressures REAL(r8), INTENT(in) :: pint (nCols,kMax+1) ! Interface pressures REAL(r8), INTENT(in) :: pmid (nCols,kMax) ! Midpoint pressures REAL(r8), INTENT(in) :: pdel (nCols,kMax) ! layer thickness REAL(r8), INTENT(in) :: rpdel(nCols,kMax) ! inverse of layer thickness REAL(r8), INTENT(in) :: t (nCols,kMax) ! temperature REAL(r8), INTENT(in) :: q (nCols,kMax) ! specific humidity REAL(r8), INTENT(in) :: rair ! Gas constant for dry air REAL(r8), INTENT(in) :: gravit ! Acceleration of gravity REAL(r8), INTENT(in) :: zvir ! rh2o/rair - 1 ! Output arguments REAL(r8), INTENT(out) :: zi(nCols,kMax+1) ! Height above surface at interfaces REAL(r8), INTENT(out) :: zm(nCols,kMax) ! Geopotential height at mid level ! !---------------------------Local variables----------------------------- ! LOGICAL :: fvdyn ! finite volume dynamics INTEGER :: i,k ! Lon, level indices REAL(r8) :: hkk(nCols) ! diagonal element of hydrostatic matrix REAL(r8) :: hkl(nCols) ! off-diagonal element REAL(r8) :: rog ! Rair / gravit REAL(r8) :: tv ! virtual temperature REAL(r8) :: tvfac ! Tv/T ! !----------------------------------------------------------------------- ! rog = rair/gravit ! Set dynamics flag fvdyn = .FALSE.!dycore_is ('LR') ! The surface height is zero by definition. DO i = 1,ncol zi(i,kMax+1) = 0.0_r8 END DO ! Compute zi, zm from bottom up. ! Note, zi(i,k) is the interface above zm(i,k) DO k = kMax, 1, -1 ! First set hydrostatic elements consistent with dynamics IF (fvdyn) THEN DO i = 1,ncol hkl(i) = piln(i,k+1) - piln(i,k) hkk(i) = 1.0_r8 - pint(i,k) * hkl(i) * rpdel(i,k) END DO ELSE DO i = 1,ncol hkl(i) = pdel(i,k) / pmid(i,k) hkk(i) = 0.5_r8 * hkl(i) END DO END IF ! Now compute tv, zm, zi DO i = 1,ncol tvfac = 1.0_r8 + zvir * q(i,k) tv = t(i,k) * tvfac zm(i,k) = zi(i,k+1) + rog * tv * hkk(i) zi(i,k) = zi(i,k+1) + rog * tv * hkl(i) END DO END DO RETURN END SUBROUTINE geopotential_t END MODULE GwddSchemeCPTEC