Index: LMDZ6/trunk/libf/phylmd/lmdz_call_gwd.f90 =================================================================== --- LMDZ6/trunk/libf/phylmd/lmdz_call_gwd.f90 (revision 6066) +++ LMDZ6/trunk/libf/phylmd/lmdz_call_gwd.f90 (revision 6068) @@ -20,5 +20,5 @@ !=================================================================== - SUBROUTINE call_gwd(klon, klev, nbsrf, is_ter, is_lic, abortphy, flag_inhib_tend, itap, JD_cur, JD_ref, JH_cur, & + SUBROUTINE call_gwd(klon, klev, nbsrf, is_ter, is_lic, is_ave, abortphy, flag_inhib_tend, itap, JD_cur, JD_ref, JH_cur, & pctsrf, is_sequential, phys_tstep, cell_area, longitude_deg, latitude_deg, pphis, & zstd, zpic, zmea, zval, zsig, zgam, zthe, pplay, paprs, presnivs, rain_fall, snow_fall, & @@ -61,5 +61,5 @@ REAL, INTENT(IN) :: JD_ref ! start day of the run REAL, INTENT(IN) :: JH_cur ! time of the day in seconds - INTEGER, INTENT(IN) :: is_ter, is_lic ! indices for land and landice subsurfaces + INTEGER, INTENT(IN) :: is_ter, is_lic, is_ave ! indices for land and landice subsurfaces and mesh-averaged LOGICAL, INTENT(IN) :: is_sequential ! sequential or parallel model REAL, INTENT(IN) :: phys_tstep ! time step [s] @@ -488,5 +488,5 @@ forall (k=1:klev) exner(:, k) = (pplay(:, k)/paprs(:, 1))**rkappa - CALL tend_to_tke(phys_tstep, klon, klev, nbsrf, paprs, exner, t_seri, u_seri, v_seri, dtadd, duadd, dvadd, pctsrf, tke) + CALL tend_to_tke(phys_tstep, klon, klev, nbsrf, is_ave, paprs, exner, t_seri, u_seri, v_seri, dtadd, duadd, dvadd, pctsrf, tke) ! Prevent pbl_tke_w from becoming negative Index: LMDZ6/trunk/libf/phylmd/lmdz_gwd_tendtotke.f90 =================================================================== --- LMDZ6/trunk/libf/phylmd/lmdz_gwd_tendtotke.f90 (revision 6066) +++ LMDZ6/trunk/libf/phylmd/lmdz_gwd_tendtotke.f90 (revision 6068) @@ -3,9 +3,9 @@ !******************* ! -! Subroutine that adds a tendency on the TKE created by the +! Subroutine that adds a tendency on the TKE created by the ! fluxes of momentum retrieved from the wind speed tendencies -! of the physics. Currently we only account for the contribution +! of the physics. Currently we only account for the contribution ! of the orographic gravity wave drag -! +! ! The basic concept is the following: ! the TKE equation writes de/dt = -u'w' du/dz -v'w' dv/dz +g/theta dtheta/dz +...... @@ -15,5 +15,5 @@ ! scheme, for instance: orographic gravity wave drag.... These contributions ! need to be accounted for. -! we explicitely calculate the fluxes, integrating the wind speed +! we explicitely calculate the fluxes, integrating the wind speed ! tendency from the top of the atmosphere ! @@ -33,118 +33,113 @@ !$gpum horizontal klon MODULE lmdz_gwd_tendtotke - PRIVATE + PRIVATE - PUBLIC tend_to_tke + PUBLIC tend_to_tke - CONTAINS +CONTAINS - SUBROUTINE tend_to_tke(dt,klon,klev,nbsrf,plev,exner,temp,windu,windv,dt_a,du_a,dv_a,pctsrf,tke) + SUBROUTINE tend_to_tke(dt, klon, klev, nbsrf, is_ave, plev, exner, temp, windu, windv, dt_a, du_a, dv_a, pctsrf, tke) -USE lmdz_gwd_ini, ONLY: RG, RCPD + USE lmdz_gwd_ini, ONLY: RG, RCPD -IMPLICIT NONE - + IMPLICIT NONE ! Declarations !============== - ! Inputs !------- - REAL, INTENT(IN) :: dt ! Time step [s] - INTEGER, INTENT(IN) :: klon,klev ! horizontal and vertical dimensions - INTEGER, INTENT(IN) :: nbsrf ! number of subsurfaces - REAL, DIMENSION(klon,klev+1), INTENT(IN) :: plev ! inter-layer pressure [Pa] - REAL, DIMENSION(klon,klev), INTENT(IN) :: temp ! temperature [K], grid-cell average or for a one subsurface - REAL, DIMENSION(klon,klev), INTENT(IN) :: windu ! zonal wind [m/s], grid-cell average or for a one subsurface - REAL, DIMENSION(klon,klev), INTENT(IN) :: windv ! meridonal wind [m/s], grid-cell average or for a one subsurface - REAL, DIMENSION(klon,klev), INTENT(IN) :: exner ! Fonction d'Exner = T/theta - REAL, DIMENSION(klon,klev), INTENT(IN) :: dt_a ! Temperature tendency [K], grid-cell average or for a one subsurface - REAL, DIMENSION(klon,klev), INTENT(IN) :: du_a ! Zonal wind speed tendency [m/s], grid-cell average or for a one subsurface - REAL, DIMENSION(klon,klev), INTENT(IN) :: dv_a ! Meridional wind speed tendency [m/s], grid-cell average or for a one subsurface - REAL, DIMENSION(klon,nbsrf), INTENT(IN) :: pctsrf ! fraction of each subsurface [0-1] + REAL, INTENT(IN) :: dt ! Time step [s] + INTEGER, INTENT(IN) :: klon, klev ! horizontal and vertical dimensions + INTEGER, INTENT(IN) :: nbsrf ! number of subsurfaces + INTEGER, INTENT(IN) :: is_ave ! index for mesh-averaged variables + REAL, DIMENSION(klon, klev + 1), INTENT(IN) :: plev ! inter-layer pressure [Pa] + REAL, DIMENSION(klon, klev), INTENT(IN) :: temp ! temperature [K], grid-cell average or for a one subsurface + REAL, DIMENSION(klon, klev), INTENT(IN) :: windu ! zonal wind [m/s], grid-cell average or for a one subsurface + REAL, DIMENSION(klon, klev), INTENT(IN) :: windv ! meridonal wind [m/s], grid-cell average or for a one subsurface + REAL, DIMENSION(klon, klev), INTENT(IN) :: exner ! Fonction d'Exner = T/theta + REAL, DIMENSION(klon, klev), INTENT(IN) :: dt_a ! Temperature tendency [K], grid-cell average or for a one subsurface + REAL, DIMENSION(klon, klev), INTENT(IN) :: du_a ! Zonal wind speed tendency [m/s], grid-cell average or for a one subsurface + REAL, DIMENSION(klon, klev), INTENT(IN) :: dv_a ! Meridional wind speed tendency [m/s], grid-cell average or for a one subsurface + REAL, DIMENSION(klon, nbsrf), INTENT(IN) :: pctsrf ! fraction of each subsurface [0-1] ! Inputs/Outputs !--------------- - REAL, DIMENSION(klon,klev+1,nbsrf+1), INTENT(INOUT) :: tke ! Turbulent Kinetic energy [m2/s2], grid-cell average or for a subsurface - + REAL, DIMENSION(klon, klev + 1, nbsrf + 1), INTENT(INOUT) :: tke ! Turbulent Kinetic energy [m2/s2], grid-cell average or for a subsurface ! Local !------- - - INTEGER i,k,isrf ! indices - REAL masse(klon,klev) ! mass in the layers [kg/m2] - REAL unsmasse(klon,klev+1) ! linear mass in the layers [kg/m2] - REAL flux_rhotw(klon,klev+1) ! flux massique de tempe. pot. rho*u'*theta' - REAL flux_rhouw(klon,klev+1) ! flux massique de quantit?? de mouvement rho*u'*w' [kg/m/s2] - REAL flux_rhovw(klon,klev+1) ! flux massique de quantit?? de mouvement rho*v'*w' [kg/m/s2] - REAL tendt(klon,klev) ! new temperature tke tendency [m2/s2/s] - REAL tendu(klon,klev) ! new zonal tke tendency [m2/s2/s] - REAL tendv(klon,klev) ! new meridonal tke tendency [m2/s2/s] - - - + INTEGER i, k, isrf ! indices + REAL masse(klon, klev) ! mass in the layers [kg/m2] + REAL unsmasse(klon, klev + 1) ! linear mass in the layers [kg/m2] + REAL flux_rhotw(klon, klev + 1) ! flux massique de tempe. pot. rho*u'*theta' + REAL flux_rhouw(klon, klev + 1) ! flux massique de quantit?? de mouvement rho*u'*w' [kg/m/s2] + REAL flux_rhovw(klon, klev + 1) ! flux massique de quantit?? de mouvement rho*v'*w' [kg/m/s2] + REAL tendt(klon, klev) ! new temperature tke tendency [m2/s2/s] + REAL tendu(klon, klev) ! new zonal tke tendency [m2/s2/s] + REAL tendv(klon, klev) ! new meridonal tke tendency [m2/s2/s] ! First calculations: !===================== - unsmasse(:,:)=0. - DO k=1,klev - masse(:,k)=(plev(:,k)-plev(:,k+1))/RG - unsmasse(:,k)=unsmasse(:,k)+0.5/masse(:,k) - unsmasse(:,k+1)=unsmasse(:,k+1)+0.5/masse(:,k) + unsmasse(:, :) = 0. + DO k = 1, klev + masse(:, k) = (plev(:, k) - plev(:, k + 1))/RG + unsmasse(:, k) = unsmasse(:, k) + 0.5/masse(:, k) + unsmasse(:, k + 1) = unsmasse(:, k + 1) + 0.5/masse(:, k) END DO - tendu(:,:)=0.0 - tendv(:,:)=0.0 + tendu(:, :) = 0.0 + tendv(:, :) = 0.0 ! Method 1: Calculation of fluxes using a downward integration !============================================================ - - ! Flux calculation - flux_rhotw(:,klev+1)=0. - flux_rhouw(:,klev+1)=0. - flux_rhovw(:,klev+1)=0. + flux_rhotw(:, klev + 1) = 0. + flux_rhouw(:, klev + 1) = 0. + flux_rhovw(:, klev + 1) = 0. - DO k=klev,1,-1 - flux_rhotw(:,k)=flux_rhotw(:,k+1)+masse(:,k)*dt_a(:,k)/exner(:,k) - flux_rhouw(:,k)=flux_rhouw(:,k+1)+masse(:,k)*du_a(:,k) - flux_rhovw(:,k)=flux_rhovw(:,k+1)+masse(:,k)*dv_a(:,k) - ENDDO - + DO k = klev, 1, -1 + flux_rhotw(:, k) = flux_rhotw(:, k + 1) + masse(:, k)*dt_a(:, k)/exner(:, k) + flux_rhouw(:, k) = flux_rhouw(:, k + 1) + masse(:, k)*du_a(:, k) + flux_rhovw(:, k) = flux_rhovw(:, k + 1) + masse(:, k)*dv_a(:, k) + END DO ! TKE update: - DO k=2,klev - tendt(:,k)=-flux_rhotw(:,k)*(exner(:,k)-exner(:,k-1))*unsmasse(:,k)*RCPD - tendu(:,k)=-flux_rhouw(:,k)*(windu(:,k)-windu(:,k-1))*unsmasse(:,k) - tendv(:,k)=-flux_rhovw(:,k)*(windv(:,k)-windv(:,k-1))*unsmasse(:,k) - ENDDO - tendt(:,1)=-flux_rhotw(:,1)*(exner(:,1)-1.)*unsmasse(:,1)*RCPD - tendu(:,1)=-1.*flux_rhouw(:,1)*windu(:,1)*unsmasse(:,1) - tendv(:,1)=-1.*flux_rhovw(:,1)*windv(:,1)*unsmasse(:,1) + DO k = 2, klev + tendt(:, k) = -flux_rhotw(:, k)*(exner(:, k) - exner(:, k - 1))*unsmasse(:, k)*RCPD + tendu(:, k) = -flux_rhouw(:, k)*(windu(:, k) - windu(:, k - 1))*unsmasse(:, k) + tendv(:, k) = -flux_rhovw(:, k)*(windv(:, k) - windv(:, k - 1))*unsmasse(:, k) + END DO + tendt(:, 1) = -flux_rhotw(:, 1)*(exner(:, 1) - 1.)*unsmasse(:, 1)*RCPD + tendu(:, 1) = -1.*flux_rhouw(:, 1)*windu(:, 1)*unsmasse(:, 1) + tendv(:, 1) = -1.*flux_rhovw(:, 1)*windv(:, 1)*unsmasse(:, 1) + DO isrf = 1, nbsrf + DO k = 1, klev + DO i = 1, klon + IF (pctsrf(i, isrf) > 0.) THEN + tke(i, k, isrf) = tke(i, k, isrf) + tendu(i, k) + tendv(i, k) + tendt(i, k) + tke(i, k, isrf) = max(tke(i, k, isrf), 1.e-10) + END IF + END DO + END DO + END DO - DO isrf=1,nbsrf - DO k=1,klev - DO i=1,klon - IF (pctsrf(i,isrf)>0.) THEN - tke(i,k,isrf)= tke(i,k,isrf)+tendu(i,k)+tendv(i,k)+tendt(i,k) - tke(i,k,isrf)= max(tke(i,k,isrf),1.e-10) - ENDIF - ENDDO - ENDDO - ENDDO +! Recompute the mesh-averaged TKE + tke(:, :, is_ave) = 0. + DO isrf = 1, nbsrf + DO k = 1, klev + DO i = 1, klon + tke(i, k, is_ave) = tke(i, k, is_ave) + pctsrf(i, isrf)*tke(i, k, isrf) + END DO + END DO + END DO -! The tendency on mean TKE (nbsrf+1) must be added as well -! or mesh-averaged TKE should be updated here - - - - END SUBROUTINE tend_to_tke + END SUBROUTINE tend_to_tke END MODULE lmdz_gwd_tendtotke Index: LMDZ6/trunk/libf/phylmd/physiq_mod.F90 =================================================================== --- LMDZ6/trunk/libf/phylmd/physiq_mod.F90 (revision 6066) +++ LMDZ6/trunk/libf/phylmd/physiq_mod.F90 (revision 6068) @@ -4724,5 +4724,5 @@ ! - CALL call_gwd(klon,klev,nbsrf,is_ter,is_lic,abortphy,flag_inhib_tend,itap, JD_cur, JD_ref, JH_cur, & + CALL call_gwd(klon,klev,nbsrf,is_ter,is_lic,is_ave,abortphy,flag_inhib_tend,itap, JD_cur, JD_ref, JH_cur, & pctsrf,is_sequential,phys_tstep,cell_area,longitude_deg,latitude_deg,pphis, & zstd,zpic,zmea,zval,zsig,zgam,zthe,pplay,paprs,presnivs,rain_fall,snow_fall, &