! $Id: top_bound.F90 5159 2024-08-02 19:58:25Z fairhead $ SUBROUTINE top_bound(vcov, ucov, teta, masse, dt) USE comconst_mod, ONLY: iflag_top_bound, mode_top_bound, & tau_top_bound USE comvert_mod, ONLY: presnivs, preff, scaleheight USE lmdz_iniprint, ONLY: lunout, prt_level USE lmdz_comdissipn, ONLY: tetaudiv, tetaurot, tetah, cdivu, crot, cdivh USE lmdz_comgeom2 USE lmdz_dimensions, ONLY: iim, jjm, llm, ndm USE lmdz_paramet IMPLICIT NONE ! ! .. DISSIPATION LINEAIRE A HAUT NIVEAU, RUN MESO, ! F. LOTT DEC. 2006 ! ( 10/12/06 ) !======================================================================= ! Auteur: F. LOTT ! ------- ! Objet: ! ------ ! Dissipation linéaire (ex top_bound de la physique) !======================================================================= ! top_bound sponge layer model: ! Quenching is modeled as: A(t)=Am+A0*exp(-lambda*t) ! where Am is the zonal average of the field (or zero), and lambda the inverse ! of the characteristic quenching/relaxation time scale ! Thus, assuming Am to be time-independent, field at time t+dt is given by: ! A(t+dt)=A(t)-(A(t)-Am)*(1-exp(-lambda*t)) ! Moreover lambda can be a function of model level (see below), and relaxation ! can be toward the average zonal field or just zero (see below). ! NB: top_bound sponge is only called from leapfrog if ok_strato=.TRUE. ! sponge parameters: (loaded/set in conf_gcm.F ; stored in comconst_mod) ! iflag_top_bound=0 for no sponge ! iflag_top_bound=1 for sponge over 4 topmost layers ! iflag_top_bound=2 for sponge from top to ~1% of top layer pressure ! mode_top_bound=0: no relaxation ! mode_top_bound=1: u and v relax towards 0 ! mode_top_bound=2: u and v relax towards their zonal mean ! mode_top_bound=3: u,v and pot. temp. relax towards their zonal mean ! tau_top_bound : inverse of charactericstic relaxation time scale at ! the topmost layer (Hz) ! Arguments: ! ---------- REAL, INTENT(INOUT) :: ucov(iip1, jjp1, llm) ! covariant zonal wind REAL, INTENT(INOUT) :: vcov(iip1, jjm, llm) ! covariant meridional wind REAL, INTENT(INOUT) :: teta(iip1, jjp1, llm) ! potential temperature REAL, INTENT(IN) :: masse(iip1, jjp1, llm) ! mass of atmosphere REAL, INTENT(IN) :: dt ! time step (s) of sponge model ! Local: ! ------ REAL :: massebx(iip1, jjp1, llm), masseby(iip1, jjm, llm), zm REAL :: uzon(jjp1, llm), vzon(jjm, llm), tzon(jjp1, llm) INTEGER :: i REAL, SAVE :: rdamp(llm) ! quenching coefficient REAL, save :: lambda(llm) ! inverse or quenching time scale (Hz) LOGICAL, SAVE :: first = .TRUE. INTEGER :: j, l IF (iflag_top_bound==0) return IF (first) THEN IF (iflag_top_bound==1) THEN ! sponge quenching over the topmost 4 atmospheric layers lambda(:) = 0. lambda(llm) = tau_top_bound lambda(llm - 1) = tau_top_bound / 2. lambda(llm - 2) = tau_top_bound / 4. lambda(llm - 3) = tau_top_bound / 8. ELSE IF (iflag_top_bound==2) THEN ! sponge quenching over topmost layers down to pressures which are ! higher than 100 times the topmost layer pressure lambda(:) = tau_top_bound & * max(presnivs(llm) / presnivs(:) - 0.01, 0.) endif ! quenching coefficient rdamp(:) ! rdamp(:)=dt*lambda(:) ! Explicit Euler approx. rdamp(:) = 1. - exp(-lambda(:) * dt) WRITE(lunout, *)'TOP_BOUND mode', mode_top_bound WRITE(lunout, *)'Sponge layer coefficients' WRITE(lunout, *)'p (Pa) z(km) tau(s) 1./tau (Hz)' DO l = 1, llm IF (rdamp(l)/=0.) THEN WRITE(lunout, '(6(1pe12.4,1x))') & presnivs(l), log(preff / presnivs(l)) * scaleheight, & 1. / lambda(l), lambda(l) endif enddo first = .FALSE. ENDIF ! of if (first) CALL massbar(masse, massebx, masseby) ! compute zonal average of vcov and u IF (mode_top_bound>=2) THEN DO l = 1, llm DO j = 1, jjm vzon(j, l) = 0. zm = 0. DO i = 1, iim ! NB: we can work using vcov zonal mean rather than v since the ! cv coefficient (which relates the two) only varies with latitudes vzon(j, l) = vzon(j, l) + vcov(i, j, l) * masseby(i, j, l) zm = zm + masseby(i, j, l) enddo vzon(j, l) = vzon(j, l) / zm enddo enddo DO l = 1, llm DO j = 2, jjm ! excluding poles uzon(j, l) = 0. zm = 0. DO i = 1, iim uzon(j, l) = uzon(j, l) + massebx(i, j, l) * ucov(i, j, l) / cu(i, j) zm = zm + massebx(i, j, l) enddo uzon(j, l) = uzon(j, l) / zm enddo enddo else ! ucov and vcov will relax towards 0 vzon(:, :) = 0. uzon(:, :) = 0. ENDIF ! of if (mode_top_bound.ge.2) ! compute zonal average of potential temperature, if necessary IF (mode_top_bound>=3) THEN DO l = 1, llm DO j = 2, jjm ! excluding poles zm = 0. tzon(j, l) = 0. DO i = 1, iim tzon(j, l) = tzon(j, l) + teta(i, j, l) * masse(i, j, l) zm = zm + masse(i, j, l) enddo tzon(j, l) = tzon(j, l) / zm enddo enddo ENDIF ! of if (mode_top_bound.ge.3) IF (mode_top_bound>=1) THEN ! Apply sponge quenching on vcov: DO l = 1, llm DO i = 1, iip1 DO j = 1, jjm vcov(i, j, l) = vcov(i, j, l) & - rdamp(l) * (vcov(i, j, l) - vzon(j, l)) enddo enddo enddo ! Apply sponge quenching on ucov: DO l = 1, llm DO i = 1, iip1 DO j = 2, jjm ! excluding poles ucov(i, j, l) = ucov(i, j, l) & - rdamp(l) * (ucov(i, j, l) - cu(i, j) * uzon(j, l)) enddo enddo enddo ENDIF ! of if (mode_top_bound.ge.1) IF (mode_top_bound>=3) THEN ! Apply sponge quenching on teta: DO l = 1, llm DO i = 1, iip1 DO j = 2, jjm ! excluding poles teta(i, j, l) = teta(i, j, l) & - rdamp(l) * (teta(i, j, l) - tzon(j, l)) enddo enddo enddo ENDIF ! of if (mode_top_bound.ge.3) END SUBROUTINE top_bound