Index: LMDZ6/trunk/libf/phylmd/ener_conserv_mod.f90 =================================================================== --- LMDZ6/trunk/libf/phylmd/ener_conserv_mod.f90 (revision 6131) +++ LMDZ6/trunk/libf/phylmd/ener_conserv_mod.f90 (revision 6132) @@ -36,10 +36,10 @@ USE phys_local_var_mod, ONLY : d_t_eva,d_t_lsc,d_q_eva,d_q_lsc USE phys_local_var_mod, ONLY : d_u_oro,d_v_oro,d_u_lif,d_v_lif -USE phys_local_var_mod, ONLY : du_gwd_hines,dv_gwd_hines,dv_gwd_front,dv_gwd_rando +USE phys_local_var_mod, ONLY : dv_gwd_front,dv_gwd_rando USE phys_state_var_mod, ONLY : du_gwd_front,du_gwd_rando USE phys_output_var_mod, ONLY : bils_ec,bils_ech,bils_tke,bils_kinetic,bils_enthalp,bils_latent,bils_diss USE add_phys_tend_mod, ONLY : fl_cor_ebil USE infotrac_phy, ONLY: nqtot -USE lmdz_gwd_ini, ONLY : ok_gwd_rando, ok_hines +USE lmdz_gwd_ini, ONLY : ok_gwd_rando USE yomcst_mod_h @@ -152,9 +152,5 @@ d_u(:,:)=d_u_vdf(:,:)+d_u_ajs(:,:)+d_u_con(:,:)+d_u_oro(:,:)+d_u_lif(:,:) d_v(:,:)=d_v_vdf(:,:)+d_v_ajs(:,:)+d_v_con(:,:)+d_v_oro(:,:)+d_v_lif(:,:) - IF (ok_hines) THEN - d_u_vdf(:,:)=d_u_vdf(:,:)+du_gwd_hines(:,:) - d_v_vdf(:,:)=d_v_vdf(:,:)+dv_gwd_hines(:,:) - ENDIF - IF (.not. ok_hines .and. ok_gwd_rando) THEN + IF (ok_gwd_rando) THEN d_u_vdf(:,:)=d_u_vdf(:,:)+du_gwd_front(:,:) d_v_vdf(:,:)=d_v_vdf(:,:)+dv_gwd_front(:,:) Index: LMDZ6/trunk/libf/phylmd/lmdz_call_gwd.F90 =================================================================== --- LMDZ6/trunk/libf/phylmd/lmdz_call_gwd.F90 (revision 6131) +++ LMDZ6/trunk/libf/phylmd/lmdz_call_gwd.F90 (revision 6132) @@ -26,10 +26,9 @@ d_u_oro, d_v_oro, d_t_oro, zustr_oro, zvstr_oro, & d_u_lif, d_v_lif, d_t_lif, zustr_lif, zvstr_lif, & - d_u_hines, d_v_hines, d_t_hines, zustr_hines, zvstr_hines, & d_u_front, d_v_front, zustr_front, zvstr_front, & d_u_precip, d_v_precip, zustr_precip, zvstr_precip, & east_gwstress, west_gwstress, aam, torsfc) - USE lmdz_gwd_ini, ONLY: ok_strato, ok_orodr, ok_orolf, ok_hines, ok_gwd_rando + USE lmdz_gwd_ini, ONLY: ok_strato, ok_orodr, ok_orolf, ok_gwd_rando USE lmdz_gwd_ini, ONLY: addtkeoro, smallscales_tkeoro, alphatkeoro USE lmdz_gwd_ini, ONLY: zpmm_orodr_t, zstd_orodr_t, zpmm_orolf_t, nm_oro_t @@ -38,5 +37,4 @@ USE lmdz_gwd_ogwd, ONLY: drag_noro_strato, lift_noro_strato USE lmdz_gwd_ogwd_old, ONLY: drag_noro, lift_noro - USE lmdz_gwd_hines, ONLY: hines_gwd USE lmdz_gwd_front, ONLY: acama_gwd_rando USE lmdz_gwd_precip, ONLY: flott_gwd_rando @@ -111,8 +109,4 @@ REAL, DIMENSION(klon, klev), INTENT(OUT) :: d_t_lif ! t increment due to sso lift [K] - REAL, DIMENSION(klon, klev), INTENT(OUT) :: d_u_hines ! u increment due to Hines param [m/s] - REAL, DIMENSION(klon, klev), INTENT(OUT) :: d_v_hines ! v increment due to Hines param [m/s] - REAL, DIMENSION(klon, klev), INTENT(OUT) :: d_t_hines ! t increment due to Hines param [K] - REAL, DIMENSION(klon), INTENT(OUT) :: zustr_oro ! vertically integrated zonal stress due to sso drag [kg.m/s/m2] REAL, DIMENSION(klon), INTENT(OUT) :: zvstr_oro ! vertically integrated meridional stress due to sso drag [kg.m/s/m2] @@ -120,7 +114,4 @@ REAL, DIMENSION(klon), INTENT(OUT) :: zustr_lif ! vertically integrated zonal stress due to sso lift [kg.m/s/m2] REAL, DIMENSION(klon), INTENT(OUT) :: zvstr_lif ! vertically integrated meridional stress due to sso lift [kg.m/s/m2] - - REAL, DIMENSION(klon), INTENT(OUT) :: zustr_hines ! vertically integrated zonal stress due to Hines param [kg.m/s/m2] - REAL, DIMENSION(klon), INTENT(OUT) :: zvstr_hines ! vertically integrated meridional stress due to Hines param [kg.m/s/m2] REAL, DIMENSION(klon), INTENT(OUT) :: zustr_front ! vertically integrated zonal stress due to fronts gwd param [kg.m/s/m2] @@ -168,14 +159,7 @@ d_t_lif(:, :) = 0. - d_u_hines(:, :) = 0. - d_v_hines(:, :) = 0. - d_t_hines(:, :) = 0. - ! DO NOT set tendencies associated with front and precip to 0 ! as they are inout variables because ! the parameterizations considers some "process memory" - - zustr_hines(:) = 0. - zvstr_hines(:) = 0. zustr_front(:) = 0. @@ -322,39 +306,6 @@ !-------------------------------------- - ! old approach following Hines 1997 - IF (ok_hines) THEN - ! this routine is now out of date. Do not use it except for specific purposes - CALL hines_gwd(klon, klev, phys_tstep, paprs, pplay, latitude_deg, t_seri, & - u_seri, v_seri, zustr_hines, zvstr_hines, d_t_hines, & - d_u_hines, d_v_hines) - - ! Add tendencies - !----------------------------------------------------------------------- - -#ifdef ISO - STOP 'appel a add_phys_tend non traite dans call_gwd' -#else - CALL add_phys_tend(d_u_hines, d_v_hines, d_t_hines, dq0, dql0, & - dqi0, dqbs0, paprs, 'hin', abortphy, flag_inhib_tend, itap, 0) -#endif - - CALL prt_enerbil('hin', itap) - - ! Local diagnostics - !----------------------------------------------------------------------- - DO k = 1, klev - zustr_hines(:) = zustr_hines(:) + d_u_hines(:, k)/phys_tstep & - *(paprs(:, k) - paprs(:, k + 1))/rg - zvstr_hines(:) = zvstr_hines(:) + d_v_hines(:, k)/phys_tstep & - *(paprs(:, k) - paprs(:, k + 1))/rg - END DO - - east_gwstress(:, :) = 0. - west_gwstress(:, :) = 0. - - END IF - ! stochastic approach following De la Camara & Lott 2015 - IF (.NOT. ok_hines .AND. ok_gwd_rando) THEN + IF (ok_gwd_rando) THEN CALL acama_gwd_rando(klon, klev, phys_tstep, pplay, presnivs, latitude_deg, t_seri, u_seri, & Index: LMDZ6/trunk/libf/phylmd/lmdz_gwd_hines.f90 =================================================================== --- LMDZ6/trunk/libf/phylmd/lmdz_gwd_hines.f90 (revision 6131) +++ (revision ) @@ -1,1993 +1,0 @@ -! $Id: hines_gwd.f90 5309 2024-11-01 11:39:44Z abarral $ -!$gpum nogpu -MODULE lmdz_gwd_hines - - IMPLICIT NONE - -CONTAINS - - SUBROUTINE hines_gwd(klon, klev, dtime, paphm1x, papm1x, rlat, tx, ux, vx, & - zustrhi, zvstrhi, d_t_hin, d_u_hin, d_v_hin) - - ! ######################################################################## - ! Parametrization of the momentum flux deposition due to a broad band - ! spectrum of gravity waves, following Hines (1997a,b), as coded by - ! McLANDRESS (1995). Modified by McFARLANE and MANZINI (1995-1997) - ! MAECHAM model stand alone version - ! - ! contact: Francois Lott, flott@lmd.ens.fr - ! - ! References: - ! Hines, C. O., 1997, Doppler-spread parameterization of gravity wave momentum - ! deposit in the middle atmosphere. Part 2: Broad and quasi-monochromatic - ! spectra and implementation, J. Atmos. Solar-Terrest. Phys., 59, 387–400 - ! ######################################################################## - - USE lmdz_gwd_ini, ONLY: RPI, RG, RD, RCPD - - IMPLICIT NONE - - INTEGER nazmth - PARAMETER(nazmth=8) - - ! INPUT ARGUMENTS. - ! ----- ---------- - - ! - 2D - ! PAPHM1 : HALF LEVEL PRESSURE (T-DT) - ! PAPM1 : FULL LEVEL PRESSURE (T-DT) - ! PTM1 : TEMPERATURE (T-DT) - ! PUM1 : ZONAL WIND (T-DT) - ! PVM1 : MERIDIONAL WIND (T-DT) - - ! REFERENCE. - ! ---------- - ! SEE MODEL DOCUMENTATION - - ! AUTHOR. - ! ------- - - ! N. MCFARLANE DKRZ-HAMBURG MAY 1995 - ! STAND ALONE E. MANZINI MPI-HAMBURG FEBRUARY 1997 - - ! BASED ON A COMBINATION OF THE OROGRAPHIC SCHEME BY N.MCFARLANE 1987 - ! AND THE HINES SCHEME AS CODED BY C. MCLANDRESS 1995. - - ! ym INTEGER KLEVM1 - - REAL paphm1(klon, klev + 1), papm1(klon, klev) - REAL ptm1(klon, klev), pum1(klon, klev), pvm1(klon, klev) - REAL prflux(klon) - ! 1 - ! 1 - ! 1 - REAL rlat(klon), coslat(klon) - - REAL th(klon, klev), utendgw(klon, klev), vtendgw(klon, klev), & - pressg(klon), uhs(klon, klev), vhs(klon, klev), zpr(klon) - - ! * VERTICAL POSITIONING ARRAYS. - - REAL sgj(klon, klev), shj(klon, klev), shxkj(klon, klev), dsgj(klon, klev) - - ! * LOGICAL SWITCHES TO CONTROL ROOF DRAG, ENVELOP GW DRAG AND - ! * HINES' DOPPLER SPREADING EXTROWAVE GW DRAG. - ! * LOZPR IS TRUE FOR ZPR ENHANCEMENT - - ! * WORK ARRAYS. - - REAL m_alpha(klon, klev, nazmth), v_alpha(klon, klev, nazmth), & - sigma_alpha(klon, klev, nazmth), sigsqh_alpha(klon, klev, nazmth), & - drag_u(klon, klev), drag_v(klon, klev), flux_u(klon, klev), & - flux_v(klon, klev), heat(klon, klev), diffco(klon, klev), & - bvfreq(klon, klev), density(klon, klev), sigma_t(klon, klev), & - visc_mol(klon, klev), alt(klon, klev), sigsqmcw(klon, klev, nazmth), & - sigmatm(klon, klev), ak_alpha(klon, nazmth), k_alpha(klon, nazmth), & - mmin_alpha(klon, nazmth), i_alpha(klon, nazmth), rmswind(klon), & - bvfbot(klon), densbot(klon) - REAL smoothr1(klon, klev), smoothr2(klon, klev) - REAL sigalpmc(klon, klev, nazmth) - REAL f2mod(klon, klev) - - ! * THES ARE THE INPUT PARAMETERS FOR HINES ROUTINE AND - ! * ARE SPECIFIED IN ROUTINE HINES_SETUP. SINCE THIS IS CALLED - ! * ONLY AT FIRST CALL TO THIS ROUTINE THESE VARIABLES MUST BE SAVED - ! * FOR USE AT SUBSEQUENT CALLS. THIS CAN BE AVOIDED BY CALLING - ! * HINES_SETUP IN MAIN PROGRAM AND PASSING THE PARAMETERS AS - ! * SUBROUTINE ARGUEMENTS. - - REAL rmscon - INTEGER nmessg, iprint, ilrms - INTEGER ifl - - INTEGER naz, icutoff, nsmax, iheatcal - REAL slope, f1, f2, f3, f5, f6, kstar(klon), alt_cutoff, smco - - ! PROVIDED AS INPUT - - INTEGER klon, klev - - REAL dtime - REAL paphm1x(klon, klev + 1), papm1x(klon, klev) - REAL ux(klon, klev), vx(klon, klev), tx(klon, klev) - - ! VARIABLES FOR OUTPUT - - REAL d_t_hin(klon, klev), d_u_hin(klon, klev), d_v_hin(klon, klev) - REAL zustrhi(klon), zvstrhi(klon) - - ! * LOGICAL SWITCHES TO CONTROL PRECIP ENHANCEMENT AND - ! * HINES' DOPPLER SPREADING EXTROWAVE GW DRAG. - ! * LOZPR IS TRUE FOR ZPR ENHANCEMENT - - LOGICAL lozpr, lorms(klon) - - ! LOCAL PARAMETERS TO MAKE THINGS WORK (TEMPORARY VARIABLE) - - REAL rhoh2o, zpcons, rgocp, zlat, dttdsf, ratio, hscal - INTEGER i, j, l, jl, jk, le, lref, lrefp, levbot - - ! DATA PARAMETERS NEEDED, EXPLAINED LATER - - REAL v0, vmin, dmpscal, taufac, hmin, apibt, cpart, fcrit - REAL pcrit, pcons - INTEGER iplev, ierror - - ! PRINT *,' IT IS STARTED HINES GOING ON...' - - ! * COMPUTATIONAL CONSTANTS. - ! ------------- ---------- - - d_t_hin(:, :) = 0. - - rhoh2o = 1000. - zpcons = (1000.*86400.)/rhoh2o - ! ym KLEVM1=KLEV-1 - - DO jl = 1, klon - paphm1(jl, 1) = paphm1x(jl, klev + 1) - DO jk = 1, klev - le = klev + 1 - jk - paphm1(jl, jk + 1) = paphm1x(jl, le) - papm1(jl, jk) = papm1x(jl, le) - ptm1(jl, jk) = tx(jl, le) - pum1(jl, jk) = ux(jl, le) - pvm1(jl, jk) = vx(jl, le) - END DO - END DO - - ! Define constants and arrays needed for the ccc/mam gwd scheme - ! *Constants: - - rgocp = rd/rcpd - lrefp = klev - 1 - lref = klev - 2 - ! 1 - ! 1 *Arrays - ! 1 - DO jk = 1, klev - DO jl = 1, klon - shj(jl, jk) = papm1(jl, jk)/paphm1(jl, klev + 1) - sgj(jl, jk) = papm1(jl, jk)/paphm1(jl, klev + 1) - dsgj(jl, jk) = (paphm1(jl, jk + 1) - paphm1(jl, jk))/paphm1(jl, klev + 1) - shxkj(jl, jk) = (papm1(jl, jk)/paphm1(jl, klev + 1))**rgocp - th(jl, jk) = ptm1(jl, jk) - END DO - END DO - - ! C - DO jl = 1, klon - pressg(jl) = paphm1(jl, klev + 1) - END DO - - DO jl = 1, klon - prflux(jl) = 0.0 - zpr(jl) = zpcons*prflux(jl) - zlat = (rlat(jl)/180.)*rpi - coslat(jl) = cos(zlat) - END DO - - ! /######################################################################### - ! / - ! / - - ! * AUG. 14/95 - C. MCLANDRESS. - ! * SEP. 95 N. MCFARLANE. - - ! * THIS ROUTINE CALCULATES THE HORIZONTAL WIND TENDENCIES - ! * DUE TO MCFARLANE'S OROGRAPHIC GW DRAG SCHEME, HINES' - ! * DOPPLER SPREAD SCHEME FOR "EXTROWAVES" AND ADDS ON - ! * ROOF DRAG. IT IS BASED ON THE ROUTINE GWDFLX8. - - ! * LREFP IS THE INDEX OF THE MODEL LEVEL BELOW THE REFERENCE LEVEL - ! * I/O ARRAYS PASSED FROM MAIN. - ! * (PRESSG = SURFACE PRESSURE) - - ! * CONSTANTS VALUES DEFINED IN DATA STATEMENT ARE : - ! * VMIN = MIMINUM WIND IN THE DIRECTION OF REFERENCE LEVEL - ! * WIND BEFORE WE CONSIDER BREAKING TO HAVE OCCURED. - ! * DMPSCAL = DAMPING TIME FOR GW DRAG IN SECONDS. - ! * TAUFAC = 1/(LENGTH SCALE). - ! * HMIN = MIMINUM ENVELOPE HEIGHT REQUIRED TO PRODUCE GW DRAG. - ! * V0 = VALUE OF WIND THAT APPROXIMATES ZERO. - - DATA vmin/5.0/, v0/1.E-10/, taufac/5.E-6/, hmin/40000./, dmpscal/6.5E+6/, & - apibt/1.5708/, cpart/0.7/, fcrit/1./ - - ! * HINES EXTROWAVE GWD CONSTANTS DEFINED IN DATA STATEMENT ARE: - ! * RMSCON = ROOT MEAN SQUARE GRAVITY WAVE WIND AT LOWEST LEVEL (M/S). - ! * NMESSG = UNIT NUMBER FOR PRINTED MESSAGES. - ! * IPRINT = 1 TO DO PRINT OUT SOME HINES ARRAYS. - ! * IFL = FIRST CALL FLAG TO HINES_SETUP ("SAVE" IT) - ! * PCRIT = CRITICAL VALUE OF ZPR (MM/D) - ! * IPLEV = LEVEL OF APPLICATION OF PRCIT - ! * PCONS = FACTOR OF ZPR ENHANCEMENT - - DATA pcrit/5./, pcons/4.75/ - - iplev = lrefp - 1 - - DATA rmscon/1.00/iprint/2/, nmessg/6/ - DATA ifl/0/ - - lozpr = .FALSE. - - ! ----------------------------------------------------------------------- - - ! * SET ERROR FLAG - - ierror = 0 - - ! * SPECIFY VARIOUS PARAMETERS FOR HINES ROUTINE AT VERY FIRST CALL. - ! * (NOTE THAT ARRAY K_ALPHA IS SPECIFIED SO MAKE SURE THAT - ! * IT IS NOT OVERWRITTEN LATER ON). - - CALL hines_setup(naz, slope, f1, f2, f3, f5, f6, kstar, icutoff, & - alt_cutoff, smco, nsmax, iheatcal, k_alpha, ierror, nmessg, klon, nazmth, & - coslat) - IF (ierror /= 0) GO TO 999 - - ! * START GWD CALCULATIONS. - - lref = lrefp - 1 - - DO j = 1, nazmth - DO l = 1, klev - DO i = 1, klon - sigsqmcw(i, l, j) = 0. - END DO - END DO - END DO - - ! * INITIALIZE NECESSARY ARRAYS. - - DO l = 1, klev - DO i = 1, klon - utendgw(i, l) = 0. - vtendgw(i, l) = 0. - - uhs(i, l) = 0. - vhs(i, l) = 0. - - END DO - END DO - - ! * IF USING HINES SCHEME THEN CALCULATE B V FREQUENCY AT ALL POINTS - ! * AND SMOOTH BVFREQ. - - DO l = 2, klev - DO i = 1, klon - dttdsf = (th(i, l)/shxkj(i, l) - th(i, l - 1)/shxkj(i, l - 1))/ & - (shj(i, l) - shj(i, l - 1)) - dttdsf = min(dttdsf, -5./sgj(i, l)) - bvfreq(i, l) = sqrt(-dttdsf*sgj(i, l)*(sgj(i, l)**rgocp)/rd)*rg/ptm1(i, l & - ) - END DO - END DO - DO l = 1, klev - DO i = 1, klon - IF (l == 1) THEN - bvfreq(i, l) = bvfreq(i, l + 1) - END IF - IF (l > 1) THEN - ratio = 5.*log(sgj(i, l)/sgj(i, l - 1)) - bvfreq(i, l) = (bvfreq(i, l - 1) + ratio*bvfreq(i, l))/(1.+ratio) - END IF - END DO - END DO - - ! * CALCULATE GW DRAG DUE TO HINES' EXTROWAVES - ! * SET MOLECULAR VISCOSITY TO A VERY SMALL VALUE. - ! * IF THE MODEL TOP IS GREATER THAN 100 KM THEN THE ACTUAL - ! * VISCOSITY COEFFICIENT COULD BE SPECIFIED HERE. - - DO l = 1, klev - DO i = 1, klon - visc_mol(i, l) = 1.5E-5 - drag_u(i, l) = 0. - drag_v(i, l) = 0. - flux_u(i, l) = 0. - flux_v(i, l) = 0. - heat(i, l) = 0. - diffco(i, l) = 0. - END DO - END DO - - ! * ALTITUDE AND DENSITY AT BOTTOM. - - DO i = 1, klon - hscal = rd*ptm1(i, klev)/rg - density(i, klev) = sgj(i, klev)*pressg(i)/(rg*hscal) - alt(i, klev) = 0. - END DO - - ! * ALTITUDE AND DENSITY AT REMAINING LEVELS. - - DO l = klev - 1, 1, -1 - DO i = 1, klon - hscal = rd*ptm1(i, l)/rg - alt(i, l) = alt(i, l + 1) + hscal*dsgj(i, l)/sgj(i, l) - density(i, l) = sgj(i, l)*pressg(i)/(rg*hscal) - END DO - END DO - - ! * INITIALIZE SWITCHES FOR HINES GWD CALCULATION - - ilrms = 0 - - DO i = 1, klon - lorms(i) = .FALSE. - END DO - - ! * DEFILE BOTTOM LAUNCH LEVEL - - levbot = iplev - - ! * BACKGROUND WIND MINUS VALUE AT BOTTOM LAUNCH LEVEL. - - DO l = 1, levbot - DO i = 1, klon - uhs(i, l) = pum1(i, l) - pum1(i, levbot) - vhs(i, l) = pvm1(i, l) - pvm1(i, levbot) - END DO - END DO - - ! * SPECIFY ROOT MEAN SQUARE WIND AT BOTTOM LAUNCH LEVEL. - - DO i = 1, klon - rmswind(i) = rmscon - END DO - - IF (lozpr) THEN - DO i = 1, klon - IF (zpr(i) > pcrit) THEN - rmswind(i) = rmscon + ((zpr(i) - pcrit)/zpr(i))*pcons - END IF - END DO - END IF - - DO i = 1, klon - IF (rmswind(i) > 0.0) THEN - ilrms = ilrms + 1 - lorms(i) = .TRUE. - END IF - END DO - - ! * CALCULATE GWD (NOTE THAT DIFFUSION COEFFICIENT AND - ! * HEATING RATE ONLY CALCULATED IF IHEATCAL = 1). - - IF (ilrms > 0) THEN - - CALL hines_extro0(drag_u, drag_v, heat, diffco, flux_u, flux_v, uhs, vhs, & - bvfreq, density, visc_mol, alt, rmswind, k_alpha, m_alpha, v_alpha, & - sigma_alpha, sigsqh_alpha, ak_alpha, mmin_alpha, i_alpha, sigma_t, & - densbot, bvfbot, 1, iheatcal, icutoff, iprint, nsmax, smco, alt_cutoff, & - kstar, slope, f1, f2, f3, f5, f6, naz, sigsqmcw, sigmatm, 1, klon, & - 1, levbot, klon, klev, nazmth, lorms, smoothr1, smoothr2, sigalpmc, & - f2mod) - - ! * ADD ON HINES' GWD TENDENCIES TO OROGRAPHIC TENDENCIES AND - ! * APPLY HINES' GW DRAG ON (UROW,VROW) WORK ARRAYS. - - DO l = 1, klev - DO i = 1, klon - utendgw(i, l) = utendgw(i, l) + drag_u(i, l) - vtendgw(i, l) = vtendgw(i, l) + drag_v(i, l) - END DO - END DO - - ! * END OF HINES CALCULATIONS. - - END IF - - ! ----------------------------------------------------------------------- - - DO jl = 1, klon - zustrhi(jl) = flux_u(jl, 1) - zvstrhi(jl) = flux_v(jl, 1) - DO jk = 1, klev - le = klev - jk + 1 - d_u_hin(jl, jk) = utendgw(jl, le)*dtime - d_v_hin(jl, jk) = vtendgw(jl, le)*dtime - END DO - END DO - - ! PRINT *,'UTENDGW:',UTENDGW - - ! PRINT *,' HINES HAS BEEN COMPLETED (LONG ISNT IT...)' - - RETURN -999 CONTINUE - - ! * IF ERROR DETECTED THEN ABORT. - - WRITE (nmessg, 6000) - WRITE (nmessg, 6010) ierror -6000 FORMAT(/' EXECUTION ABORTED IN GWDOREXV') -6010 FORMAT(' ERROR FLAG =', I4) - - RETURN - END SUBROUTINE hines_gwd -! / -! / - - SUBROUTINE hines_extro0(drag_u, drag_v, heat, diffco, flux_u, flux_v, vel_u, & - vel_v, bvfreq, density, visc_mol, alt, rmswind, k_alpha, m_alpha, & - v_alpha, sigma_alpha, sigsqh_alpha, ak_alpha, mmin_alpha, i_alpha, & - sigma_t, densb, bvfb, iorder, iheatcal, icutoff, iprint, nsmax, smco, & - alt_cutoff, kstar, slope, f1, f2, f3, f5, f6, naz, sigsqmcw, sigmatm, & - il1, il2, lev1, lev2, klons, klevs, nazmth, lorms, smoothr1, smoothr2, & - sigalpmc, f2mod) - - IMPLICIT NONE - - ! Main routine for Hines' "extrowave" gravity wave parameterization based - ! on Hines' Doppler spread theory. This routine calculates zonal - ! and meridional components of gravity wave drag, heating rates - ! and diffusion coefficient on a longitude by altitude grid. - ! No "mythical" lower boundary region calculation is made so it - ! is assumed that lowest level winds are weak (i.e, approximately zero). - - ! Aug. 13/95 - C. McLandress - ! SEPT. /95 - N.McFarlane - - ! Modifications: - - ! Output arguements: - - ! * DRAG_U = zonal component of gravity wave drag (m/s^2). - ! * DRAG_V = meridional component of gravity wave drag (m/s^2). - ! * HEAT = gravity wave heating (K/sec). - ! * DIFFCO = diffusion coefficient (m^2/sec) - ! * FLUX_U = zonal component of vertical momentum flux (Pascals) - ! * FLUX_V = meridional component of vertical momentum flux (Pascals) - - ! Input arguements: - - ! * VEL_U = background zonal wind component (m/s). - ! * VEL_V = background meridional wind component (m/s). - ! * BVFREQ = background Brunt Vassala frequency (radians/sec). - ! * DENSITY = background density (kg/m^3) - ! * VISC_MOL = molecular viscosity (m^2/s) - ! * ALT = altitude of momentum, density, buoyancy levels (m) - ! * (NOTE: levels ordered so that ALT(I,1) > ALT(I,2), etc.) - ! * RMSWIND = root mean square gravity wave wind at lowest level (m/s). - ! * K_ALPHA = horizontal wavenumber of each azimuth (1/m). - ! * IORDER = 1 means vertical levels are indexed from top down - ! * (i.e., highest level indexed 1 and lowest level klevS); - ! * .NE. 1 highest level is index klevS. - ! * IHEATCAL = 1 to calculate heating rates and diffusion coefficient. - ! * IPRINT = 1 to print out various arrays. - ! * ICUTOFF = 1 to exponentially damp GWD, heating and diffusion - ! * arrays above ALT_CUTOFF; otherwise arrays not modified. - ! * ALT_CUTOFF = altitude in meters above which exponential decay applied. - ! * SMCO = smoothing factor used to smooth cutoff vertical - ! * wavenumbers and total rms winds in vertical direction - ! * before calculating drag or heating - ! * (SMCO >= 1 ==> 1:SMCO:1 stencil used). - ! * NSMAX = number of times smoother applied ( >= 1), - ! * = 0 means no smoothing performed. - ! * KSTAR = typical gravity wave horizontal wavenumber (1/m). - ! * SLOPE = slope of incident vertical wavenumber spectrum - ! * (SLOPE must equal 1., 1.5 or 2.). - ! * F1 to F6 = Hines's fudge factors (F4 not needed since used for - ! * vertical flux of vertical momentum). - ! * NAZ = actual number of horizontal azimuths used. - ! * IL1 = first longitudinal index to use (IL1 >= 1). - ! * IL2 = last longitudinal index to use (IL1 <= IL2 <= klonS). - ! * LEV1 = index of first level for drag calculation. - ! * LEV2 = index of last level for drag calculation - ! * (i.e., LEV1 < LEV2 <= klevS). - ! * klonS = number of longitudes. - ! * klevS = number of vertical levels. - ! * NAZMTH = azimuthal array dimension (NAZMTH >= NAZ). - - ! Work arrays. - - ! * M_ALPHA = cutoff vertical wavenumber (1/m). - ! * V_ALPHA = wind component at each azimuth (m/s) and if IHEATCAL=1 - ! * holds vertical derivative of cutoff wavenumber. - ! * SIGMA_ALPHA = total rms wind in each azimuth (m/s). - ! * SIGSQH_ALPHA = portion of wind variance from waves having wave - ! * normals in the alpha azimuth (m/s). - ! * SIGMA_T = total rms horizontal wind (m/s). - ! * AK_ALPHA = spectral amplitude factor at each azimuth - ! * (i.e.,{AjKj}) in m^4/s^2. - ! * I_ALPHA = Hines' integral. - ! * MMIN_ALPHA = minimum value of cutoff wavenumber. - ! * DENSB = background density at bottom level. - ! * BVFB = buoyancy frequency at bottom level and - ! * work array for ICUTOFF = 1. - - ! * LORMS = .TRUE. for drag computation - - INTEGER naz, klons, klevs, nazmth, il1, il2, lev1, lev2 - INTEGER icutoff, nsmax, iorder, iheatcal, iprint - REAL kstar(klons), f1, f2, f3, f5, f6, slope - REAL alt_cutoff, smco - REAL drag_u(klons, klevs), drag_v(klons, klevs) - REAL heat(klons, klevs), diffco(klons, klevs) - REAL flux_u(klons, klevs), flux_v(klons, klevs) - REAL vel_u(klons, klevs), vel_v(klons, klevs) - REAL bvfreq(klons, klevs), density(klons, klevs) - REAL visc_mol(klons, klevs), alt(klons, klevs) - REAL rmswind(klons), bvfb(klons), densb(klons) - REAL sigma_t(klons, klevs), sigsqmcw(klons, klevs, nazmth) - REAL sigma_alpha(klons, klevs, nazmth), sigmatm(klons, klevs) - REAL sigsqh_alpha(klons, klevs, nazmth) - REAL m_alpha(klons, klevs, nazmth), v_alpha(klons, klevs, nazmth) - REAL ak_alpha(klons, nazmth), k_alpha(klons, nazmth) - REAL mmin_alpha(klons, nazmth), i_alpha(klons, nazmth) - REAL smoothr1(klons, klevs), smoothr2(klons, klevs) - REAL sigalpmc(klons, klevs, nazmth) - REAL f2mod(klons, klevs) - - LOGICAL lorms(klons) - - ! Internal variables. - - INTEGER levbot, levtop, i, n, l, lev1p, lev2m - INTEGER ilprt1, ilprt2 - ! ----------------------------------------------------------------------- - - ! PRINT *,' IN HINES_EXTRO0' - lev1p = lev1 + 1 - lev2m = lev2 - 1 - - ! Index of lowest altitude level (bottom of drag calculation). - - levbot = lev2 - levtop = lev1 - IF (iorder /= 1) THEN - WRITE (6, 1) -1 FORMAT(2X, ' error: IORDER NOT ONE! ') - END IF - - ! Buoyancy and density at bottom level. - - DO i = il1, il2 - bvfb(i) = bvfreq(i, levbot) - densb(i) = density(i, levbot) - END DO - - ! initialize some variables - - DO n = 1, naz - DO l = lev1, lev2 - DO i = il1, il2 - m_alpha(i, l, n) = 0.0 - END DO - END DO - END DO - DO l = lev1, lev2 - DO i = il1, il2 - sigma_t(i, l) = 0.0 - END DO - END DO - DO n = 1, naz - DO i = il1, il2 - i_alpha(i, n) = 0.0 - END DO - END DO - - ! Compute azimuthal wind components from zonal and meridional winds. - - CALL hines_wind(v_alpha, vel_u, vel_v, naz, il1, il2, lev1, lev2, klons, & - klevs, nazmth) - - ! Calculate cutoff vertical wavenumber and velocity variances. - - CALL hines_wavnum(m_alpha, sigma_alpha, sigsqh_alpha, sigma_t, ak_alpha, & - v_alpha, visc_mol, density, densb, bvfreq, bvfb, rmswind, i_alpha, & - mmin_alpha, kstar, slope, f1, f2, f3, naz, levbot, levtop, il1, il2, & - klons, klevs, nazmth, sigsqmcw, sigmatm, lorms, sigalpmc, f2mod) - ! Smooth cutoff wavenumbers and total rms velocity in the vertical - ! direction NSMAX times, using FLUX_U as temporary work array. - - IF (nsmax > 0) THEN - DO n = 1, naz - DO l = lev1, lev2 - DO i = il1, il2 - smoothr1(i, l) = m_alpha(i, l, n) - END DO - END DO - CALL vert_smooth(smoothr1, smoothr2, smco, nsmax, il1, il2, lev1, lev2, & - klons, klevs) - DO l = lev1, lev2 - DO i = il1, il2 - m_alpha(i, l, n) = smoothr1(i, l) - END DO - END DO - END DO - CALL vert_smooth(sigma_t, smoothr2, smco, nsmax, il1, il2, lev1, lev2, & - klons, klevs) - END IF - - ! Calculate zonal and meridional components of the - ! momentum flux and drag. - - CALL hines_flux(flux_u, flux_v, drag_u, drag_v, alt, density, densb, & - m_alpha, ak_alpha, k_alpha, slope, naz, il1, il2, lev1, lev2, klons, & - klevs, nazmth, lorms) - - ! Cutoff drag above ALT_CUTOFF, using BVFB as temporary work array. - - IF (icutoff == 1) THEN - CALL hines_exp(drag_u, bvfb, alt, alt_cutoff, iorder, il1, il2, lev1, & - lev2, klons, klevs) - CALL hines_exp(drag_v, bvfb, alt, alt_cutoff, iorder, il1, il2, lev1, & - lev2, klons, klevs) - END IF - - ! Print out various arrays for diagnostic purposes. - - IF (iprint == 1) THEN - ilprt1 = 15 - ilprt2 = 16 - CALL hines_print(flux_u, flux_v, drag_u, drag_v, alt, sigma_t, & - sigma_alpha, v_alpha, m_alpha, 1, 1, 6, ilprt1, ilprt2, lev1, lev2, & - naz, klons, klevs, nazmth) - END IF - - ! If not calculating heating rate and diffusion coefficient then finished. - - IF (iheatcal /= 1) RETURN - - ! Calculate vertical derivative of cutoff wavenumber (store - ! in array V_ALPHA) using centered differences at interior gridpoints - ! and one-sided differences at first and last levels. - - DO n = 1, naz - DO l = lev1p, lev2m - DO i = il1, il2 - v_alpha(i, l, n) = (m_alpha(i, l + 1, n) - m_alpha(i, l - 1, n))/ & - (alt(i, l + 1) - alt(i, l - 1)) - END DO - END DO - DO i = il1, il2 - v_alpha(i, lev1, n) = (m_alpha(i, lev1p, n) - m_alpha(i, lev1, n))/ & - (alt(i, lev1p) - alt(i, lev1)) - END DO - DO i = il1, il2 - v_alpha(i, lev2, n) = (m_alpha(i, lev2, n) - m_alpha(i, lev2m, n))/ & - (alt(i, lev2) - alt(i, lev2m)) - END DO - END DO - - ! Heating rate and diffusion coefficient. - - CALL hines_heat(heat, diffco, m_alpha, v_alpha, ak_alpha, k_alpha, bvfreq, & - density, densb, sigma_t, visc_mol, kstar, slope, f2, f3, f5, f6, naz, & - il1, il2, lev1, lev2, klons, klevs, nazmth) - - ! Finished. - - RETURN - ! ----------------------------------------------------------------------- - END SUBROUTINE hines_extro0 - - SUBROUTINE hines_wavnum(m_alpha, sigma_alpha, sigsqh_alpha, sigma_t, & - ak_alpha, v_alpha, visc_mol, density, densb, bvfreq, bvfb, rms_wind, & - i_alpha, mmin_alpha, kstar, slope, f1, f2, f3, naz, levbot, levtop, il1, & - il2, klons, klevs, nazmth, sigsqmcw, sigmatm, lorms, sigalpmc, f2mod) - IMPLICIT NONE - ! This routine calculates the cutoff vertical wavenumber and velocity - ! variances on a longitude by altitude grid for the Hines' Doppler - ! spread gravity wave drag parameterization scheme. - ! NOTE: (1) only values of four or eight can be used for # azimuths (NAZ). - ! (2) only values of 1.0, 1.5 or 2.0 can be used for slope (SLOPE). - - ! Aug. 10/95 - C. McLandress - - ! Output arguements: - - ! * M_ALPHA = cutoff wavenumber at each azimuth (1/m). - ! * SIGMA_ALPHA = total rms wind in each azimuth (m/s). - ! * SIGSQH_ALPHA = portion of wind variance from waves having wave - ! * normals in the alpha azimuth (m/s). - ! * SIGMA_T = total rms horizontal wind (m/s). - ! * AK_ALPHA = spectral amplitude factor at each azimuth - ! * (i.e.,{AjKj}) in m^4/s^2. - - ! Input arguements: - - ! * V_ALPHA = wind component at each azimuth (m/s). - ! * VISC_MOL = molecular viscosity (m^2/s) - ! * DENSITY = background density (kg/m^3). - ! * DENSB = background density at model bottom (kg/m^3). - ! * BVFREQ = background Brunt Vassala frequency (radians/sec). - ! * BVFB = background Brunt Vassala frequency at model bottom. - ! * RMS_WIND = root mean square gravity wave wind at lowest level (m/s). - ! * KSTAR = typical gravity wave horizontal wavenumber (1/m). - ! * SLOPE = slope of incident vertical wavenumber spectrum - ! * (SLOPE = 1., 1.5 or 2.). - ! * F1,F2,F3 = Hines's fudge factors. - ! * NAZ = actual number of horizontal azimuths used (4 or 8). - ! * LEVBOT = index of lowest vertical level. - ! * LEVTOP = index of highest vertical level - ! * (NOTE: if LEVTOP < LEVBOT then level index - ! * increases from top down). - ! * IL1 = first longitudinal index to use (IL1 >= 1). - ! * IL2 = last longitudinal index to use (IL1 <= IL2 <= klonS). - ! * klonS = number of longitudes. - ! * klevS = number of vertical levels. - ! * NAZMTH = azimuthal array dimension (NAZMTH >= NAZ). - - ! * LORMS = .TRUE. for drag computation - - ! Input work arrays: - - ! * I_ALPHA = Hines' integral at a single level. - ! * MMIN_ALPHA = minimum value of cutoff wavenumber. - - INTEGER naz, levbot, levtop, il1, il2, klons, klevs, nazmth - REAL slope, kstar(klons), f1, f2, f3, f2mfac - REAL m_alpha(klons, klevs, nazmth) - REAL sigma_alpha(klons, klevs, nazmth) - REAL sigalpmc(klons, klevs, nazmth) - REAL sigsqh_alpha(klons, klevs, nazmth) - REAL sigsqmcw(klons, klevs, nazmth) - REAL sigma_t(klons, klevs) - REAL sigmatm(klons, klevs) - REAL ak_alpha(klons, nazmth) - REAL v_alpha(klons, klevs, nazmth) - REAL visc_mol(klons, klevs) - REAL f2mod(klons, klevs) - REAL density(klons, klevs), densb(klons) - REAL bvfreq(klons, klevs), bvfb(klons), rms_wind(klons) - REAL i_alpha(klons, nazmth), mmin_alpha(klons, nazmth) - - LOGICAL lorms(klons) - - ! Internal variables. - - INTEGER i, l, n, lstart, lend, lincr, lbelow - REAL m_sub_m_turb, m_sub_m_mol, m_trial - REAL visc, visc_min, azfac, sp1 - - ! c REAL N_OVER_M(1000), SIGFAC(1000) - - REAL n_over_m(klons), sigfac(klons) - DATA visc_min/1.E-10/ - ! ----------------------------------------------------------------------- - - ! PRINT *,'IN HINES_WAVNUM' - sp1 = slope + 1. - - ! Indices of levels to process. - - IF (levbot > levtop) THEN - lstart = levbot - 1 - lend = levtop - lincr = -1 - ELSE - WRITE (6, 1) -1 FORMAT(2X, ' error: IORDER NOT ONE! ') - END IF - - ! Use horizontal isotropy to calculate azimuthal variances at bottom level. - - azfac = 1./real(naz) - DO n = 1, naz - DO i = il1, il2 - sigsqh_alpha(i, levbot, n) = azfac*rms_wind(i)**2 - END DO - END DO - - ! Velocity variances at bottom level. - - CALL hines_sigma(sigma_t, sigma_alpha, sigsqh_alpha, naz, levbot, il1, il2, & - klons, klevs, nazmth) - - CALL hines_sigma(sigmatm, sigalpmc, sigsqmcw, naz, levbot, il1, il2, klons, & - klevs, nazmth) - - ! Calculate cutoff wavenumber and spectral amplitude factor - ! at bottom level where it is assumed that background winds vanish - ! and also initialize minimum value of cutoff wavnumber. - - DO n = 1, naz - DO i = il1, il2 - IF (lorms(i)) THEN - m_alpha(i, levbot, n) = bvfb(i)/(f1*sigma_alpha(i, levbot, n) + f2* & - sigma_t(i, levbot)) - ak_alpha(i, n) = sigsqh_alpha(i, levbot, n)/ & - (m_alpha(i, levbot, n)**sp1/sp1) - mmin_alpha(i, n) = m_alpha(i, levbot, n) - END IF - END DO - END DO - - ! Calculate quantities from the bottom upwards, - ! starting one level above bottom. - - DO l = lstart, lend, lincr - - ! Level beneath present level. - - lbelow = l - lincr - - ! Calculate N/m_M where m_M is maximum permissible value of the vertical - ! wavenumber (i.e., m > m_M are obliterated) and N is buoyancy frequency. - ! m_M is taken as the smaller of the instability-induced - ! wavenumber (M_SUB_M_TURB) and that imposed by molecular viscosity - ! (M_SUB_M_MOL). Since variance at this level is not yet known - ! use value at level below. - - DO i = il1, il2 - IF (lorms(i)) THEN - - f2mfac = sigmatm(i, lbelow)**2 - f2mod(i, lbelow) = 1.+2.*f2mfac/(f2mfac + sigma_t(i, lbelow)**2) - - visc = amax1(visc_mol(i, l), visc_min) - m_sub_m_turb = bvfreq(i, l)/(f2*f2mod(i, lbelow)*sigma_t(i, lbelow)) - m_sub_m_mol = (bvfreq(i, l)*kstar(i)/visc)**0.33333333/f3 - IF (m_sub_m_turb < m_sub_m_mol) THEN - n_over_m(i) = f2*f2mod(i, lbelow)*sigma_t(i, lbelow) - ELSE - n_over_m(i) = bvfreq(i, l)/m_sub_m_mol - END IF - END IF - END DO - - ! Calculate cutoff wavenumber at this level. - - DO n = 1, naz - DO i = il1, il2 - IF (lorms(i)) THEN - - ! Calculate trial value (since variance at this level is not yet - ! known - ! use value at level below). If trial value is negative or if it - ! exceeds - ! minimum value (not permitted) then set it to minimum value. - - m_trial = bvfb(i)/(f1*(sigma_alpha(i, lbelow, n) + sigalpmc(i, lbelow, & - n)) + n_over_m(i) + v_alpha(i, l, n)) - IF (m_trial <= 0. .OR. m_trial > mmin_alpha(i, n)) THEN - m_trial = mmin_alpha(i, n) - END IF - m_alpha(i, l, n) = m_trial - - ! Reset minimum value of cutoff wavenumber if necessary. - - IF (m_alpha(i, l, n) < mmin_alpha(i, n)) THEN - mmin_alpha(i, n) = m_alpha(i, l, n) - END IF - - END IF - END DO - END DO - - ! Calculate the Hines integral at this level. - - CALL hines_intgrl(i_alpha, v_alpha, m_alpha, bvfb, slope, naz, l, il1, & - il2, klons, klevs, nazmth, lorms) - - ! Calculate the velocity variances at this level. - - DO i = il1, il2 - sigfac(i) = densb(i)/density(i, l)*bvfreq(i, l)/bvfb(i) - END DO - DO n = 1, naz - DO i = il1, il2 - sigsqh_alpha(i, l, n) = sigfac(i)*ak_alpha(i, n)*i_alpha(i, n) - END DO - END DO - CALL hines_sigma(sigma_t, sigma_alpha, sigsqh_alpha, naz, l, il1, il2, & - klons, klevs, nazmth) - - CALL hines_sigma(sigmatm, sigalpmc, sigsqmcw, naz, l, il1, il2, klons, & - klevs, nazmth) - - ! End of level loop. - - END DO - - RETURN - ! ----------------------------------------------------------------------- - END SUBROUTINE hines_wavnum - - SUBROUTINE hines_wind(v_alpha, vel_u, vel_v, naz, il1, il2, lev1, lev2, & - klons, klevs, nazmth) - IMPLICIT NONE - ! This routine calculates the azimuthal horizontal background wind - ! components - ! on a longitude by altitude grid for the case of 4 or 8 azimuths for - ! the Hines' Doppler spread GWD parameterization scheme. - - ! Aug. 7/95 - C. McLandress - - ! Output arguement: - - ! * V_ALPHA = background wind component at each azimuth (m/s). - ! * (note: first azimuth is in eastward direction - ! * and rotate in counterclockwise direction.) - - ! Input arguements: - - ! * VEL_U = background zonal wind component (m/s). - ! * VEL_V = background meridional wind component (m/s). - ! * NAZ = actual number of horizontal azimuths used (must be 4 or 8). - ! * IL1 = first longitudinal index to use (IL1 >= 1). - ! * IL2 = last longitudinal index to use (IL1 <= IL2 <= klonS). - ! * LEV1 = first altitude level to use (LEV1 >=1). - ! * LEV2 = last altitude level to use (LEV1 < LEV2 <= klevS). - ! * klonS = number of longitudes. - ! * klevS = number of vertical levels. - ! * NAZMTH = azimuthal array dimension (NAZMTH >= NAZ). - - ! Constants in DATA statements. - - ! * COS45 = cosine of 45 degrees. - ! * UMIN = minimum allowable value for zonal or meridional - ! * wind component (m/s). - - ! Subroutine arguements. - - INTEGER naz, il1, il2, lev1, lev2 - INTEGER klons, klevs, nazmth - REAL v_alpha(klons, klevs, nazmth) - REAL vel_u(klons, klevs), vel_v(klons, klevs) - - ! Internal variables. - - INTEGER i, l - REAL u, v, cos45, umin - - DATA cos45/0.7071068/ - DATA umin/0.001/ - ! ----------------------------------------------------------------------- - - ! Case with 4 azimuths. - - ! PRINT *,'IN HINES_WIND' - IF (naz == 4) THEN - DO l = lev1, lev2 - DO i = il1, il2 - u = vel_u(i, l) - v = vel_v(i, l) - IF (abs(u) < umin) u = umin - IF (abs(v) < umin) v = umin - v_alpha(i, l, 1) = u - v_alpha(i, l, 2) = v - v_alpha(i, l, 3) = -u - v_alpha(i, l, 4) = -v - END DO - END DO - END IF - - ! Case with 8 azimuths. - - IF (naz == 8) THEN - DO l = lev1, lev2 - DO i = il1, il2 - u = vel_u(i, l) - v = vel_v(i, l) - IF (abs(u) < umin) u = umin - IF (abs(v) < umin) v = umin - v_alpha(i, l, 1) = u - v_alpha(i, l, 2) = cos45*(v + u) - v_alpha(i, l, 3) = v - v_alpha(i, l, 4) = cos45*(v - u) - v_alpha(i, l, 5) = -u - v_alpha(i, l, 6) = -v_alpha(i, l, 2) - v_alpha(i, l, 7) = -v - v_alpha(i, l, 8) = -v_alpha(i, l, 4) - END DO - END DO - END IF - - RETURN - ! ----------------------------------------------------------------------- - END SUBROUTINE hines_wind - - SUBROUTINE hines_flux(flux_u, flux_v, drag_u, drag_v, alt, density, densb, & - m_alpha, ak_alpha, k_alpha, slope, naz, il1, il2, lev1, lev2, klons, & - klevs, nazmth, lorms) - - IMPLICIT NONE - ! Calculate zonal and meridional components of the vertical flux - ! of horizontal momentum and corresponding wave drag (force per unit mass) - ! on a longitude by altitude grid for the Hines' Doppler spread - ! GWD parameterization scheme. - ! NOTE: only 4 or 8 azimuths can be used. - - ! Aug. 6/95 - C. McLandress - - ! Output arguements: - - ! * FLUX_U = zonal component of vertical momentum flux (Pascals) - ! * FLUX_V = meridional component of vertical momentum flux (Pascals) - ! * DRAG_U = zonal component of drag (m/s^2). - ! * DRAG_V = meridional component of drag (m/s^2). - - ! Input arguements: - - ! * ALT = altitudes (m). - ! * DENSITY = background density (kg/m^3). - ! * DENSB = background density at bottom level (kg/m^3). - ! * M_ALPHA = cutoff vertical wavenumber (1/m). - ! * AK_ALPHA = spectral amplitude factor (i.e., {AjKj} in m^4/s^2). - ! * K_ALPHA = horizontal wavenumber (1/m). - ! * SLOPE = slope of incident vertical wavenumber spectrum. - ! * NAZ = actual number of horizontal azimuths used (must be 4 or 8). - ! * IL1 = first longitudinal index to use (IL1 >= 1). - ! * IL2 = last longitudinal index to use (IL1 <= IL2 <= klonS). - ! * LEV1 = first altitude level to use (LEV1 >=1). - ! * LEV2 = last altitude level to use (LEV1 < LEV2 <= klevS). - ! * klonS = number of longitudes. - ! * klevS = number of vertical levels. - ! * NAZMTH = azimuthal array dimension (NAZMTH >= NAZ). - - ! * LORMS = .TRUE. for drag computation - - ! Constant in DATA statement. - - ! * COS45 = cosine of 45 degrees. - - ! Subroutine arguements. - - INTEGER naz, il1, il2, lev1, lev2 - INTEGER klons, klevs, nazmth - REAL slope - REAL flux_u(klons, klevs), flux_v(klons, klevs) - REAL drag_u(klons, klevs), drag_v(klons, klevs) - REAL alt(klons, klevs), density(klons, klevs), densb(klons) - REAL m_alpha(klons, klevs, nazmth) - REAL ak_alpha(klons, nazmth), k_alpha(klons, nazmth) - - LOGICAL lorms(klons) - - ! Internal variables. - - INTEGER i, l, lev1p, lev2m, lev2p - REAL cos45, prod2, prod4, prod6, prod8, dendz, dendz2 - DATA cos45/0.7071068/ - ! ----------------------------------------------------------------------- - - lev1p = lev1 + 1 - lev2m = lev2 - 1 - lev2p = lev2 + 1 - - ! Sum over azimuths for case where SLOPE = 1. - - IF (slope == 1.) THEN - - ! Case with 4 azimuths. - - IF (naz == 4) THEN - DO l = lev1, lev2 - DO i = il1, il2 - flux_u(i, l) = ak_alpha(i, 1)*k_alpha(i, 1)*m_alpha(i, l, 1) - & - ak_alpha(i, 3)*k_alpha(i, 3)*m_alpha(i, l, 3) - flux_v(i, l) = ak_alpha(i, 2)*k_alpha(i, 2)*m_alpha(i, l, 2) - & - ak_alpha(i, 4)*k_alpha(i, 4)*m_alpha(i, l, 4) - END DO - END DO - END IF - - ! Case with 8 azimuths. - - IF (naz == 8) THEN - DO l = lev1, lev2 - DO i = il1, il2 - prod2 = ak_alpha(i, 2)*k_alpha(i, 2)*m_alpha(i, l, 2) - prod4 = ak_alpha(i, 4)*k_alpha(i, 4)*m_alpha(i, l, 4) - prod6 = ak_alpha(i, 6)*k_alpha(i, 6)*m_alpha(i, l, 6) - prod8 = ak_alpha(i, 8)*k_alpha(i, 8)*m_alpha(i, l, 8) - flux_u(i, l) = ak_alpha(i, 1)*k_alpha(i, 1)*m_alpha(i, l, 1) - & - ak_alpha(i, 5)*k_alpha(i, 5)*m_alpha(i, l, 5) + & - cos45*(prod2 - prod4 - prod6 + prod8) - flux_v(i, l) = ak_alpha(i, 3)*k_alpha(i, 3)*m_alpha(i, l, 3) - & - ak_alpha(i, 7)*k_alpha(i, 7)*m_alpha(i, l, 7) + & - cos45*(prod2 + prod4 - prod6 - prod8) - END DO - END DO - END IF - - END IF - - ! Sum over azimuths for case where SLOPE not equal to 1. - - IF (slope /= 1.) THEN - - ! Case with 4 azimuths. - - IF (naz == 4) THEN - DO l = lev1, lev2 - DO i = il1, il2 - flux_u(i, l) = ak_alpha(i, 1)*k_alpha(i, 1)* & - m_alpha(i, l, 1)**slope - ak_alpha(i, 3)*k_alpha(i, 3)*m_alpha(i, & - l, 3)**slope - flux_v(i, l) = ak_alpha(i, 2)*k_alpha(i, 2)* & - m_alpha(i, l, 2)**slope - ak_alpha(i, 4)*k_alpha(i, 4)*m_alpha(i, & - l, 4)**slope - END DO - END DO - END IF - - ! Case with 8 azimuths. - - IF (naz == 8) THEN - DO l = lev1, lev2 - DO i = il1, il2 - prod2 = ak_alpha(i, 2)*k_alpha(i, 2)*m_alpha(i, l, 2)**slope - prod4 = ak_alpha(i, 4)*k_alpha(i, 4)*m_alpha(i, l, 4)**slope - prod6 = ak_alpha(i, 6)*k_alpha(i, 6)*m_alpha(i, l, 6)**slope - prod8 = ak_alpha(i, 8)*k_alpha(i, 8)*m_alpha(i, l, 8)**slope - flux_u(i, l) = ak_alpha(i, 1)*k_alpha(i, 1)* & - m_alpha(i, l, 1)**slope - ak_alpha(i, 5)*k_alpha(i, 5)*m_alpha(i, & - l, 5)**slope + cos45*(prod2 - prod4 - prod6 + prod8) - flux_v(i, l) = ak_alpha(i, 3)*k_alpha(i, 3)* & - m_alpha(i, l, 3)**slope - ak_alpha(i, 7)*k_alpha(i, 7)*m_alpha(i, & - l, 7)**slope + cos45*(prod2 + prod4 - prod6 - prod8) - END DO - END DO - END IF - - END IF - - ! Calculate flux from sum. - - DO l = lev1, lev2 - DO i = il1, il2 - flux_u(i, l) = flux_u(i, l)*densb(i)/slope - flux_v(i, l) = flux_v(i, l)*densb(i)/slope - END DO - END DO - - ! Calculate drag at intermediate levels using centered differences - - DO l = lev1p, lev2m - DO i = il1, il2 - IF (lorms(i)) THEN - ! cc DENDZ2 = DENSITY(I,L) * ( ALT(I,L+1) - ALT(I,L-1) ) - dendz2 = density(i, l)*(alt(i, l - 1) - alt(i, l)) - ! cc DRAG_U(I,L) = - ( FLUX_U(I,L+1) - FLUX_U(I,L-1) ) / DENDZ2 - drag_u(i, l) = -(flux_u(i, l - 1) - flux_u(i, l))/dendz2 - ! cc DRAG_V(I,L) = - ( FLUX_V(I,L+1) - FLUX_V(I,L-1) ) / DENDZ2 - drag_v(i, l) = -(flux_v(i, l - 1) - flux_v(i, l))/dendz2 - - END IF - END DO - END DO - - ! Drag at first and last levels using one-side differences. - - DO i = il1, il2 - IF (lorms(i)) THEN - dendz = density(i, lev1)*(alt(i, lev1) - alt(i, lev1p)) - drag_u(i, lev1) = flux_u(i, lev1)/dendz - drag_v(i, lev1) = flux_v(i, lev1)/dendz - END IF - END DO - DO i = il1, il2 - IF (lorms(i)) THEN - dendz = density(i, lev2)*(alt(i, lev2m) - alt(i, lev2)) - drag_u(i, lev2) = -(flux_u(i, lev2m) - flux_u(i, lev2))/dendz - drag_v(i, lev2) = -(flux_v(i, lev2m) - flux_v(i, lev2))/dendz - END IF - END DO - IF (klevs > lev2) THEN - DO i = il1, il2 - IF (lorms(i)) THEN - dendz = density(i, lev2p)*(alt(i, lev2) - alt(i, lev2p)) - drag_u(i, lev2p) = -flux_u(i, lev2)/dendz - drag_v(i, lev2p) = -flux_v(i, lev2)/dendz - END IF - END DO - END IF - - RETURN - ! ----------------------------------------------------------------------- - END SUBROUTINE hines_flux - - SUBROUTINE hines_heat(heat, diffco, m_alpha, dmdz_alpha, ak_alpha, k_alpha, & - bvfreq, density, densb, sigma_t, visc_mol, kstar, slope, f2, f3, f5, f6, & - naz, il1, il2, lev1, lev2, klons, klevs, nazmth) - IMPLICIT NONE - ! This routine calculates the gravity wave induced heating and - ! diffusion coefficient on a longitude by altitude grid for - ! the Hines' Doppler spread gravity wave drag parameterization scheme. - - ! Aug. 6/95 - C. McLandress - - ! Output arguements: - - ! * HEAT = gravity wave heating (K/sec). - ! * DIFFCO = diffusion coefficient (m^2/sec) - - ! Input arguements: - - ! * M_ALPHA = cutoff vertical wavenumber (1/m). - ! * DMDZ_ALPHA = vertical derivative of cutoff wavenumber. - ! * AK_ALPHA = spectral amplitude factor of each azimuth - ! (i.e., {AjKj} in m^4/s^2). - ! * K_ALPHA = horizontal wavenumber of each azimuth (1/m). - ! * BVFREQ = background Brunt Vassala frequency (rad/sec). - ! * DENSITY = background density (kg/m^3). - ! * DENSB = background density at bottom level (kg/m^3). - ! * SIGMA_T = total rms horizontal wind (m/s). - ! * VISC_MOL = molecular viscosity (m^2/s). - ! * KSTAR = typical gravity wave horizontal wavenumber (1/m). - ! * SLOPE = slope of incident vertical wavenumber spectrum. - ! * F2,F3,F5,F6 = Hines's fudge factors. - ! * NAZ = actual number of horizontal azimuths used. - ! * IL1 = first longitudinal index to use (IL1 >= 1). - ! * IL2 = last longitudinal index to use (IL1 <= IL2 <= klonS). - ! * LEV1 = first altitude level to use (LEV1 >=1). - ! * LEV2 = last altitude level to use (LEV1 < LEV2 <= klevS). - ! * klonS = number of longitudes. - ! * klevS = number of vertical levels. - ! * NAZMTH = azimuthal array dimension (NAZMTH >= NAZ). - - INTEGER naz, il1, il2, lev1, lev2, klons, klevs, nazmth - REAL kstar(klons), slope, f2, f3, f5, f6 - REAL heat(klons, klevs), diffco(klons, klevs) - REAL m_alpha(klons, klevs, nazmth), dmdz_alpha(klons, klevs, nazmth) - REAL ak_alpha(klons, nazmth), k_alpha(klons, nazmth) - REAL bvfreq(klons, klevs), density(klons, klevs), densb(klons) - REAL sigma_t(klons, klevs), visc_mol(klons, klevs) - - ! Internal variables. - - INTEGER i, l, n - REAL m_sub_m_turb, m_sub_m_mol, m_sub_m, heatng - REAL visc, visc_min, cpgas, sm1 - - ! specific heat at constant pressure - - DATA cpgas/1004./ - - ! minimum permissible viscosity - - DATA visc_min/1.E-10/ - ! ----------------------------------------------------------------------- - - ! Initialize heating array. - - DO l = 1, klevs - DO i = 1, klons - heat(i, l) = 0. - END DO - END DO - - ! Perform sum over azimuths for case where SLOPE = 1. - - IF (slope == 1.) THEN - DO n = 1, naz - DO l = lev1, lev2 - DO i = il1, il2 - heat(i, l) = heat(i, l) + ak_alpha(i, n)*k_alpha(i, n)*dmdz_alpha(i & - , l, n) - END DO - END DO - END DO - END IF - - ! Perform sum over azimuths for case where SLOPE not 1. - - IF (slope /= 1.) THEN - sm1 = slope - 1. - DO n = 1, naz - DO l = lev1, lev2 - DO i = il1, il2 - heat(i, l) = heat(i, l) + ak_alpha(i, n)*k_alpha(i, n)*m_alpha(i, l & - , n)**sm1*dmdz_alpha(i, l, n) - END DO - END DO - END DO - END IF - - ! Heating and diffusion. - - DO l = lev1, lev2 - DO i = il1, il2 - - ! Maximum permissible value of cutoff wavenumber is the smaller - ! of the instability-induced wavenumber (M_SUB_M_TURB) and - ! that imposed by molecular viscosity (M_SUB_M_MOL). - - visc = amax1(visc_mol(i, l), visc_min) - m_sub_m_turb = bvfreq(i, l)/(f2*sigma_t(i, l)) - m_sub_m_mol = (bvfreq(i, l)*kstar(i)/visc)**0.33333333/f3 - m_sub_m = amin1(m_sub_m_turb, m_sub_m_mol) - - heatng = -heat(i, l)*f5*bvfreq(i, l)/m_sub_m*densb(i)/density(i, l) - diffco(i, l) = f6*heatng**0.33333333/m_sub_m**1.33333333 - heat(i, l) = heatng/cpgas - - END DO - END DO - - RETURN - ! ----------------------------------------------------------------------- - END SUBROUTINE hines_heat - - SUBROUTINE hines_sigma(sigma_t, sigma_alpha, sigsqh_alpha, naz, lev, il1, & - il2, klons, klevs, nazmth) - IMPLICIT NONE - ! This routine calculates the total rms and azimuthal rms horizontal - ! velocities at a given level on a longitude by altitude grid for - ! the Hines' Doppler spread GWD parameterization scheme. - ! NOTE: only four or eight azimuths can be used. - - ! Aug. 7/95 - C. McLandress - - ! Output arguements: - - ! * SIGMA_T = total rms horizontal wind (m/s). - ! * SIGMA_ALPHA = total rms wind in each azimuth (m/s). - - ! Input arguements: - - ! * SIGSQH_ALPHA = portion of wind variance from waves having wave - ! * normals in the alpha azimuth (m/s). - ! * NAZ = actual number of horizontal azimuths used (must be 4 or 8). - ! * LEV = altitude level to process. - ! * IL1 = first longitudinal index to use (IL1 >= 1). - ! * IL2 = last longitudinal index to use (IL1 <= IL2 <= klonS). - ! * klonS = number of longitudes. - ! * klevS = number of vertical levels. - ! * NAZMTH = azimuthal array dimension (NAZMTH >= NAZ). - - ! Subroutine arguements. - - INTEGER lev, naz, il1, il2 - INTEGER klons, klevs, nazmth - REAL sigma_t(klons, klevs) - REAL sigma_alpha(klons, klevs, nazmth) - REAL sigsqh_alpha(klons, klevs, nazmth) - - ! Internal variables. - - INTEGER i, n - REAL sum_even, sum_odd - ! ----------------------------------------------------------------------- - - ! Calculate azimuthal rms velocity for the 4 azimuth case. - - IF (naz == 4) THEN - DO i = il1, il2 - sigma_alpha(i, lev, 1) = sqrt(sigsqh_alpha(i, lev, 1) + sigsqh_alpha(i, lev, & - 3)) - sigma_alpha(i, lev, 2) = sqrt(sigsqh_alpha(i, lev, 2) + sigsqh_alpha(i, lev, & - 4)) - sigma_alpha(i, lev, 3) = sigma_alpha(i, lev, 1) - sigma_alpha(i, lev, 4) = sigma_alpha(i, lev, 2) - END DO - END IF - - ! Calculate azimuthal rms velocity for the 8 azimuth case. - - IF (naz == 8) THEN - DO i = il1, il2 - sum_odd = (sigsqh_alpha(i, lev, 1) + sigsqh_alpha(i, lev, 3) + & - sigsqh_alpha(i, lev, 5) + sigsqh_alpha(i, lev, 7))/2. - sum_even = (sigsqh_alpha(i, lev, 2) + sigsqh_alpha(i, lev, 4) + & - sigsqh_alpha(i, lev, 6) + sigsqh_alpha(i, lev, 8))/2. - sigma_alpha(i, lev, 1) = sqrt(sigsqh_alpha(i, lev, 1) + sigsqh_alpha(i, lev, & - 5) + sum_even) - sigma_alpha(i, lev, 2) = sqrt(sigsqh_alpha(i, lev, 2) + sigsqh_alpha(i, lev, & - 6) + sum_odd) - sigma_alpha(i, lev, 3) = sqrt(sigsqh_alpha(i, lev, 3) + sigsqh_alpha(i, lev, & - 7) + sum_even) - sigma_alpha(i, lev, 4) = sqrt(sigsqh_alpha(i, lev, 4) + sigsqh_alpha(i, lev, & - 8) + sum_odd) - sigma_alpha(i, lev, 5) = sigma_alpha(i, lev, 1) - sigma_alpha(i, lev, 6) = sigma_alpha(i, lev, 2) - sigma_alpha(i, lev, 7) = sigma_alpha(i, lev, 3) - sigma_alpha(i, lev, 8) = sigma_alpha(i, lev, 4) - END DO - END IF - - ! Calculate total rms velocity. - - DO i = il1, il2 - sigma_t(i, lev) = 0. - END DO - DO n = 1, naz - DO i = il1, il2 - sigma_t(i, lev) = sigma_t(i, lev) + sigsqh_alpha(i, lev, n) - END DO - END DO - DO i = il1, il2 - sigma_t(i, lev) = sqrt(sigma_t(i, lev)) - END DO - - RETURN - ! ----------------------------------------------------------------------- - END SUBROUTINE hines_sigma - - SUBROUTINE hines_intgrl(i_alpha, v_alpha, m_alpha, bvfb, slope, naz, lev, & - il1, il2, klons, klevs, nazmth, lorms) - IMPLICIT NONE - ! This routine calculates the vertical wavenumber integral - ! for a single vertical level at each azimuth on a longitude grid - ! for the Hines' Doppler spread GWD parameterization scheme. - ! NOTE: (1) only spectral slopes of 1, 1.5 or 2 are permitted. - ! (2) the integral is written in terms of the product QM - ! which by construction is always less than 1. Series - ! solutions are used for small |QM| and analytical solutions - ! for remaining values. - - ! Aug. 8/95 - C. McLandress - - ! Output arguement: - - ! * I_ALPHA = Hines' integral. - - ! Input arguements: - - ! * V_ALPHA = azimuthal wind component (m/s). - ! * M_ALPHA = azimuthal cutoff vertical wavenumber (1/m). - ! * BVFB = background Brunt Vassala frequency at model bottom. - ! * SLOPE = slope of initial vertical wavenumber spectrum - ! * (must use SLOPE = 1., 1.5 or 2.) - ! * NAZ = actual number of horizontal azimuths used. - ! * LEV = altitude level to process. - ! * IL1 = first longitudinal index to use (IL1 >= 1). - ! * IL2 = last longitudinal index to use (IL1 <= IL2 <= klonS). - ! * klonS = number of longitudes. - ! * klevS = number of vertical levels. - ! * NAZMTH = azimuthal array dimension (NAZMTH >= NAZ). - - ! * LORMS = .TRUE. for drag computation - - ! Constants in DATA statements: - - ! * QMIN = minimum value of Q_ALPHA (avoids indeterminant form of integral) - ! * QM_MIN = minimum value of Q_ALPHA * M_ALPHA (used to avoid numerical - ! * problems). - - INTEGER lev, naz, il1, il2, klons, klevs, nazmth - REAL i_alpha(klons, nazmth) - REAL v_alpha(klons, klevs, nazmth) - REAL m_alpha(klons, klevs, nazmth) - REAL bvfb(klons), slope - - LOGICAL lorms(klons) - - ! Internal variables. - - INTEGER i, n - REAL q_alpha, qm, sqrtqm, q_min, qm_min - - DATA q_min/1.0/, qm_min/0.01/ - ! ----------------------------------------------------------------------- - - ! For integer value SLOPE = 1. - - IF (slope == 1.) THEN - - DO n = 1, naz - DO i = il1, il2 - IF (lorms(i)) THEN - - q_alpha = v_alpha(i, lev, n)/bvfb(i) - qm = q_alpha*m_alpha(i, lev, n) - - ! If |QM| is small then use first 4 terms series of Taylor series - ! expansion of integral in order to avoid indeterminate form of - ! integral, - ! otherwise use analytical form of integral. - - IF (abs(q_alpha) < q_min .OR. abs(qm) < qm_min) THEN - IF (q_alpha == 0.) THEN - i_alpha(i, n) = m_alpha(i, lev, n)**2/2. - ELSE - i_alpha(i, n) = (qm**2/2.+qm**3/3.+qm**4/4.+qm**5/5.)/ & - q_alpha**2 - END IF - ELSE - i_alpha(i, n) = -(alog(1.-qm) + qm)/q_alpha**2 - END IF - - END IF - END DO - END DO - - END IF - - ! For integer value SLOPE = 2. - - IF (slope == 2.) THEN - - DO n = 1, naz - DO i = il1, il2 - IF (lorms(i)) THEN - - q_alpha = v_alpha(i, lev, n)/bvfb(i) - qm = q_alpha*m_alpha(i, lev, n) - - ! If |QM| is small then use first 4 terms series of Taylor series - ! expansion of integral in order to avoid indeterminate form of - ! integral, - ! otherwise use analytical form of integral. - - IF (abs(q_alpha) < q_min .OR. abs(qm) < qm_min) THEN - IF (q_alpha == 0.) THEN - i_alpha(i, n) = m_alpha(i, lev, n)**3/3. - ELSE - i_alpha(i, n) = (qm**3/3.+qm**4/4.+qm**5/5.+qm**6/6.)/ & - q_alpha**3 - END IF - ELSE - i_alpha(i, n) = -(alog(1.-qm) + qm + qm**2/2.)/q_alpha**3 - END IF - - END IF - END DO - END DO - - END IF - - ! For real value SLOPE = 1.5 - - IF (slope == 1.5) THEN - - DO n = 1, naz - DO i = il1, il2 - IF (lorms(i)) THEN - - q_alpha = v_alpha(i, lev, n)/bvfb(i) - qm = q_alpha*m_alpha(i, lev, n) - - ! If |QM| is small then use first 4 terms series of Taylor series - ! expansion of integral in order to avoid indeterminate form of - ! integral, - ! otherwise use analytical form of integral. - - IF (abs(q_alpha) < q_min .OR. abs(qm) < qm_min) THEN - IF (q_alpha == 0.) THEN - i_alpha(i, n) = m_alpha(i, lev, n)**2.5/2.5 - ELSE - i_alpha(i, n) = (qm/2.5 + qm**2/3.5 + qm**3/4.5 + qm**4/5.5)* & - m_alpha(i, lev, n)**1.5/q_alpha - END IF - ELSE - qm = abs(qm) - sqrtqm = sqrt(qm) - IF (q_alpha >= 0.) THEN - i_alpha(i, n) = (alog((1.+sqrtqm)/(1.-sqrtqm)) - 2.*sqrtqm*(1.+qm & - /3.))/q_alpha**2.5 - ELSE - i_alpha(i, n) = 2.*(atan(sqrtqm) + sqrtqm*(qm/3.-1.))/ & - abs(q_alpha)**2.5 - END IF - END IF - - END IF - END DO - END DO - - END IF - - ! If integral is negative (which in principal should not happen) then - ! print a message and some info since execution will abort when calculating - ! the variances. - - ! DO 80 N = 1,NAZ - ! DO 70 I = IL1,IL2 - ! IF (I_ALPHA(I,N).LT.0.) THEN - ! WRITE (6,*) - ! WRITE (6,*) '******************************' - ! WRITE (6,*) 'Hines integral I_ALPHA < 0 ' - ! WRITE (6,*) ' longitude I=',I - ! WRITE (6,*) ' azimuth N=',N - ! WRITE (6,*) ' level LEV=',LEV - ! WRITE (6,*) ' I_ALPHA =',I_ALPHA(I,N) - ! WRITE (6,*) ' V_ALPHA =',V_ALPHA(I,LEV,N) - ! WRITE (6,*) ' M_ALPHA =',M_ALPHA(I,LEV,N) - ! WRITE (6,*) ' Q_ALPHA =',V_ALPHA(I,LEV,N) / BVFB(I) - ! WRITE (6,*) ' QM =',V_ALPHA(I,LEV,N) / BVFB(I) - ! ^ * M_ALPHA(I,LEV,N) - ! WRITE (6,*) '******************************' - ! END IF - ! 70 CONTINUE - ! 80 CONTINUE - - RETURN - ! ----------------------------------------------------------------------- - END SUBROUTINE hines_intgrl - - SUBROUTINE hines_setup(naz, slope, f1, f2, f3, f5, f6, kstar, icutoff, & - alt_cutoff, smco, nsmax, iheatcal, k_alpha, ierror, nmessg, klons, & - nazmth, coslat) - IMPLICIT NONE - ! This routine specifies various parameters needed for the - ! the Hines' Doppler spread gravity wave drag parameterization scheme. - - ! Aug. 8/95 - C. McLandress - - ! Output arguements: - - ! * NAZ = actual number of horizontal azimuths used - ! * (code set up presently for only NAZ = 4 or 8). - ! * SLOPE = slope of incident vertical wavenumber spectrum - ! * (code set up presently for SLOPE 1., 1.5 or 2.). - ! * F1 = "fudge factor" used in calculation of trial value of - ! * azimuthal cutoff wavenumber M_ALPHA (1.2 <= F1 <= 1.9). - ! * F2 = "fudge factor" used in calculation of maximum - ! * permissible instabiliy-induced cutoff wavenumber - ! * M_SUB_M_TURB (0.1 <= F2 <= 1.4). - ! * F3 = "fudge factor" used in calculation of maximum - ! * permissible molecular viscosity-induced cutoff wavenumber - ! * M_SUB_M_MOL (0.1 <= F2 <= 1.4). - ! * F5 = "fudge factor" used in calculation of heating rate - ! * (1 <= F5 <= 3). - ! * F6 = "fudge factor" used in calculation of turbulent - ! * diffusivity coefficient. - ! * KSTAR = typical gravity wave horizontal wavenumber (1/m) - ! * used in calculation of M_SUB_M_TURB. - ! * ICUTOFF = 1 to exponentially damp off GWD, heating and diffusion - ! * arrays above ALT_CUTOFF; otherwise arrays not modified. - ! * ALT_CUTOFF = altitude in meters above which exponential decay applied. - ! * SMCO = smoother used to smooth cutoff vertical wavenumbers - ! * and total rms winds before calculating drag or heating. - ! * (==> a 1:SMCO:1 stencil used; SMCO >= 1.). - ! * NSMAX = number of times smoother applied ( >= 1), - ! * = 0 means no smoothing performed. - ! * IHEATCAL = 1 to calculate heating rates and diffusion coefficient. - ! * = 0 means only drag and flux calculated. - ! * K_ALPHA = horizontal wavenumber of each azimuth (1/m) which - ! * is set here to KSTAR. - ! * IERROR = error flag. - ! * = 0 no errors. - ! * = 10 ==> NAZ > NAZMTH - ! * = 20 ==> invalid number of azimuths (NAZ must be 4 or 8). - ! * = 30 ==> invalid slope (SLOPE must be 1., 1.5 or 2.). - ! * = 40 ==> invalid smoother (SMCO must be >= 1.) - - ! Input arguements: - - ! * NMESSG = output unit number where messages to be printed. - ! * klonS = number of longitudes. - ! * NAZMTH = azimuthal array dimension (NAZMTH >= NAZ). - - INTEGER naz, klons, nazmth, iheatcal, icutoff - INTEGER nmessg, nsmax, ierror - REAL kstar(klons), slope, f1, f2, f3, f5, f6, alt_cutoff, smco - REAL k_alpha(klons, nazmth), coslat(klons) - REAL ksmin, ksmax - - ! Internal variables. - - INTEGER i, n - ! ----------------------------------------------------------------------- - - ! Specify constants. - - naz = 8 - slope = 1. - f1 = 1.5 - f2 = 0.3 - f3 = 1.0 - f5 = 3.0 - f6 = 1.0 - ksmin = 1.E-5 - ksmax = 1.E-4 - DO i = 1, klons - kstar(i) = ksmin/(coslat(i) + (ksmin/ksmax)) - END DO - icutoff = 1 - alt_cutoff = 105.E3 - smco = 2.0 - ! SMCO = 1.0 - nsmax = 5 - ! NSMAX = 2 - iheatcal = 0 - - ! Print information to output file. - - ! WRITE (NMESSG,6000) - ! 6000 FORMAT (/' Subroutine HINES_SETUP:') - ! WRITE (NMESSG,*) ' SLOPE = ', SLOPE - ! WRITE (NMESSG,*) ' NAZ = ', NAZ - ! WRITE (NMESSG,*) ' F1,F2,F3 = ', F1, F2, F3 - ! WRITE (NMESSG,*) ' F5,F6 = ', F5, F6 - ! WRITE (NMESSG,*) ' KSTAR = ', KSTAR - ! > ,' COSLAT = ', COSLAT - ! IF (ICUTOFF .EQ. 1) THEN - ! WRITE (NMESSG,*) ' Drag exponentially damped above ', - ! & ALT_CUTOFF/1.E3 - ! END IF - ! IF (NSMAX.LT.1 ) THEN - ! WRITE (NMESSG,*) ' No smoothing of cutoff wavenumbers, etc' - ! ELSE - ! WRITE (NMESSG,*) ' Cutoff wavenumbers and sig_t smoothed:' - ! WRITE (NMESSG,*) ' SMCO =', SMCO - ! WRITE (NMESSG,*) ' NSMAX =', NSMAX - ! END IF - - ! Check that things are setup correctly and log error if not - - ierror = 0 - IF (naz > nazmth) ierror = 10 - IF (naz /= 4 .AND. naz /= 8) ierror = 20 - IF (slope /= 1. .AND. slope /= 1.5 .AND. slope /= 2.) ierror = 30 - IF (smco < 1.) ierror = 40 - - ! Use single value for azimuthal-dependent horizontal wavenumber. - - DO n = 1, naz - DO i = 1, klons - k_alpha(i, n) = kstar(i) - END DO - END DO - - RETURN - ! ----------------------------------------------------------------------- - END SUBROUTINE hines_setup - - SUBROUTINE hines_print(flux_u, flux_v, drag_u, drag_v, alt, sigma_t, & - sigma_alpha, v_alpha, m_alpha, iu_print, iv_print, nmessg, ilprt1, & - ilprt2, levprt1, levprt2, naz, klons, klevs, nazmth) - IMPLICIT NONE - ! Print out altitude profiles of various quantities from - ! Hines' Doppler spread gravity wave drag parameterization scheme. - ! (NOTE: only for NAZ = 4 or 8). - - ! Aug. 8/95 - C. McLandress - - ! Input arguements: - - ! * IU_PRINT = 1 to print out values in east-west direction. - ! * IV_PRINT = 1 to print out values in north-south direction. - ! * NMESSG = unit number for printed output. - ! * ILPRT1 = first longitudinal index to print. - ! * ILPRT2 = last longitudinal index to print. - ! * LEVPRT1 = first altitude level to print. - ! * LEVPRT2 = last altitude level to print. - - INTEGER naz, ilprt1, ilprt2, levprt1, levprt2 - INTEGER klons, klevs, nazmth - INTEGER iu_print, iv_print, nmessg - REAL flux_u(klons, klevs), flux_v(klons, klevs) - REAL drag_u(klons, klevs), drag_v(klons, klevs) - REAL alt(klons, klevs), sigma_t(klons, klevs) - REAL sigma_alpha(klons, klevs, nazmth) - REAL v_alpha(klons, klevs, nazmth), m_alpha(klons, klevs, nazmth) - - ! Internal variables. - - INTEGER n_east, n_west, n_north, n_south - INTEGER i, l - ! ----------------------------------------------------------------------- - - ! Azimuthal indices of cardinal directions. - - n_east = 1 - IF (naz == 4) THEN - n_west = 3 - n_north = 2 - n_south = 4 - ELSE IF (naz == 8) THEN - n_west = 5 - n_north = 3 - n_south = 7 - END IF - - ! Print out values for range of longitudes. - - DO i = ilprt1, ilprt2 - - ! Print east-west wind, sigmas, cutoff wavenumbers, flux and drag. - - IF (iu_print == 1) THEN - WRITE (nmessg, *) - WRITE (nmessg, 6001) i - WRITE (nmessg, 6005) -6001 FORMAT('Hines GW (east-west) at longitude I =', I3) -6005 FORMAT(15X, ' U ', 2X, 'sig_E', 2X, 'sig_T', 3X, 'm_E', 4X, 'm_W', 4X, & - 'fluxU', 5X, 'gwdU') - DO l = levprt1, levprt2 - WRITE (nmessg, 6701) alt(i, l)/1.E3, v_alpha(i, l, n_east), & - sigma_alpha(i, l, n_east), sigma_t(i, l), & - m_alpha(i, l, n_east)*1.E3, m_alpha(i, l, n_west)*1.E3, & - flux_u(i, l)*1.E5, drag_u(i, l)*24.*3600. - END DO -6701 FORMAT(' z=', F7.2, 1X, 3F7.1, 2F7.3, F9.4, F9.3) - END IF - - ! Print north-south winds, sigmas, cutoff wavenumbers, flux and drag. - - IF (iv_print == 1) THEN - WRITE (nmessg, *) - WRITE (nmessg, 6002) i -6002 FORMAT('Hines GW (north-south) at longitude I =', I3) - WRITE (nmessg, 6006) -6006 FORMAT(15X, ' V ', 2X, 'sig_N', 2X, 'sig_T', 3X, 'm_N', 4X, 'm_S', 4X, & - 'fluxV', 5X, 'gwdV') - DO l = levprt1, levprt2 - WRITE (nmessg, 6701) alt(i, l)/1.E3, v_alpha(i, l, n_north), & - sigma_alpha(i, l, n_north), sigma_t(i, l), & - m_alpha(i, l, n_north)*1.E3, m_alpha(i, l, n_south)*1.E3, & - flux_v(i, l)*1.E5, drag_v(i, l)*24.*3600. - END DO - END IF - - END DO - - RETURN - ! ----------------------------------------------------------------------- - END SUBROUTINE hines_print - - SUBROUTINE hines_exp(data, data_zmax, alt, alt_exp, iorder, il1, il2, lev1, & - lev2, klons, klevs) - IMPLICIT NONE - ! This routine exponentially damps a longitude by altitude array - ! of data above a specified altitude. - - ! Aug. 13/95 - C. McLandress - - ! Output arguements: - - ! * DATA = modified data array. - - ! Input arguements: - - ! * DATA = original data array. - ! * ALT = altitudes. - ! * ALT_EXP = altitude above which exponential decay applied. - ! * IORDER = 1 means vertical levels are indexed from top down - ! * (i.e., highest level indexed 1 and lowest level klevS); - ! * .NE. 1 highest level is index klevS. - ! * IL1 = first longitudinal index to use (IL1 >= 1). - ! * IL2 = last longitudinal index to use (IL1 <= IL2 <= klonS). - ! * LEV1 = first altitude level to use (LEV1 >=1). - ! * LEV2 = last altitude level to use (LEV1 < LEV2 <= klevS). - ! * klonS = number of longitudes. - ! * klevS = number of vertical - - ! Input work arrays: - - ! * DATA_ZMAX = data values just above altitude ALT_EXP. - - INTEGER iorder, il1, il2, lev1, lev2, klons, klevs - REAL alt_exp - REAL data(klons, klevs), data_zmax(klons), alt(klons, klevs) - - ! Internal variables. - - INTEGER levbot, levtop, lincr, i, l - REAL hscale - DATA hscale/5.E3/ - ! ----------------------------------------------------------------------- - - ! Index of lowest altitude level (bottom of drag calculation). - - levbot = lev2 - levtop = lev1 - lincr = 1 - IF (iorder /= 1) THEN - levbot = lev1 - levtop = lev2 - lincr = -1 - END IF - - ! Data values at first level above ALT_EXP. - - DO i = il1, il2 - DO l = levtop, levbot, lincr - IF (alt(i, l) >= alt_exp) THEN - data_zmax(i) = data(i, l) - END IF - END DO - END DO - - ! Exponentially damp field above ALT_EXP to model top at L=1. - - DO l = 1, lev2 - DO i = il1, il2 - IF (alt(i, l) >= alt_exp) THEN - data(i, l) = data_zmax(i)*exp((alt_exp - alt(i, l))/hscale) - END IF - END DO - END DO - - RETURN - ! ----------------------------------------------------------------------- - END SUBROUTINE hines_exp - - SUBROUTINE vert_smooth(data, work, coeff, nsmooth, il1, il2, lev1, lev2, & - klons, klevs) - IMPLICIT NONE - ! Smooth a longitude by altitude array in the vertical over a - ! specified number of levels using a three point smoother. - - ! NOTE: input array DATA is modified on output! - - ! Aug. 3/95 - C. McLandress - - ! Output arguement: - - ! * DATA = smoothed array (on output). - - ! Input arguements: - - ! * DATA = unsmoothed array of data (on input). - ! * WORK = work array of same dimension as DATA. - ! * COEFF = smoothing coefficient for a 1:COEFF:1 stencil. - ! * (e.g., COEFF = 2 will result in a smoother which - ! * weights the level L gridpoint by two and the two - ! * adjecent levels (L+1 and L-1) by one). - ! * NSMOOTH = number of times to smooth in vertical. - ! * (e.g., NSMOOTH=1 means smoothed only once, - ! * NSMOOTH=2 means smoothing repeated twice, etc.) - ! * IL1 = first longitudinal index to use (IL1 >= 1). - ! * IL2 = last longitudinal index to use (IL1 <= IL2 <= klonS). - ! * LEV1 = first altitude level to use (LEV1 >=1). - ! * LEV2 = last altitude level to use (LEV1 < LEV2 <= klevS). - ! * klonS = number of longitudes. - ! * klevS = number of vertical levels. - - ! Subroutine arguements. - - INTEGER nsmooth, il1, il2, lev1, lev2, klons, klevs - REAL coeff - REAL data(klons, klevs), work(klons, klevs) - - ! Internal variables. - - INTEGER i, l, ns, lev1p, lev2m - REAL sum_wts - ! ----------------------------------------------------------------------- - - ! Calculate sum of weights. - - sum_wts = coeff + 2. - - lev1p = lev1 + 1 - lev2m = lev2 - 1 - - ! Smooth NSMOOTH times - - DO ns = 1, nsmooth - - ! Copy data into work array. - - DO l = lev1, lev2 - DO i = il1, il2 - work(i, l) = data(i, l) - END DO - END DO - - ! Smooth array WORK in vertical direction and put into DATA. - - DO l = lev1p, lev2m - DO i = il1, il2 - data(i, l) = (work(i, l + 1) + coeff*work(i, l) + work(i, l - 1))/sum_wts - END DO - END DO - - END DO - - RETURN - ! ----------------------------------------------------------------------- - END SUBROUTINE vert_smooth - -END MODULE lmdz_gwd_hines - Index: LMDZ6/trunk/libf/phylmd/lmdz_gwd_ini.f90 =================================================================== --- LMDZ6/trunk/libf/phylmd/lmdz_gwd_ini.f90 (revision 6131) +++ LMDZ6/trunk/libf/phylmd/lmdz_gwd_ini.f90 (revision 6132) @@ -78,7 +78,4 @@ ! for non-orographic gravity wave drag - LOGICAL, SAVE, PROTECTED :: ok_hines = .false. ! activate non-orographic gravity wave drag from Hines - !$OMP THREADPRIVATE(ok_hines) - LOGICAL, SAVE, PROTECTED :: ok_gwd_rando = .false. ! activate non-orographic stochastic gravity wave drag params !$OMP THREADPRIVATE(ok_gwd_rando) @@ -177,5 +174,4 @@ CALL getin_p('ok_orodr', ok_orodr) CALL getin_p('ok_orolf', ok_orolf) - CALL getin_p('ok_hines', ok_hines) CALL getin_p('ok_gwd_rando', ok_gwd_rando) CALL getin_p('zpmm_orodr_t', zpmm_orodr_t) @@ -200,5 +196,4 @@ WRITE (lunout, *) 'gwd_ini, ok_orodr:', ok_orodr WRITE (lunout, *) 'gwd_ini, ok_orolf:', ok_orolf - WRITE (lunout, *) 'gwd_ini, ok_hines:', ok_hines WRITE (lunout, *) 'gwd_ini, ok_gwd_rando:', ok_gwd_rando WRITE (lunout, *) 'gwd_ini, zpmm_orodr_t:', zpmm_orodr_t @@ -546,9 +541,7 @@ END IF - IF (.NOT. ok_hines) THEN - IF (NW + 4*(NA - 1) + NA >= KLEV) THEN - abort_message = 'NW+3*NA>=KLEV Problem to generate waves associated with fronts' - CALL abort_physic(modname, abort_message, 1) - END IF + IF (NW + 4*(NA - 1) + NA >= KLEV) THEN + abort_message = 'NW+3*NA>=KLEV Problem to generate waves associated with fronts' + CALL abort_physic(modname, abort_message, 1) END IF Index: LMDZ6/trunk/libf/phylmd/phyetat0_mod.f90 =================================================================== --- LMDZ6/trunk/libf/phylmd/phyetat0_mod.f90 (revision 6131) +++ LMDZ6/trunk/libf/phylmd/phyetat0_mod.f90 (revision 6132) @@ -66,5 +66,5 @@ USE alpale_mod USE compbl_mod_h - USE lmdz_gwd_ini, ONLY: ok_hines, ok_gwd_rando, ok_orodr, ok_orolf + USE lmdz_gwd_ini, ONLY: ok_gwd_rando, ok_orodr, ok_orolf IMPLICIT NONE @@ -671,5 +671,5 @@ IF (ok_gwd_rando) found = & phyetat0_get(du_gwd_rando,"du_gwd_rando","du_gwd_rando",0.) - IF (.NOT. ok_hines .AND. ok_gwd_rando) found & + IF (ok_gwd_rando) found & = phyetat0_get(du_gwd_front,"du_gwd_front","du_gwd_front",0.) Index: LMDZ6/trunk/libf/phylmd/phys_local_var_mod.F90 =================================================================== --- LMDZ6/trunk/libf/phylmd/phys_local_var_mod.F90 (revision 6131) +++ LMDZ6/trunk/libf/phylmd/phys_local_var_mod.F90 (revision 6132) @@ -119,8 +119,4 @@ !$OMP THREADPRIVATE(d_u_lif, d_v_lif) ! Tendances Ondes de G non oro (runs strato). - REAL, SAVE, ALLOCATABLE :: du_gwd_hines(:,:) - !$OMP THREADPRIVATE(du_gwd_hines) - REAL, SAVE, ALLOCATABLE :: dv_gwd_hines(:,:) - !$OMP THREADPRIVATE(dv_gwd_hines) REAL, SAVE, ALLOCATABLE :: dv_gwd_rando(:,:) !$OMP THREADPRIVATE(dv_gwd_rando) @@ -131,6 +127,4 @@ REAL, SAVE, ALLOCATABLE :: west_gwstress(:,:) !$OMP THREADPRIVATE(west_gwstress) - REAL, SAVE, ALLOCATABLE :: dt_gwd_hines(:,:) - !$OMP THREADPRIVATE(dt_gwd_hines) ! tendance due a l'oxydation du methane REAL, SAVE, ALLOCATABLE :: d_q_ch4(:,:) @@ -981,10 +975,8 @@ ALLOCATE(topsw0_aero(klon,naero_grp), solsw0_aero(klon,naero_grp)) ALLOCATE(topswcf_aero(klon,3), solswcf_aero(klon,3)) - ALLOCATE(du_gwd_hines(klon,klev),dv_gwd_hines(klon,klev)) ALLOCATE(dv_gwd_rando(klon,klev),dv_gwd_front(klon,klev)) ALLOCATE(east_gwstress(klon,klev),west_gwstress(klon,klev)) east_gwstress(:,:)=0. !ym missing init west_gwstress(:,:)=0. !ym missing init - ALLOCATE(dt_gwd_hines(klon,klev)) ALLOCATE(d_q_ch4(klon,klev)) ALLOCATE(stratomask(klon,klev)) @@ -1503,5 +1495,4 @@ DEALLOCATE(load_tmp9) DEALLOCATE(load_tmp10) - DEALLOCATE(du_gwd_hines,dv_gwd_hines,dt_gwd_hines) DEALLOCATE(d_q_ch4) DEALLOCATE(dv_gwd_rando,dv_gwd_front) Index: LMDZ6/trunk/libf/phylmd/phys_output_ctrlout_mod.F90 =================================================================== --- LMDZ6/trunk/libf/phylmd/phys_output_ctrlout_mod.F90 (revision 6131) +++ LMDZ6/trunk/libf/phylmd/phys_output_ctrlout_mod.F90 (revision 6132) @@ -2402,8 +2402,4 @@ TYPE(ctrl_out), SAVE :: o_dvlif = ctrl_out((/ 4, 10, 10, 10, 10, 10, 11, 11, 11, 11/), & 'dvlif', 'Orography dV', 'm/s2', (/ ('', i=1, 10) /)) - TYPE(ctrl_out), SAVE :: o_du_gwd_hines = ctrl_out((/ 4, 10, 10, 10, 10, 10, 11, 11, 11, 11/), & - 'du_gwd_hines', 'Hines GWD dU', 'm/s2', (/ ('', i=1, 10) /)) - TYPE(ctrl_out), SAVE :: o_dv_gwd_hines = ctrl_out((/ 4, 10, 10, 10, 10, 10, 11, 11, 11, 11/), & - 'dv_gwd_hines', 'Hines GWD dV', 'm/s2', (/ ('', i=1, 10) /)) TYPE(ctrl_out), SAVE :: o_du_gwd_front = ctrl_out((/ 4, 10, 10, 10, 10, 10, 11, 11, 11, 11/), & 'du_gwd_front', 'Fronts GWD dU', 'm/s2', (/ ('', i=1, 10) /)) @@ -2418,6 +2414,4 @@ TYPE(ctrl_out), SAVE :: o_dtlif = ctrl_out((/ 4, 10, 10, 10, 10, 10, 11, 11, 11, 11/), & 'dtlif', 'Orography dT', 'K/s', (/ ('', i=1, 10) /)) - TYPE(ctrl_out), SAVE :: o_dthin = ctrl_out((/ 4, 10, 10, 10, 10, 10, 11, 11, 11, 11/), & - 'dthin', 'Hines GWD dT', 'K/s', (/ ('', i=1, 10) /)) TYPE(ctrl_out), SAVE :: o_dqch4 = ctrl_out((/ 4, 10, 10, 10, 10, 10, 11, 11, 11, 11/), & 'dqch4', 'H2O due to CH4 oxidation & photolysis', '(kg/kg)/s', (/ ('', i=1, 10) /)) @@ -2429,10 +2423,4 @@ = ctrl_out((/ 4, 10, 10, 10, 10, 10, 11, 11, 11, 11/), 'dv_gwd_rando', & "Random gravity waves dV/dt", "m/s2", (/ ('', i=1, 10) /)) - type(ctrl_out), save:: o_ustr_gwd_hines & - = ctrl_out((/ 4, 10, 10, 10, 10, 10, 11, 11, 11, 11/), 'ustr_gwd_hines', & - "zonal wind stress Hines gravity waves", "Pa", (/ ('', i=1, 10) /)) - type(ctrl_out), save:: o_vstr_gwd_hines & - = ctrl_out((/ 4, 10, 10, 10, 10, 10, 11, 11, 11, 11/), 'vstr_gwd_hines', & - "meridional wind stress Hines gravity waves", "Pa", (/ ('', i=1, 10) /)) type(ctrl_out), save:: o_ustr_gwd_front & = ctrl_out((/ 4, 10, 10, 10, 10, 10, 11, 11, 11, 11/), 'ustr_gwd_front', & Index: LMDZ6/trunk/libf/phylmd/phys_output_var_mod.f90 =================================================================== --- LMDZ6/trunk/libf/phylmd/phys_output_var_mod.f90 (revision 6131) +++ LMDZ6/trunk/libf/phylmd/phys_output_var_mod.f90 (revision 6132) @@ -142,8 +142,6 @@ LOGICAL, SAVE :: vars_defined = .FALSE. ! ug PAS THREADPRIVATE ET C'EST NORMAL - REAL, allocatable:: zustr_gwd_hines(:), zvstr_gwd_hines(:) ! (klon) REAL, allocatable:: zustr_gwd_front(:), zvstr_gwd_front(:) ! (klon) REAL, allocatable:: zustr_gwd_rando(:), zvstr_gwd_rando(:) ! (klon) - !$OMP THREADPRIVATE(zustr_gwd_hines, zvstr_gwd_hines) !$OMP THREADPRIVATE(zustr_gwd_front, zvstr_gwd_front) !$OMP THREADPRIVATE(zustr_gwd_rando, zvstr_gwd_rando) @@ -242,5 +240,4 @@ ! stress from non orographic gravity wave drag params - allocate(zustr_gwd_hines(klon), zvstr_gwd_hines(klon)) allocate(zustr_gwd_front(klon), zvstr_gwd_front(klon)) allocate(zustr_gwd_rando(klon), zvstr_gwd_rando(klon)) @@ -295,5 +292,4 @@ ! stress from non orographic gravity wave drag params - deallocate(zustr_gwd_hines, zvstr_gwd_hines) deallocate(zustr_gwd_front, zvstr_gwd_front) deallocate(zustr_gwd_rando, zvstr_gwd_rando) Index: LMDZ6/trunk/libf/phylmd/phys_output_write_mod.F90 =================================================================== --- LMDZ6/trunk/libf/phylmd/phys_output_write_mod.F90 (revision 6131) +++ LMDZ6/trunk/libf/phylmd/phys_output_write_mod.F90 (revision 6132) @@ -33,5 +33,5 @@ USE phystokenc_mod, ONLY: offline USE water_int_mod, ONLY: water_int - USE lmdz_gwd_ini, ONLY: ok_orodr, ok_orolf, ok_hines, ok_gwd_rando + USE lmdz_gwd_ini, ONLY: ok_orodr, ok_orolf, ok_gwd_rando USE conf_phys_m, ONLY : aerosol_couple, ok_ade, ok_aie, ok_volcan, & flag_aerosol, flag_aerosol_strat, ok_cdnc @@ -180,5 +180,5 @@ o_duvdf, o_dvvdf, o_duoro, o_dvoro, & o_dtoro, o_dulif, o_dvlif, o_dtlif, & - o_du_gwd_hines, o_dv_gwd_hines, o_dthin, o_dqch4, o_rsu, & + o_dqch4, o_rsu, & o_du_gwd_front, o_dv_gwd_front, & o_east_gwstress, o_west_gwstress, & @@ -202,5 +202,5 @@ o_dtr_sat, o_dtr_uscav, o_dtr_wet_cv, o_dtr_wet, & o_trac_cum, o_du_gwd_rando, o_dv_gwd_rando, & - o_ustr_gwd_hines,o_vstr_gwd_hines,o_ustr_gwd_rando,o_vstr_gwd_rando, & + o_ustr_gwd_rando,o_vstr_gwd_rando, & o_ustr_gwd_front,o_vstr_gwd_front, & o_sens_prec_liq_oce, o_sens_prec_liq_sic, & @@ -439,5 +439,5 @@ zw2, fraca, zmax_th, d_q_ajsb, d_t_ec, d_u_vdf, & d_v_vdf, d_u_oro, d_v_oro, d_t_oro, d_u_lif, & - d_v_lif, d_t_lif, du_gwd_hines, dv_gwd_hines, dt_gwd_hines, & + d_v_lif, d_t_lif, & dv_gwd_rando, dv_gwd_front, & east_gwstress, west_gwstress, & @@ -486,5 +486,5 @@ bils_ec,bils_ech, bils_tke, bils_kinetic, bils_latent, bils_enthalp, & itau_con, nfiles, clef_files, nid_files, dryaod_diag, & - zustr_gwd_hines, zvstr_gwd_hines,zustr_gwd_rando, zvstr_gwd_rando, & + zustr_gwd_rando, zvstr_gwd_rando, & zustr_gwd_front, zvstr_gwd_front, sza_o, & sens_prec_liq_o, sens_prec_sol_o, lat_prec_liq_o, lat_prec_sol_o, & @@ -2785,18 +2785,5 @@ ENDIF - IF (ok_hines) THEN - IF (vars_defined) zx_tmp_fi3d=du_gwd_hines/pdtphys - CALL histwrite_phy(o_du_gwd_hines, zx_tmp_fi3d) - - IF (vars_defined) zx_tmp_fi3d= dv_gwd_hines/pdtphys - CALL histwrite_phy(o_dv_gwd_hines, zx_tmp_fi3d) - - IF (vars_defined) zx_tmp_fi3d(1:klon,1:klev)=dt_gwd_hines(1:klon,1:klev)/pdtphys - CALL histwrite_phy(o_dthin, zx_tmp_fi3d) - CALL histwrite_phy(o_ustr_gwd_hines, zustr_gwd_hines) - CALL histwrite_phy(o_vstr_gwd_hines, zvstr_gwd_hines) - ENDIF - - IF (.not. ok_hines .and. ok_gwd_rando) THEN + IF (ok_gwd_rando) THEN IF (vars_defined) zx_tmp_fi3d=du_gwd_front / pdtphys CALL histwrite_phy(o_du_gwd_front, zx_tmp_fi3d) Index: LMDZ6/trunk/libf/phylmd/physiq_mod.F90 =================================================================== --- LMDZ6/trunk/libf/phylmd/physiq_mod.F90 (revision 6131) +++ LMDZ6/trunk/libf/phylmd/physiq_mod.F90 (revision 6132) @@ -254,5 +254,4 @@ d_t_ec, & ! - du_gwd_hines,dv_gwd_hines,dt_gwd_hines, & dv_gwd_rando,dv_gwd_front, & east_gwstress,west_gwstress, & @@ -4685,5 +4684,4 @@ d_u_oro, d_v_oro, d_t_oro, zustr_oro, zvstr_oro, & d_u_lif, d_v_lif, d_t_lif, zustr_lif, zvstr_lif, & - du_gwd_hines, dv_gwd_hines, dt_gwd_hines, zustr_gwd_hines, zvstr_gwd_hines, & du_gwd_front, dv_gwd_front, zustr_gwd_front, zvstr_gwd_front, & du_gwd_rando, dv_gwd_rando, zustr_gwd_rando, zvstr_gwd_rando, & Index: LMDZ6/trunk/libf/phylmdiso/lmdz_gwd_hines.f90 =================================================================== --- LMDZ6/trunk/libf/phylmdiso/lmdz_gwd_hines.f90 (revision 6131) +++ (revision ) @@ -1,1 +1,0 @@ -link ../phylmd/lmdz_gwd_hines.f90 Index: LMDZ6/trunk/libf/phylmdiso/phyetat0_mod.F90 =================================================================== --- LMDZ6/trunk/libf/phylmdiso/phyetat0_mod.F90 (revision 6131) +++ LMDZ6/trunk/libf/phylmdiso/phyetat0_mod.F90 (revision 6132) @@ -78,5 +78,5 @@ USE alpale_mod USE compbl_mod_h - USE lmdz_gwd_ini, ONLY: ok_hines, ok_gwd_rando, ok_orodr, ok_orolf + USE lmdz_gwd_ini, ONLY: ok_gwd_rando, ok_orodr, ok_orolf IMPLICIT NONE @@ -675,5 +675,5 @@ IF (ok_gwd_rando) found = & phyetat0_get(du_gwd_rando,"du_gwd_rando","du_gwd_rando",0.) - IF (.NOT. ok_hines .AND. ok_gwd_rando) found = & + IF (ok_gwd_rando) found = & phyetat0_get(du_gwd_front,"du_gwd_front","du_gwd_front",0.) Index: LMDZ6/trunk/libf/phylmdiso/phyredem.F90 =================================================================== --- LMDZ6/trunk/libf/phylmdiso/phyredem.F90 (revision 6131) +++ LMDZ6/trunk/libf/phylmdiso/phyredem.F90 (revision 6132) @@ -13,5 +13,5 @@ USE clesphys_mod_h USE dimphy, ONLY: klon, klev - USE lmdz_gwd_ini, ONLY: ok_gwd_rando, ok_hines, ok_orodr, ok_orolf + USE lmdz_gwd_ini, ONLY: ok_gwd_rando, ok_orodr, ok_orolf USE fonte_neige_mod, ONLY : fonte_neige_final USE pbl_surface_mod, ONLY : pbl_surface_final @@ -461,5 +461,5 @@ "tendency on zonal wind due to flott gravity waves", du_gwd_rando) - IF (.NOT. ok_hines .AND. ok_gwd_rando) CALL put_field(pass,"du_gwd_front", & + IF (ok_gwd_rando) CALL put_field(pass,"du_gwd_front", & "tendency on zonal wind due to acama gravity waves", du_gwd_front) Index: LMDZ6/trunk/libf/phylmdiso/physiq_mod.F90 =================================================================== --- LMDZ6/trunk/libf/phylmdiso/physiq_mod.F90 (revision 6131) +++ LMDZ6/trunk/libf/phylmdiso/physiq_mod.F90 (revision 6132) @@ -264,7 +264,4 @@ d_t_lif,d_u_lif,d_v_lif, & d_t_ec, & - ! - ! du_gwd_hines,dv_gwd_hines,d_t_hin, & FH a enlever avant 2026/05 d_t_hin pas dans phylmdiso - du_gwd_hines,dv_gwd_hines,dt_gwd_hines, & dv_gwd_rando,dv_gwd_front, & east_gwstress,west_gwstress, & @@ -6226,5 +6223,4 @@ d_u_oro, d_v_oro, d_t_oro, zustr_oro, zvstr_oro, & d_u_lif, d_v_lif, d_t_lif, zustr_lif, zvstr_lif, & - du_gwd_hines, dv_gwd_hines, dt_gwd_hines, zustr_gwd_hines, zvstr_gwd_hines, & du_gwd_front, dv_gwd_front, zustr_gwd_front, zvstr_gwd_front, & du_gwd_rando, dv_gwd_rando, zustr_gwd_rando, zvstr_gwd_rando, &