Changeset 5246 for LMDZ6/trunk/libf/dyn3dmem/top_bound_loc.f90
- Timestamp:
- Oct 21, 2024, 2:58:45 PM (22 hours ago)
- File:
-
- 1 moved
-
LMDZ6/trunk/libf/dyn3dmem/top_bound_loc.f90 (moved) (moved from LMDZ6/trunk/libf/dyn3dmem/top_bound_loc.F) (1 diff)
Legend:
- Unmodified
- Added
- Removed
-
LMDZ6/trunk/libf/dyn3dmem/top_bound_loc.f90
r5245 r5246 2 2 ! $Id: $ 3 3 ! 4 5 6 USE comconst_mod, ONLY: iflag_top_bound, mode_top_bound,7 &tau_top_bound8 9 10 11 c 12 13 14 15 16 17 c.. DISSIPATION LINEAIRE A HAUT NIVEAU, RUN MESO,18 CF. LOTT DEC. 200619 c( 10/12/06 )20 21 c=======================================================================22 c 23 c Auteur: F. LOTT 24 c-------25 c 26 cObjet:27 c------28 c 29 cDissipation linéaire (ex top_bound de la physique)30 c 31 c=======================================================================32 33 ! top_bound sponge layer model:34 ! Quenching is modeled as: A(t)=Am+A0*exp(-lambda*t)35 ! where Am is the zonal average of the field (or zero), and lambda the inverse36 ! of the characteristic quenching/relaxation time scale37 ! Thus, assuming Am to be time-independent, field at time t+dt is given by:38 ! A(t+dt)=A(t)-(A(t)-Am)*(1-exp(-lambda*t))39 ! Moreover lambda can be a function of model level (see below), and relaxation40 ! can be toward the average zonal field or just zero (see below).41 42 ! NB: top_bound sponge is only called from leapfrog if ok_strato=.true.43 44 ! sponge parameters: (loaded/set in conf_gcm.F ; stored in comconst_mod)45 ! iflag_top_bound=0 for no sponge46 ! iflag_top_bound=1 for sponge over 4 topmost layers47 ! iflag_top_bound=2 for sponge from top to ~1% of top layer pressure48 ! mode_top_bound=0: no relaxation49 ! mode_top_bound=1: u and v relax towards 050 ! mode_top_bound=2: u and v relax towards their zonal mean51 ! mode_top_bound=3: u,v and pot. temp. relax towards their zonal mean52 ! tau_top_bound : inverse of charactericstic relaxation time scale at53 !the topmost layer (Hz)54 55 56 57 58 59 cArguments:60 c----------61 62 63 64 65 real,intent(in) :: masse(iip1,jjb_u:jje_u,llm) ! mass of atmosphere66 67 68 !REAL dv(iip1,jjb_v:jje_v,llm),du(iip1,jjb_u:jje_u,llm)69 !REAL dh(iip1,jjb_u:jje_u,llm)70 71 cLocal:72 c------73 REALmassebx(iip1,jjb_u:jje_u,llm),masseby(iip1,jjb_v:jje_v,llm)74 REALzm75 REALuzon(jjb_u:jje_u,llm),vzon(jjb_v:jje_v,llm)76 REALtzon(jjb_u:jje_u,llm)77 78 integeri79 80 81 82 INTEGERj,l,jjb,jje83 84 85 86 87 88 c$OMP BARRIER89 c$OMP MASTER90 91 ! sponge quenching over the topmost 4 atmospheric layers92 93 94 95 96 97 98 ! sponge quenching over topmost layers down to pressures which are99 ! higher than 100 times the topmost layer pressure100 lambda(:)=tau_top_bound101 s*max(presnivs(llm)/presnivs(:)-0.01,0.)102 103 104 ! quenching coefficient rdamp(:)105 !rdamp(:)=dt*lambda(:) ! Explicit Euler approx.106 107 108 109 110 111 112 113 write(lunout,'(6(1pe12.4,1x))')114 & presnivs(l),log(preff/presnivs(l))*scaleheight,115 &1./lambda(l),lambda(l)116 117 118 119 c$OMP END MASTER120 c$OMP BARRIER121 122 123 124 125 126 127 128 129 130 131 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK) 132 133 134 135 136 137 ! NB: we can work using vcov zonal mean rather than v since the138 ! cv coefficient (which relates the two) only varies with latitudes 139 140 141 142 143 144 145 c$OMP END DO NOWAIT 146 147 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK) 148 149 150 151 c$OMP END DO NOWAIT152 153 154 155 156 157 158 159 160 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK) 161 162 163 164 165 166 167 168 169 170 171 172 c$OMP END DO NOWAIT173 174 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK) 175 176 177 178 c$OMP END DO NOWAIT179 180 181 182 183 184 185 186 187 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK) 188 189 190 191 192 193 194 195 196 197 198 199 c$OMP END DO NOWAIT200 201 202 203 204 205 206 207 208 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK) 209 210 211 212 vcov(i,j,l)=vcov(i,j,l)213 &-rdamp(l)*(vcov(i,j,l)-vzon(j,l))214 215 216 217 c$OMP END DO NOWAIT218 219 220 221 222 223 224 225 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK) 226 227 228 229 ucov(i,j,l)=ucov(i,j,l)230 &-rdamp(l)*(ucov(i,j,l)-cu(i,j)*uzon(j,l))231 232 233 234 c$OMP END DO NOWAIT235 236 237 if (mode_top_bound.ge.3) then238 239 240 241 242 243 244 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK) 245 246 247 248 teta(i,j,l)=teta(i,j,l)249 &-rdamp(l)*(teta(i,j,l)-tzon(j,l))250 251 252 253 c$OMP END DO NOWAIT254 255 256 END 4 SUBROUTINE top_bound_loc(vcov,ucov,teta,masse,dt) 5 USE parallel_lmdz 6 USE comconst_mod, ONLY: iflag_top_bound, mode_top_bound, & 7 tau_top_bound 8 USE comvert_mod, ONLY: presnivs, preff, scaleheight 9 10 IMPLICIT NONE 11 ! 12 include "dimensions.h" 13 include "paramet.h" 14 include "comgeom2.h" 15 16 17 ! .. DISSIPATION LINEAIRE A HAUT NIVEAU, RUN MESO, 18 ! F. LOTT DEC. 2006 19 ! ( 10/12/06 ) 20 21 !======================================================================= 22 ! 23 ! Auteur: F. LOTT 24 ! ------- 25 ! 26 ! Objet: 27 ! ------ 28 ! 29 ! Dissipation linéaire (ex top_bound de la physique) 30 ! 31 !======================================================================= 32 33 ! top_bound sponge layer model: 34 ! Quenching is modeled as: A(t)=Am+A0*exp(-lambda*t) 35 ! where Am is the zonal average of the field (or zero), and lambda the inverse 36 ! of the characteristic quenching/relaxation time scale 37 ! Thus, assuming Am to be time-independent, field at time t+dt is given by: 38 ! A(t+dt)=A(t)-(A(t)-Am)*(1-exp(-lambda*t)) 39 ! Moreover lambda can be a function of model level (see below), and relaxation 40 ! can be toward the average zonal field or just zero (see below). 41 42 ! NB: top_bound sponge is only called from leapfrog if ok_strato=.true. 43 44 ! sponge parameters: (loaded/set in conf_gcm.F ; stored in comconst_mod) 45 ! iflag_top_bound=0 for no sponge 46 ! iflag_top_bound=1 for sponge over 4 topmost layers 47 ! iflag_top_bound=2 for sponge from top to ~1% of top layer pressure 48 ! mode_top_bound=0: no relaxation 49 ! mode_top_bound=1: u and v relax towards 0 50 ! mode_top_bound=2: u and v relax towards their zonal mean 51 ! mode_top_bound=3: u,v and pot. temp. relax towards their zonal mean 52 ! tau_top_bound : inverse of charactericstic relaxation time scale at 53 ! the topmost layer (Hz) 54 55 56 INCLUDE "comdissipn.h" 57 INCLUDE "iniprint.h" 58 59 ! Arguments: 60 ! ---------- 61 62 real,intent(inout) :: ucov(iip1,jjb_u:jje_u,llm) ! covariant zonal wind 63 real,intent(inout) :: vcov(iip1,jjb_v:jje_v,llm) ! covariant meridional wind 64 real,intent(inout) :: teta(iip1,jjb_u:jje_u,llm) ! potential temperature 65 real,intent(in) :: masse(iip1,jjb_u:jje_u,llm) ! mass of atmosphere 66 real,intent(in) :: dt ! time step (s) of sponge model 67 68 ! REAL dv(iip1,jjb_v:jje_v,llm),du(iip1,jjb_u:jje_u,llm) 69 ! REAL dh(iip1,jjb_u:jje_u,llm) 70 71 ! Local: 72 ! ------ 73 REAL :: massebx(iip1,jjb_u:jje_u,llm),masseby(iip1,jjb_v:jje_v,llm) 74 REAL :: zm 75 REAL :: uzon(jjb_u:jje_u,llm),vzon(jjb_v:jje_v,llm) 76 REAL :: tzon(jjb_u:jje_u,llm) 77 78 integer :: i 79 REAL,SAVE :: rdamp(llm) 80 real,save :: lambda(llm) ! inverse or quenching time scale (Hz) 81 LOGICAL,SAVE :: first=.true. 82 INTEGER :: j,l,jjb,jje 83 84 85 if (iflag_top_bound == 0) return 86 87 if (first) then 88 !$OMP BARRIER 89 !$OMP MASTER 90 if (iflag_top_bound == 1) then 91 ! sponge quenching over the topmost 4 atmospheric layers 92 lambda(:)=0. 93 lambda(llm)=tau_top_bound 94 lambda(llm-1)=tau_top_bound/2. 95 lambda(llm-2)=tau_top_bound/4. 96 lambda(llm-3)=tau_top_bound/8. 97 else if (iflag_top_bound == 2) then 98 ! sponge quenching over topmost layers down to pressures which are 99 ! higher than 100 times the topmost layer pressure 100 lambda(:)=tau_top_bound & 101 *max(presnivs(llm)/presnivs(:)-0.01,0.) 102 endif 103 104 ! quenching coefficient rdamp(:) 105 ! rdamp(:)=dt*lambda(:) ! Explicit Euler approx. 106 rdamp(:)=1.-exp(-lambda(:)*dt) 107 108 write(lunout,*)'TOP_BOUND mode',mode_top_bound 109 write(lunout,*)'Sponge layer coefficients' 110 write(lunout,*)'p (Pa) z(km) tau(s) 1./tau (Hz)' 111 do l=1,llm 112 if (rdamp(l).ne.0.) then 113 write(lunout,'(6(1pe12.4,1x))') & 114 presnivs(l),log(preff/presnivs(l))*scaleheight, & 115 1./lambda(l),lambda(l) 116 endif 117 enddo 118 first=.false. 119 !$OMP END MASTER 120 !$OMP BARRIER 121 endif ! of if (first) 122 123 124 CALL massbar_loc(masse,massebx,masseby) 125 126 ! ! compute zonal average of vcov (or set it to zero) 127 if (mode_top_bound.ge.2) then 128 jjb=jj_begin 129 jje=jj_end 130 IF (pole_sud) jje=jj_end-1 131 !$OMP DO SCHEDULE(STATIC,OMP_CHUNK) 132 do l=1,llm 133 do j=jjb,jje 134 zm=0. 135 vzon(j,l)=0 136 do i=1,iim 137 ! NB: we can work using vcov zonal mean rather than v since the 138 ! cv coefficient (which relates the two) only varies with latitudes 139 vzon(j,l)=vzon(j,l)+vcov(i,j,l)*masseby(i,j,l) 140 zm=zm+masseby(i,j,l) 141 enddo 142 vzon(j,l)=vzon(j,l)/zm 143 enddo 144 enddo 145 !$OMP END DO NOWAIT 146 else 147 !$OMP DO SCHEDULE(STATIC,OMP_CHUNK) 148 do l=1,llm 149 vzon(:,l)=0. 150 enddo 151 !$OMP END DO NOWAIT 152 endif ! of if (mode_top_bound.ge.2) 153 154 ! ! compute zonal average of u (or set it to zero) 155 if (mode_top_bound.ge.2) then 156 jjb=jj_begin 157 jje=jj_end 158 IF (pole_nord) jjb=jj_begin+1 159 IF (pole_sud) jje=jj_end-1 160 !$OMP DO SCHEDULE(STATIC,OMP_CHUNK) 161 do l=1,llm 162 do j=jjb,jje 163 uzon(j,l)=0. 164 zm=0. 165 do i=1,iim 166 uzon(j,l)=uzon(j,l)+massebx(i,j,l)*ucov(i,j,l)/cu(i,j) 167 zm=zm+massebx(i,j,l) 168 enddo 169 uzon(j,l)=uzon(j,l)/zm 170 enddo 171 enddo 172 !$OMP END DO NOWAIT 173 else 174 !$OMP DO SCHEDULE(STATIC,OMP_CHUNK) 175 do l=1,llm 176 uzon(:,l)=0. 177 enddo 178 !$OMP END DO NOWAIT 179 endif ! of if (mode_top_bound.ge.2) 180 181 ! ! compute zonal average of potential temperature, if necessary 182 if (mode_top_bound.ge.3) then 183 jjb=jj_begin 184 jje=jj_end 185 IF (pole_nord) jjb=jj_begin+1 186 IF (pole_sud) jje=jj_end-1 187 !$OMP DO SCHEDULE(STATIC,OMP_CHUNK) 188 do l=1,llm 189 do j=jjb,jje 190 zm=0. 191 tzon(j,l)=0. 192 do i=1,iim 193 tzon(j,l)=tzon(j,l)+teta(i,j,l)*masse(i,j,l) 194 zm=zm+masse(i,j,l) 195 enddo 196 tzon(j,l)=tzon(j,l)/zm 197 enddo 198 enddo 199 !$OMP END DO NOWAIT 200 endif ! of if (mode_top_bound.ge.3) 201 202 if (mode_top_bound.ge.1) then 203 ! ! Apply sponge quenching on vcov: 204 jjb=jj_begin 205 jje=jj_end 206 IF (pole_sud) jje=jj_end-1 207 208 !$OMP DO SCHEDULE(STATIC,OMP_CHUNK) 209 do l=1,llm 210 do j=jjb,jje 211 do i=1,iip1 212 vcov(i,j,l)=vcov(i,j,l) & 213 -rdamp(l)*(vcov(i,j,l)-vzon(j,l)) 214 enddo 215 enddo 216 enddo 217 !$OMP END DO NOWAIT 218 219 ! ! Apply sponge quenching on ucov: 220 jjb=jj_begin 221 jje=jj_end 222 IF (pole_nord) jjb=jj_begin+1 223 IF (pole_sud) jje=jj_end-1 224 225 !$OMP DO SCHEDULE(STATIC,OMP_CHUNK) 226 do l=1,llm 227 do j=jjb,jje 228 do i=1,iip1 229 ucov(i,j,l)=ucov(i,j,l) & 230 -rdamp(l)*(ucov(i,j,l)-cu(i,j)*uzon(j,l)) 231 enddo 232 enddo 233 enddo 234 !$OMP END DO NOWAIT 235 endif ! of if (mode_top_bound.ge.1) 236 237 if (mode_top_bound.ge.3) then 238 ! ! Apply sponge quenching on teta: 239 jjb=jj_begin 240 jje=jj_end 241 IF (pole_nord) jjb=jj_begin+1 242 IF (pole_sud) jje=jj_end-1 243 244 !$OMP DO SCHEDULE(STATIC,OMP_CHUNK) 245 do l=1,llm 246 do j=jjb,jje 247 do i=1,iip1 248 teta(i,j,l)=teta(i,j,l) & 249 -rdamp(l)*(teta(i,j,l)-tzon(j,l)) 250 enddo 251 enddo 252 enddo 253 !$OMP END DO NOWAIT 254 endif ! of if (mode_top_bond.ge.3) 255 256 END SUBROUTINE top_bound_loc
Note: See TracChangeset
for help on using the changeset viewer.
