- Timestamp:
- Jul 23, 2024, 7:14:34 PM (8 weeks ago)
- File:
-
- 1 moved
Legend:
- Unmodified
- Added
- Removed
-
LMDZ6/branches/Amaury_dev/libf/dyn3d_common/pentes_ini.f90
r5104 r5105 2 2 ! $Header$ 3 3 4 SUBROUTINE pentes_ini (q,w,masse,pbaru,pbarv,mode) 5 6 USE comconst_mod, ONLY: pi, dtvr 7 8 IMPLICIT NONE 9 10 c======================================================================= 11 c Adaptation LMDZ: A.Armengaud (LGGE) 12 c ---------------- 13 c 14 c ******************************************************************** 15 c Transport des traceurs par la methode des pentes 16 c ******************************************************************** 17 c Reference possible : Russel. G.L., Lerner J.A.: 18 c A new Finite-Differencing Scheme for Traceur Transport 19 c Equation , Journal of Applied Meteorology, pp 1483-1498,dec. 81 20 c ******************************************************************** 21 c q,w,masse,pbaru et pbarv 22 c sont des arguments d'entree pour le s-pg .... 23 c 24 c======================================================================= 25 26 27 include "dimensions.h" 28 include "paramet.h" 29 include "comgeom2.h" 30 31 c Arguments: 32 c ---------- 33 integer mode 34 REAL pbaru( ip1jmp1,llm ),pbarv( ip1jm,llm ) 35 REAL q( iip1,jjp1,llm,0:3) 36 REAL w( ip1jmp1,llm ) 37 REAL masse( iip1,jjp1,llm) 38 c Local: 39 c ------ 40 LOGICAL limit 41 REAL sm ( iip1,jjp1, llm ) 42 REAL s0( iip1,jjp1,llm ), sx( iip1,jjp1,llm ) 43 REAL sy( iip1,jjp1,llm ), sz( iip1,jjp1,llm ) 44 real masn,mass,zz 45 INTEGER i,j,l,iq 46 47 c modif Fred 24 03 96 48 49 real sinlon(iip1),sinlondlon(iip1) 50 real coslon(iip1),coslondlon(iip1) 51 save sinlon,coslon,sinlondlon,coslondlon 52 real dyn1,dyn2,dys1,dys2 53 real qpn,qps,dqzpn,dqzps 54 real smn,sms,s0n,s0s,sxn(iip1),sxs(iip1) 55 real qmin,zq,pente_max 56 c 57 REAL SSUM 58 integer ismax,ismin,lati,latf 59 EXTERNAL SSUM, ismin,ismax 60 logical first 61 save first 62 c fin modif 63 64 c EXTERNAL masskg 65 EXTERNAL advx 66 EXTERNAL advy 67 EXTERNAL advz 68 69 c modif Fred 24 03 96 70 data first/.TRUE./ 71 72 limit = .TRUE. 73 pente_max=2 74 c if (mode.eq.1.or.mode.eq.3) then 75 c if (mode.eq.1) then 76 if (mode>=1) then 77 lati=2 78 latf=jjm 79 else 80 lati=1 81 latf=jjp1 82 endif 83 84 qmin=0.4995 85 qmin=0. 86 if(first) then 87 PRINT*,'SCHEMA AMONT NOUVEAU' 88 first=.FALSE. 89 do i=2,iip1 90 coslon(i)=cos(rlonv(i)) 91 sinlon(i)=sin(rlonv(i)) 92 coslondlon(i)=coslon(i)*(rlonu(i)-rlonu(i-1))/pi 93 sinlondlon(i)=sinlon(i)*(rlonu(i)-rlonu(i-1))/pi 94 PRINT*,coslondlon(i),sinlondlon(i) 95 enddo 96 coslon(1)=coslon(iip1) 97 coslondlon(1)=coslondlon(iip1) 98 sinlon(1)=sinlon(iip1) 99 sinlondlon(1)=sinlondlon(iip1) 100 PRINT*,'sum sinlondlon ',ssum(iim,sinlondlon,1)/sinlondlon(1) 101 PRINT*,'sum coslondlon ',ssum(iim,coslondlon,1)/coslondlon(1) 102 DO l = 1,llm 103 DO j = 1,jjp1 104 DO i = 1,iip1 105 q ( i,j,l,1 )=0. 106 q ( i,j,l,2 )=0. 107 q ( i,j,l,3 )=0. 108 ENDDO 109 ENDDO 110 ENDDO 111 112 endif 113 c Fin modif Fred 114 115 c *** q contient les qqtes de traceur avant l'advection 116 117 c *** Affectation des tableaux S a partir de Q 118 c *** Rem : utilisation de SCOPY ulterieurement 119 120 DO l = 1,llm 121 DO j = 1,jjp1 122 DO i = 1,iip1 123 s0( i,j,llm+1-l ) = q ( i,j,l,0 ) 124 sx( i,j,llm+1-l ) = q ( i,j,l,1 ) 125 sy( i,j,llm+1-l ) = q ( i,j,l,2 ) 126 sz( i,j,llm+1-l ) = q ( i,j,l,3 ) 127 ENDDO 128 ENDDO 4 SUBROUTINE pentes_ini (q,w,masse,pbaru,pbarv,mode) 5 6 USE comconst_mod, ONLY: pi, dtvr 7 8 IMPLICIT NONE 9 10 !======================================================================= 11 ! Adaptation LMDZ: A.Armengaud (LGGE) 12 ! ---------------- 13 ! 14 ! ******************************************************************** 15 ! Transport des traceurs par la methode des pentes 16 ! ******************************************************************** 17 ! Reference possible : Russel. G.L., Lerner J.A.: 18 ! A new Finite-Differencing Scheme for Traceur Transport 19 ! Equation , Journal of Applied Meteorology, pp 1483-1498,dec. 81 20 ! ******************************************************************** 21 ! q,w,masse,pbaru et pbarv 22 ! sont des arguments d'entree pour le s-pg .... 23 ! 24 !======================================================================= 25 26 27 include "dimensions.h" 28 include "paramet.h" 29 include "comgeom2.h" 30 31 ! Arguments: 32 ! ---------- 33 integer :: mode 34 REAL :: pbaru( ip1jmp1,llm ),pbarv( ip1jm,llm ) 35 REAL :: q( iip1,jjp1,llm,0:3) 36 REAL :: w( ip1jmp1,llm ) 37 REAL :: masse( iip1,jjp1,llm) 38 ! Local: 39 ! ------ 40 LOGICAL :: limit 41 REAL :: sm ( iip1,jjp1, llm ) 42 REAL :: s0( iip1,jjp1,llm ), sx( iip1,jjp1,llm ) 43 REAL :: sy( iip1,jjp1,llm ), sz( iip1,jjp1,llm ) 44 real :: masn,mass,zz 45 INTEGER :: i,j,l,iq 46 47 ! modif Fred 24 03 96 48 49 real :: sinlon(iip1),sinlondlon(iip1) 50 real :: coslon(iip1),coslondlon(iip1) 51 save sinlon,coslon,sinlondlon,coslondlon 52 real :: dyn1,dyn2,dys1,dys2 53 real :: qpn,qps,dqzpn,dqzps 54 real :: smn,sms,s0n,s0s,sxn(iip1),sxs(iip1) 55 real :: qmin,zq,pente_max 56 ! 57 REAL :: SSUM 58 integer :: ismax,ismin,lati,latf 59 EXTERNAL SSUM, ismin,ismax 60 logical :: first 61 save first 62 ! fin modif 63 64 ! EXTERNAL masskg 65 EXTERNAL advx 66 EXTERNAL advy 67 EXTERNAL advz 68 69 ! modif Fred 24 03 96 70 data first/.TRUE./ 71 72 limit = .TRUE. 73 pente_max=2 74 ! if (mode.eq.1.or.mode.eq.3) then 75 ! if (mode.eq.1) then 76 if (mode>=1) then 77 lati=2 78 latf=jjm 79 else 80 lati=1 81 latf=jjp1 82 endif 83 84 qmin=0.4995 85 qmin=0. 86 if(first) then 87 PRINT*,'SCHEMA AMONT NOUVEAU' 88 first=.FALSE. 89 do i=2,iip1 90 coslon(i)=cos(rlonv(i)) 91 sinlon(i)=sin(rlonv(i)) 92 coslondlon(i)=coslon(i)*(rlonu(i)-rlonu(i-1))/pi 93 sinlondlon(i)=sinlon(i)*(rlonu(i)-rlonu(i-1))/pi 94 PRINT*,coslondlon(i),sinlondlon(i) 95 enddo 96 coslon(1)=coslon(iip1) 97 coslondlon(1)=coslondlon(iip1) 98 sinlon(1)=sinlon(iip1) 99 sinlondlon(1)=sinlondlon(iip1) 100 PRINT*,'sum sinlondlon ',ssum(iim,sinlondlon,1)/sinlondlon(1) 101 PRINT*,'sum coslondlon ',ssum(iim,coslondlon,1)/coslondlon(1) 102 DO l = 1,llm 103 DO j = 1,jjp1 104 DO i = 1,iip1 105 q ( i,j,l,1 )=0. 106 q ( i,j,l,2 )=0. 107 q ( i,j,l,3 )=0. 108 ENDDO 109 ENDDO 110 ENDDO 111 112 endif 113 ! Fin modif Fred 114 115 ! *** q contient les qqtes de traceur avant l'advection 116 117 ! *** Affectation des tableaux S a partir de Q 118 ! *** Rem : utilisation de SCOPY ulterieurement 119 120 DO l = 1,llm 121 DO j = 1,jjp1 122 DO i = 1,iip1 123 s0( i,j,llm+1-l ) = q ( i,j,l,0 ) 124 sx( i,j,llm+1-l ) = q ( i,j,l,1 ) 125 sy( i,j,llm+1-l ) = q ( i,j,l,2 ) 126 sz( i,j,llm+1-l ) = q ( i,j,l,3 ) 127 ENDDO 128 ENDDO 129 ENDDO 130 131 ! PRINT*,'----- S0 just before conversion -------' 132 ! PRINT*,'S0(16,12,1)=',s0(16,12,1) 133 ! PRINT*,'Q(16,12,1,4)=',q(16,12,1,4) 134 135 ! *** On calcule la masse d'air en kg 136 137 DO l = 1,llm 138 DO j = 1,jjp1 139 DO i = 1,iip1 140 sm ( i,j,llm+1-l)=masse( i,j,l ) 141 ENDDO 142 ENDDO 143 ENDDO 144 145 ! *** On converti les champs S en atome (resp. kg) 146 ! *** Les routines d'advection traitent les champs 147 ! *** a advecter si ces derniers sont en atome (resp. kg) 148 ! *** A optimiser !!! 149 150 DO l = 1,llm 151 DO j = 1,jjp1 152 DO i = 1,iip1 153 s0(i,j,l) = s0(i,j,l) * sm ( i,j,l ) 154 sx(i,j,l) = sx(i,j,l) * sm ( i,j,l ) 155 sy(i,j,l) = sy(i,j,l) * sm ( i,j,l ) 156 sz(i,j,l) = sz(i,j,l) * sm ( i,j,l ) 129 157 ENDDO 130 131 c PRINT*,'----- S0 just before conversion -------' 132 c PRINT*,'S0(16,12,1)=',s0(16,12,1) 133 c PRINT*,'Q(16,12,1,4)=',q(16,12,1,4) 134 135 c *** On calcule la masse d'air en kg 136 137 DO l = 1,llm 138 DO j = 1,jjp1 139 DO i = 1,iip1 140 sm ( i,j,llm+1-l)=masse( i,j,l ) 141 ENDDO 142 ENDDO 143 ENDDO 144 145 c *** On converti les champs S en atome (resp. kg) 146 c *** Les routines d'advection traitent les champs 147 c *** a advecter si ces derniers sont en atome (resp. kg) 148 c *** A optimiser !!! 149 150 DO l = 1,llm 151 DO j = 1,jjp1 152 DO i = 1,iip1 153 s0(i,j,l) = s0(i,j,l) * sm ( i,j,l ) 154 sx(i,j,l) = sx(i,j,l) * sm ( i,j,l ) 155 sy(i,j,l) = sy(i,j,l) * sm ( i,j,l ) 156 sz(i,j,l) = sz(i,j,l) * sm ( i,j,l ) 157 ENDDO 158 ENDDO 159 ENDDO 160 161 c ss0 = 0. 162 c DO l = 1,llm 163 c DO j = 1,jjp1 164 c DO i = 1,iim 165 c ss0 = ss0 + s0 ( i,j,l ) 166 c ENDDO 167 c ENDDO 168 c ENDDO 169 c PRINT*, 'valeur tot s0 avant advection=',ss0 170 171 c *** Appel des subroutines d'advection en X, en Y et en Z 172 c *** Advection avec "time-splitting" 173 174 c----------------------------------------------------------- 175 c PRINT*,'----- S0 just before ADVX -------' 176 c PRINT*,'S0(16,12,1)=',s0(16,12,1) 177 178 c----------------------------------------------------------- 179 c do l=1,llm 180 c do j=1,jjp1 181 c do i=1,iip1 182 c zq=s0(i,j,l)/sm(i,j,l) 183 c if(zq.lt.qmin) 184 c , PRINT*,'avant advx1, s0(',i,',',j,',',l,')=',zq 185 c enddo 186 c enddo 187 c enddo 188 CCC 189 if(mode==2) then 190 do l=1,llm 191 s0s=0. 192 s0n=0. 193 dyn1=0. 194 dys1=0. 195 dyn2=0. 196 dys2=0. 197 smn=0. 198 sms=0. 199 do i=1,iim 200 smn=smn+sm(i,1,l) 201 sms=sms+sm(i,jjp1,l) 202 s0n=s0n+s0(i,1,l) 203 s0s=s0s+s0(i,jjp1,l) 204 zz=sy(i,1,l)/sm(i,1,l) 205 dyn1=dyn1+sinlondlon(i)*zz 206 dyn2=dyn2+coslondlon(i)*zz 207 zz=sy(i,jjp1,l)/sm(i,jjp1,l) 208 dys1=dys1+sinlondlon(i)*zz 209 dys2=dys2+coslondlon(i)*zz 210 enddo 211 do i=1,iim 212 sy(i,1,l)=dyn1*sinlon(i)+dyn2*coslon(i) 213 sy(i,jjp1,l)=dys1*sinlon(i)+dys2*coslon(i) 214 enddo 215 do i=1,iim 216 s0(i,1,l)=s0n/smn+sy(i,1,l) 217 s0(i,jjp1,l)=s0s/sms-sy(i,jjp1,l) 218 enddo 219 220 s0(iip1,1,l)=s0(1,1,l) 221 s0(iip1,jjp1,l)=s0(1,jjp1,l) 222 223 do i=1,iim 224 sxn(i)=s0(i+1,1,l)-s0(i,1,l) 225 sxs(i)=s0(i+1,jjp1,l)-s0(i,jjp1,l) 226 c on rerentre les masses 227 enddo 228 do i=1,iim 229 sy(i,1,l)=sy(i,1,l)*sm(i,1,l) 230 sy(i,jjp1,l)=sy(i,jjp1,l)*sm(i,jjp1,l) 231 s0(i,1,l)=s0(i,1,l)*sm(i,1,l) 232 s0(i,jjp1,l)=s0(i,jjp1,l)*sm(i,jjp1,l) 233 enddo 234 sxn(iip1)=sxn(1) 235 sxs(iip1)=sxs(1) 236 do i=1,iim 237 sx(i+1,1,l)=0.25*(sxn(i)+sxn(i+1))*sm(i+1,1,l) 238 sx(i+1,jjp1,l)=0.25*(sxs(i)+sxs(i+1))*sm(i+1,jjp1,l) 239 enddo 240 s0(iip1,1,l)=s0(1,1,l) 241 s0(iip1,jjp1,l)=s0(1,jjp1,l) 242 sy(iip1,1,l)=sy(1,1,l) 243 sy(iip1,jjp1,l)=sy(1,jjp1,l) 244 sx(1,1,l)=sx(iip1,1,l) 245 sx(1,jjp1,l)=sx(iip1,jjp1,l) 246 enddo 247 endif 248 249 if (mode==4) then 250 do l=1,llm 251 do i=1,iip1 252 sx(i,1,l)=0. 253 sx(i,jjp1,l)=0. 254 sy(i,1,l)=0. 255 sy(i,jjp1,l)=0. 256 enddo 257 enddo 258 endif 259 CALL limx(s0,sx,sm,pente_max) 260 c CALL minmaxq(zq,1.e33,-1.e33,'avant advx ') 261 CALL advx( limit,.5*dtvr,pbaru,sm,s0,sx,sy,sz,lati,latf) 262 c CALL minmaxq(zq,1.e33,-1.e33,'avant advy ') 263 if (mode==4) then 264 do l=1,llm 265 do i=1,iip1 266 sx(i,1,l)=0. 267 sx(i,jjp1,l)=0. 268 sy(i,1,l)=0. 269 sy(i,jjp1,l)=0. 270 enddo 271 enddo 272 endif 273 CALL limy(s0,sy,sm,pente_max) 274 CALL advy( limit,.5*dtvr,pbarv,sm,s0,sx,sy,sz ) 275 c CALL minmaxq(zq,1.e33,-1.e33,'avant advz ') 276 do j=1,jjp1 277 do i=1,iip1 278 sz(i,j,1)=0. 279 sz(i,j,llm)=0. 280 enddo 158 ENDDO 159 ENDDO 160 161 ! ss0 = 0. 162 ! DO l = 1,llm 163 ! DO j = 1,jjp1 164 ! DO i = 1,iim 165 ! ss0 = ss0 + s0 ( i,j,l ) 166 ! ENDDO 167 ! ENDDO 168 ! ENDDO 169 ! PRINT*, 'valeur tot s0 avant advection=',ss0 170 171 ! *** Appel des subroutines d'advection en X, en Y et en Z 172 ! *** Advection avec "time-splitting" 173 174 !----------------------------------------------------------- 175 ! PRINT*,'----- S0 just before ADVX -------' 176 ! PRINT*,'S0(16,12,1)=',s0(16,12,1) 177 178 !----------------------------------------------------------- 179 ! do l=1,llm 180 ! do j=1,jjp1 181 ! do i=1,iip1 182 ! zq=s0(i,j,l)/sm(i,j,l) 183 ! if(zq.lt.qmin) 184 ! , PRINT*,'avant advx1, s0(',i,',',j,',',l,')=',zq 185 ! enddo 186 ! enddo 187 ! enddo 188 !CC 189 if(mode==2) then 190 do l=1,llm 191 s0s=0. 192 s0n=0. 193 dyn1=0. 194 dys1=0. 195 dyn2=0. 196 dys2=0. 197 smn=0. 198 sms=0. 199 do i=1,iim 200 smn=smn+sm(i,1,l) 201 sms=sms+sm(i,jjp1,l) 202 s0n=s0n+s0(i,1,l) 203 s0s=s0s+s0(i,jjp1,l) 204 zz=sy(i,1,l)/sm(i,1,l) 205 dyn1=dyn1+sinlondlon(i)*zz 206 dyn2=dyn2+coslondlon(i)*zz 207 zz=sy(i,jjp1,l)/sm(i,jjp1,l) 208 dys1=dys1+sinlondlon(i)*zz 209 dys2=dys2+coslondlon(i)*zz 210 enddo 211 do i=1,iim 212 sy(i,1,l)=dyn1*sinlon(i)+dyn2*coslon(i) 213 sy(i,jjp1,l)=dys1*sinlon(i)+dys2*coslon(i) 214 enddo 215 do i=1,iim 216 s0(i,1,l)=s0n/smn+sy(i,1,l) 217 s0(i,jjp1,l)=s0s/sms-sy(i,jjp1,l) 218 enddo 219 220 s0(iip1,1,l)=s0(1,1,l) 221 s0(iip1,jjp1,l)=s0(1,jjp1,l) 222 223 do i=1,iim 224 sxn(i)=s0(i+1,1,l)-s0(i,1,l) 225 sxs(i)=s0(i+1,jjp1,l)-s0(i,jjp1,l) 226 ! on rerentre les masses 227 enddo 228 do i=1,iim 229 sy(i,1,l)=sy(i,1,l)*sm(i,1,l) 230 sy(i,jjp1,l)=sy(i,jjp1,l)*sm(i,jjp1,l) 231 s0(i,1,l)=s0(i,1,l)*sm(i,1,l) 232 s0(i,jjp1,l)=s0(i,jjp1,l)*sm(i,jjp1,l) 233 enddo 234 sxn(iip1)=sxn(1) 235 sxs(iip1)=sxs(1) 236 do i=1,iim 237 sx(i+1,1,l)=0.25*(sxn(i)+sxn(i+1))*sm(i+1,1,l) 238 sx(i+1,jjp1,l)=0.25*(sxs(i)+sxs(i+1))*sm(i+1,jjp1,l) 239 enddo 240 s0(iip1,1,l)=s0(1,1,l) 241 s0(iip1,jjp1,l)=s0(1,jjp1,l) 242 sy(iip1,1,l)=sy(1,1,l) 243 sy(iip1,jjp1,l)=sy(1,jjp1,l) 244 sx(1,1,l)=sx(iip1,1,l) 245 sx(1,jjp1,l)=sx(iip1,jjp1,l) 246 enddo 247 endif 248 249 if (mode==4) then 250 do l=1,llm 251 do i=1,iip1 252 sx(i,1,l)=0. 253 sx(i,jjp1,l)=0. 254 sy(i,1,l)=0. 255 sy(i,jjp1,l)=0. 256 enddo 257 enddo 258 endif 259 CALL limx(s0,sx,sm,pente_max) 260 ! CALL minmaxq(zq,1.e33,-1.e33,'avant advx ') 261 CALL advx( limit,.5*dtvr,pbaru,sm,s0,sx,sy,sz,lati,latf) 262 ! CALL minmaxq(zq,1.e33,-1.e33,'avant advy ') 263 if (mode==4) then 264 do l=1,llm 265 do i=1,iip1 266 sx(i,1,l)=0. 267 sx(i,jjp1,l)=0. 268 sy(i,1,l)=0. 269 sy(i,jjp1,l)=0. 270 enddo 271 enddo 272 endif 273 CALL limy(s0,sy,sm,pente_max) 274 CALL advy( limit,.5*dtvr,pbarv,sm,s0,sx,sy,sz ) 275 ! CALL minmaxq(zq,1.e33,-1.e33,'avant advz ') 276 do j=1,jjp1 277 do i=1,iip1 278 sz(i,j,1)=0. 279 sz(i,j,llm)=0. 280 enddo 281 enddo 282 CALL limz(s0,sz,sm,pente_max) 283 CALL advz( limit,dtvr,w,sm,s0,sx,sy,sz ) 284 if (mode==4) then 285 do l=1,llm 286 do i=1,iip1 287 sx(i,1,l)=0. 288 sx(i,jjp1,l)=0. 289 sy(i,1,l)=0. 290 sy(i,jjp1,l)=0. 291 enddo 292 enddo 293 endif 294 CALL limy(s0,sy,sm,pente_max) 295 CALL advy( limit,.5*dtvr,pbarv,sm,s0,sx,sy,sz ) 296 do l=1,llm 297 do j=1,jjp1 298 sm(iip1,j,l)=sm(1,j,l) 299 s0(iip1,j,l)=s0(1,j,l) 300 sx(iip1,j,l)=sx(1,j,l) 301 sy(iip1,j,l)=sy(1,j,l) 302 sz(iip1,j,l)=sz(1,j,l) 303 enddo 304 enddo 305 306 307 ! CALL minmaxq(zq,1.e33,-1.e33,'avant advx ') 308 if (mode==4) then 309 do l=1,llm 310 do i=1,iip1 311 sx(i,1,l)=0. 312 sx(i,jjp1,l)=0. 313 sy(i,1,l)=0. 314 sy(i,jjp1,l)=0. 315 enddo 316 enddo 317 endif 318 CALL limx(s0,sx,sm,pente_max) 319 CALL advx( limit,.5*dtvr,pbaru,sm,s0,sx,sy,sz,lati,latf) 320 ! CALL minmaxq(zq,1.e33,-1.e33,'apres advx ') 321 ! do l=1,llm 322 ! do j=1,jjp1 323 ! do i=1,iip1 324 ! zq=s0(i,j,l)/sm(i,j,l) 325 ! if(zq.lt.qmin) 326 ! , PRINT*,'apres advx2, s0(',i,',',j,',',l,')=',zq 327 ! enddo 328 ! enddo 329 ! enddo 330 ! *** On repasse les S dans la variable q directement 14/10/94 331 ! On revient a des rapports de melange en divisant par la masse 332 333 ! En dehors des poles: 334 335 DO l = 1,llm 336 DO j = 1,jjp1 337 DO i = 1,iim 338 q(i,j,llm+1-l,0)=s0(i,j,l)/sm(i,j,l) 339 q(i,j,llm+1-l,1)=sx(i,j,l)/sm(i,j,l) 340 q(i,j,llm+1-l,2)=sy(i,j,l)/sm(i,j,l) 341 q(i,j,llm+1-l,3)=sz(i,j,l)/sm(i,j,l) 342 ENDDO 343 ENDDO 344 ENDDO 345 346 ! Traitements specifiques au pole 347 348 if(mode>=1) then 349 DO l=1,llm 350 ! filtrages aux poles 351 masn=ssum(iim,sm(1,1,l),1) 352 mass=ssum(iim,sm(1,jjp1,l),1) 353 qpn=ssum(iim,s0(1,1,l),1)/masn 354 qps=ssum(iim,s0(1,jjp1,l),1)/mass 355 dqzpn=ssum(iim,sz(1,1,l),1)/masn 356 dqzps=ssum(iim,sz(1,jjp1,l),1)/mass 357 do i=1,iip1 358 q( i,1,llm+1-l,3)=dqzpn 359 q( i,jjp1,llm+1-l,3)=dqzps 360 q( i,1,llm+1-l,0)=qpn 361 q( i,jjp1,llm+1-l,0)=qps 362 enddo 363 if(mode==3) then 364 dyn1=0. 365 dys1=0. 366 dyn2=0. 367 dys2=0. 368 do i=1,iim 369 dyn1=dyn1+sinlondlon(i)*sy(i,1,l)/sm(i,1,l) 370 dyn2=dyn2+coslondlon(i)*sy(i,1,l)/sm(i,1,l) 371 dys1=dys1+sinlondlon(i)*sy(i,jjp1,l)/sm(i,jjp1,l) 372 dys2=dys2+coslondlon(i)*sy(i,jjp1,l)/sm(i,jjp1,l) 373 enddo 374 do i=1,iim 375 q(i,1,llm+1-l,2)= & 376 (sinlon(i)*dyn1+coslon(i)*dyn2) 377 q(i,1,llm+1-l,0)=q(i,1,llm+1-l,0)+q(i,1,llm+1-l,2) 378 q(i,jjp1,llm+1-l,2)= & 379 (sinlon(i)*dys1+coslon(i)*dys2) 380 q(i,jjp1,llm+1-l,0)=q(i,jjp1,llm+1-l,0) & 381 -q(i,jjp1,llm+1-l,2) 382 enddo 383 endif 384 if(mode==1) then 385 ! on filtre les valeurs au bord de la "grande maille pole" 386 dyn1=0. 387 dys1=0. 388 dyn2=0. 389 dys2=0. 390 do i=1,iim 391 zz=s0(i,2,l)/sm(i,2,l)-q(i,1,llm+1-l,0) 392 dyn1=dyn1+sinlondlon(i)*zz 393 dyn2=dyn2+coslondlon(i)*zz 394 zz=q(i,jjp1,llm+1-l,0)-s0(i,jjm,l)/sm(i,jjm,l) 395 dys1=dys1+sinlondlon(i)*zz 396 dys2=dys2+coslondlon(i)*zz 397 enddo 398 do i=1,iim 399 q(i,1,llm+1-l,2)= & 400 (sinlon(i)*dyn1+coslon(i)*dyn2)/2. 401 q(i,1,llm+1-l,0)=q(i,1,llm+1-l,0)+q(i,1,llm+1-l,2) 402 q(i,jjp1,llm+1-l,2)= & 403 (sinlon(i)*dys1+coslon(i)*dys2)/2. 404 q(i,jjp1,llm+1-l,0)=q(i,jjp1,llm+1-l,0) & 405 -q(i,jjp1,llm+1-l,2) 406 enddo 407 q(iip1,1,llm+1-l,0)=q(1,1,llm+1-l,0) 408 q(iip1,jjp1,llm+1-l,0)=q(1,jjp1,llm+1-l,0) 409 410 do i=1,iim 411 sxn(i)=q(i+1,1,llm+1-l,0)-q(i,1,llm+1-l,0) 412 sxs(i)=q(i+1,jjp1,llm+1-l,0)-q(i,jjp1,llm+1-l,0) 413 enddo 414 sxn(iip1)=sxn(1) 415 sxs(iip1)=sxs(1) 416 do i=1,iim 417 q(i+1,1,llm+1-l,1)=0.25*(sxn(i)+sxn(i+1)) 418 q(i+1,jjp1,llm+1-l,1)=0.25*(sxs(i)+sxs(i+1)) 419 enddo 420 q(1,1,llm+1-l,1)=q(iip1,1,llm+1-l,1) 421 q(1,jjp1,llm+1-l,1)=q(iip1,jjp1,llm+1-l,1) 422 423 endif 424 425 ENDDO 426 endif 427 428 ! bouclage en longitude 429 do iq=0,3 430 do l=1,llm 431 do j=1,jjp1 432 q(iip1,j,l,iq)=q(1,j,l,iq) 433 enddo 434 enddo 435 enddo 436 437 ! PRINT*, ' SORTIE DE PENTES --- ca peut glisser ....' 438 439 DO l = 1,llm 440 DO j = 1,jjp1 441 DO i = 1,iip1 442 IF (q(i,j,l,0)<0.) THEN 443 ! PRINT*,'------------ BIP-----------' 444 ! PRINT*,'Q0(',i,j,l,')=',q(i,j,l,0) 445 ! PRINT*,'QX(',i,j,l,')=',q(i,j,l,1) 446 ! PRINT*,'QY(',i,j,l,')=',q(i,j,l,2) 447 ! PRINT*,'QZ(',i,j,l,')=',q(i,j,l,3) 448 ! PRINT*,' PBL EN SORTIE DE PENTES' 449 q(i,j,l,0)=0. 450 ! STOP 451 ENDIF 452 ENDDO 453 ENDDO 454 ENDDO 455 456 ! PRINT*, '-------------------------------------------' 457 458 do l=1,llm 459 do j=1,jjp1 460 do i=1,iip1 461 if(q(i,j,l,0)<qmin) & 462 PRINT*,'apres pentes, s0(',i,',',j,',',l,')=',q(i,j,l,0) 281 463 enddo 282 CALL limz(s0,sz,sm,pente_max)283 CALL advz( limit,dtvr,w,sm,s0,sx,sy,sz )284 if (mode==4) then285 do l=1,llm286 do i=1,iip1287 sx(i,1,l)=0.288 sx(i,jjp1,l)=0.289 sy(i,1,l)=0.290 sy(i,jjp1,l)=0.291 enddo292 enddo293 endif294 CALL limy(s0,sy,sm,pente_max)295 CALL advy( limit,.5*dtvr,pbarv,sm,s0,sx,sy,sz )296 do l=1,llm297 do j=1,jjp1298 sm(iip1,j,l)=sm(1,j,l)299 s0(iip1,j,l)=s0(1,j,l)300 sx(iip1,j,l)=sx(1,j,l)301 sy(iip1,j,l)=sy(1,j,l)302 sz(iip1,j,l)=sz(1,j,l)303 enddo304 enddo305 306 307 c CALL minmaxq(zq,1.e33,-1.e33,'avant advx ')308 if (mode==4) then309 do l=1,llm310 do i=1,iip1311 sx(i,1,l)=0.312 sx(i,jjp1,l)=0.313 sy(i,1,l)=0.314 sy(i,jjp1,l)=0.315 enddo316 enddo317 endif318 CALL limx(s0,sx,sm,pente_max)319 CALL advx( limit,.5*dtvr,pbaru,sm,s0,sx,sy,sz,lati,latf)320 c CALL minmaxq(zq,1.e33,-1.e33,'apres advx ')321 c do l=1,llm322 c do j=1,jjp1323 c do i=1,iip1324 c zq=s0(i,j,l)/sm(i,j,l)325 c if(zq.lt.qmin)326 c , PRINT*,'apres advx2, s0(',i,',',j,',',l,')=',zq327 c enddo328 c enddo329 c enddo330 c *** On repasse les S dans la variable q directement 14/10/94331 c On revient a des rapports de melange en divisant par la masse332 333 c En dehors des poles:334 335 DO l = 1,llm336 DO j = 1,jjp1337 DO i = 1,iim338 q(i,j,llm+1-l,0)=s0(i,j,l)/sm(i,j,l)339 q(i,j,llm+1-l,1)=sx(i,j,l)/sm(i,j,l)340 q(i,j,llm+1-l,2)=sy(i,j,l)/sm(i,j,l)341 q(i,j,llm+1-l,3)=sz(i,j,l)/sm(i,j,l)342 ENDDO343 ENDDO344 ENDDO345 346 c Traitements specifiques au pole347 348 if(mode>=1) then349 DO l=1,llm350 c filtrages aux poles351 masn=ssum(iim,sm(1,1,l),1)352 mass=ssum(iim,sm(1,jjp1,l),1)353 qpn=ssum(iim,s0(1,1,l),1)/masn354 qps=ssum(iim,s0(1,jjp1,l),1)/mass355 dqzpn=ssum(iim,sz(1,1,l),1)/masn356 dqzps=ssum(iim,sz(1,jjp1,l),1)/mass357 do i=1,iip1358 q( i,1,llm+1-l,3)=dqzpn359 q( i,jjp1,llm+1-l,3)=dqzps360 q( i,1,llm+1-l,0)=qpn361 q( i,jjp1,llm+1-l,0)=qps362 enddo363 if(mode==3) then364 dyn1=0.365 dys1=0.366 dyn2=0.367 dys2=0.368 do i=1,iim369 dyn1=dyn1+sinlondlon(i)*sy(i,1,l)/sm(i,1,l)370 dyn2=dyn2+coslondlon(i)*sy(i,1,l)/sm(i,1,l)371 dys1=dys1+sinlondlon(i)*sy(i,jjp1,l)/sm(i,jjp1,l)372 dys2=dys2+coslondlon(i)*sy(i,jjp1,l)/sm(i,jjp1,l)373 enddo374 do i=1,iim375 q(i,1,llm+1-l,2)=376 s (sinlon(i)*dyn1+coslon(i)*dyn2)377 q(i,1,llm+1-l,0)=q(i,1,llm+1-l,0)+q(i,1,llm+1-l,2)378 q(i,jjp1,llm+1-l,2)=379 s (sinlon(i)*dys1+coslon(i)*dys2)380 q(i,jjp1,llm+1-l,0)=q(i,jjp1,llm+1-l,0)381 s -q(i,jjp1,llm+1-l,2)382 enddo383 endif384 if(mode==1) then385 c on filtre les valeurs au bord de la "grande maille pole"386 dyn1=0.387 dys1=0.388 dyn2=0.389 dys2=0.390 do i=1,iim391 zz=s0(i,2,l)/sm(i,2,l)-q(i,1,llm+1-l,0)392 dyn1=dyn1+sinlondlon(i)*zz393 dyn2=dyn2+coslondlon(i)*zz394 zz=q(i,jjp1,llm+1-l,0)-s0(i,jjm,l)/sm(i,jjm,l)395 dys1=dys1+sinlondlon(i)*zz396 dys2=dys2+coslondlon(i)*zz397 enddo398 do i=1,iim399 q(i,1,llm+1-l,2)=400 s (sinlon(i)*dyn1+coslon(i)*dyn2)/2.401 q(i,1,llm+1-l,0)=q(i,1,llm+1-l,0)+q(i,1,llm+1-l,2)402 q(i,jjp1,llm+1-l,2)=403 s (sinlon(i)*dys1+coslon(i)*dys2)/2.404 q(i,jjp1,llm+1-l,0)=q(i,jjp1,llm+1-l,0)405 s -q(i,jjp1,llm+1-l,2)406 enddo407 q(iip1,1,llm+1-l,0)=q(1,1,llm+1-l,0)408 q(iip1,jjp1,llm+1-l,0)=q(1,jjp1,llm+1-l,0)409 410 do i=1,iim411 sxn(i)=q(i+1,1,llm+1-l,0)-q(i,1,llm+1-l,0)412 sxs(i)=q(i+1,jjp1,llm+1-l,0)-q(i,jjp1,llm+1-l,0)413 enddo414 sxn(iip1)=sxn(1)415 sxs(iip1)=sxs(1)416 do i=1,iim417 q(i+1,1,llm+1-l,1)=0.25*(sxn(i)+sxn(i+1))418 q(i+1,jjp1,llm+1-l,1)=0.25*(sxs(i)+sxs(i+1))419 enddo420 q(1,1,llm+1-l,1)=q(iip1,1,llm+1-l,1)421 q(1,jjp1,llm+1-l,1)=q(iip1,jjp1,llm+1-l,1)422 423 endif424 425 ENDDO426 endif427 428 c bouclage en longitude429 do iq=0,3430 do l=1,llm431 do j=1,jjp1432 q(iip1,j,l,iq)=q(1,j,l,iq)433 enddo434 enddo435 464 enddo 436 437 c PRINT*, ' SORTIE DE PENTES --- ca peut glisser ....' 438 439 DO l = 1,llm 440 DO j = 1,jjp1 441 DO i = 1,iip1 442 IF (q(i,j,l,0)<0.) THEN 443 c PRINT*,'------------ BIP-----------' 444 c PRINT*,'Q0(',i,j,l,')=',q(i,j,l,0) 445 c PRINT*,'QX(',i,j,l,')=',q(i,j,l,1) 446 c PRINT*,'QY(',i,j,l,')=',q(i,j,l,2) 447 c PRINT*,'QZ(',i,j,l,')=',q(i,j,l,3) 448 c PRINT*,' PBL EN SORTIE DE PENTES' 449 q(i,j,l,0)=0. 450 c STOP 451 ENDIF 452 ENDDO 453 ENDDO 454 ENDDO 455 456 c PRINT*, '-------------------------------------------' 457 458 do l=1,llm 459 do j=1,jjp1 460 do i=1,iip1 461 if(q(i,j,l,0)<qmin) 462 , PRINT*,'apres pentes, s0(',i,',',j,',',l,')=',q(i,j,l,0) 463 enddo 464 enddo 465 enddo 466 RETURN 467 END 468 469 470 471 472 473 474 475 476 477 478 479 465 enddo 466 RETURN 467 END SUBROUTINE pentes_ini 468 469 470 471 472 473 474 475 476 477 478 479
Note: See TracChangeset
for help on using the changeset viewer.
