! ! $Header$ ! SUBROUTINE advx(limit,dtx,pbaru,sm,s0, $ sx,sy,sz,lati,latf) IMPLICIT NONE CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C C first-order moments (FOM) advection of tracer in X direction C C C C Source : Pascal Simon (Meteo,CNRM) C C Adaptation : A.Armengaud (LGGE) juin 94 C C C C limit,dtx,pbaru,pbarv,sm,s0,sx,sy,sz C C sont des arguments d'entree pour le s-pg... C C C C sm,s0,sx,sy,sz C C sont les arguments de sortie pour le s-pg C C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C parametres principaux du modele C !----------------------------------------------------------------------- ! INCLUDE 'dimensions.h' ! ! dimensions.h contient les dimensions du modele ! ndm est tel que iim=2**ndm !----------------------------------------------------------------------- INTEGER iim,jjm,llm,ndm PARAMETER (iim= 128,jjm=96,llm=64,ndm=1) !----------------------------------------------------------------------- ! ! $Header$ ! ! ! ATTENTION!!!!: ce fichier include est compatible format fixe/format libre ! veillez n'utiliser que des ! pour les commentaires ! et bien positionner les & des lignes de continuation ! (les placer en colonne 6 et en colonne 73) ! ! !----------------------------------------------------------------------- ! INCLUDE 'paramet.h' INTEGER iip1,iip2,iip3,jjp1,llmp1,llmp2,llmm1 INTEGER kftd,ip1jm,ip1jmp1,ip1jmi1,ijp1llm INTEGER ijmllm,mvar INTEGER jcfil,jcfllm PARAMETER( iip1= iim+1,iip2=iim+2,iip3=iim+3 & & ,jjp1=jjm+1-1/jjm) PARAMETER( llmp1 = llm+1, llmp2 = llm+2, llmm1 = llm-1 ) PARAMETER( kftd = iim/2 -ndm ) PARAMETER( ip1jm = iip1*jjm, ip1jmp1= iip1*jjp1 ) PARAMETER( ip1jmi1= ip1jm - iip1 ) PARAMETER( ijp1llm= ip1jmp1 * llm, ijmllm= ip1jm * llm ) PARAMETER( mvar= ip1jmp1*( 2*llm+1) + ijmllm ) PARAMETER( jcfil=jjm/2+5, jcfllm=jcfil*llm ) !----------------------------------------------------------------------- ! ! $Id: comconst.h 1437 2010-09-30 08:29:10Z emillour $ ! !----------------------------------------------------------------------- ! INCLUDE comconst.h COMMON/comconsti/im,jm,lllm,imp1,jmp1,lllmm1,lllmp1,lcl, & & iflag_top_bound,mode_top_bound COMMON/comconstr/dtvr,daysec, & & pi,dtphys,dtdiss,rad,r,kappa,cotot,unsim,g,omeg & & ,dissip_fac_mid,dissip_fac_up,dissip_deltaz,dissip_hdelta & & ,dissip_pupstart ,tau_top_bound, & & daylen,molmass, ihf COMMON/cpdetvenus/cpp,nu_venus,t0_venus INTEGER im,jm,lllm,imp1,jmp1,lllmm1,lllmp1,lcl REAL dtvr ! dynamical time step (in s) REAL daysec !length (in s) of a standard day REAL pi ! something like 3.14159.... REAL dtphys ! (s) time step for the physics REAL dtdiss ! (s) time step for the dissipation REAL rad ! (m) radius of the planet REAL r ! Reduced Gas constant r=R/mu ! with R=8.31.. J.K-1.mol-1, mu: mol mass of atmosphere (kg/mol) REAL cpp ! Cp REAL kappa ! kappa=R/Cp REAL cotot REAL unsim ! = 1./iim REAL g ! (m/s2) gravity REAL omeg ! (rad/s) rotation rate of the planet ! Dissipation factors, for Earth model: REAL dissip_factz,dissip_zref !dissip_deltaz ! Dissipation factors, for other planets: REAL dissip_fac_mid,dissip_fac_up,dissip_deltaz,dissip_hdelta REAL dissip_pupstart INTEGER iflag_top_bound,mode_top_bound REAL tau_top_bound REAL daylen ! length of solar day, in 'standard' day length REAL molmass ! (g/mol) molar mass of the atmosphere REAL nu_venus,t0_venus ! coeffs needed for Cp(T), Venus atmosphere REAL ihf ! (W/m2) intrinsic heat flux for giant planets !----------------------------------------------------------------------- ! ! $Id: comvert.h 1654 2012-09-24 15:07:18Z aslmd $ ! !----------------------------------------------------------------------- ! INCLUDE 'comvert.h' COMMON/comvertr/ap(llm+1),bp(llm+1),presnivs(llm),dpres(llm), & & pa,preff,nivsigs(llm),nivsig(llm+1), & & aps(llm),bps(llm),scaleheight,pseudoalt(llm) common/comverti/disvert_type, pressure_exner real ap ! hybrid pressure contribution at interlayers real bp ! hybrid sigma contribution at interlayer real presnivs ! (reference) pressure at mid-layers real dpres real pa ! reference pressure (Pa) at which hybrid coordinates ! become purely pressure real preff ! reference surface pressure (Pa) real nivsigs real nivsig real aps ! hybrid pressure contribution at mid-layers real bps ! hybrid sigma contribution at mid-layers real scaleheight ! atmospheric (reference) scale height (km) real pseudoalt ! pseudo-altitude of model levels (km), based on presnivs(), ! preff and scaleheight integer disvert_type ! type of vertical discretization: ! 1: Earth (default for planet_type==earth), ! automatic generation ! 2: Planets (default for planet_type!=earth), ! using 'z2sig.def' (or 'esasig.def) file logical pressure_exner ! compute pressure inside layers using Exner function, else use mean ! of pressure values at interfaces !----------------------------------------------------------------------- C Arguments : C ----------- C dtx : frequence fictive d'appel du transport C pbaru, pbarv : flux de masse en x et y en Pa.m2.s-1 INTEGER ntra PARAMETER (ntra = 1) C ATTENTION partout ou on trouve ntra, insertion de boucle C possible dans l'avenir. REAL dtx REAL pbaru ( iip1,jjp1,llm ) C moments: SM total mass in each grid box C S0 mass of tracer in each grid box C Si 1rst order moment in i direction C REAL SM(iip1,jjp1,llm),S0(iip1,jjp1,llm,ntra) REAL sx(iip1,jjp1,llm,ntra) $ ,sy(iip1,jjp1,llm,ntra) REAL sz(iip1,jjp1,llm,ntra) C Local : C ------- C mass fluxes across the boundaries (UGRI,VGRI,WGRI) C mass fluxes in kg C declaration : REAL UGRI(iip1,jjp1,llm) C Rem : VGRI et WGRI ne sont pas utilises dans C cette subroutine ( advection en x uniquement ) C C Ti are the moments for the current latitude and level C REAL TM(iim) REAL T0(iim,ntra),TX(iim,ntra) REAL TY(iim,ntra),TZ(iim,ntra) REAL TEMPTM ! just a temporary variable C C the moments F are similarly defined and used as temporary C storage for portions of the grid boxes in transit C REAL FM(iim) REAL F0(iim,ntra),FX(iim,ntra) REAL FY(iim,ntra),FZ(iim,ntra) C C work arrays C REAL ALF(iim),ALF1(iim),ALFQ(iim),ALF1Q(iim) C REAL SMNEW(iim),UEXT(iim) C REAL sqi,sqf LOGICAL LIMIT INTEGER NUM(jjp1),LONK,NUMK INTEGER lon,lati,latf,niv INTEGER i,i2,i3,j,jv,l,k,itrac lon = iim niv = llm C *** Test de passage d'arguments ****** C ------------------------------------- DO 300 j = 1,jjp1 NUM(j) = 1 300 CONTINUE sqi = 0. sqf = 0. DO l = 1,llm DO j = 1,jjp1 DO i = 1,iim cIM 240305 sqi = sqi + S0(i,j,l,9) sqi = sqi + S0(i,j,l,ntra) ENDDO ENDDO ENDDO PRINT*,'-------- DIAG DANS ADVX - ENTREE ---------' PRINT*,'sqi=',sqi C Interface : adaptation nouveau modele C ------------------------------------- C C --------------------------------------------------------- C Conversion des flux de masses en kg/s C pbaru est en N/s d'ou : C ugri est en kg/s DO 500 l = 1,llm DO 500 j = 1,jjm+1 DO 500 i = 1,iip1 C ugri (i,j,llm+1-l) = pbaru (i,j,l) * ( dsig(l) / g ) ugri (i,j,llm+1-l) = pbaru (i,j,l) 500 CONTINUE C --------------------------------------------------------- C --------------------------------------------------------- C --------------------------------------------------------- C start here C C boucle principale sur les niveaux et les latitudes C DO 1 L=1,NIV DO 1 K=lati,latf C C initialisation C C program assumes periodic boundaries in X C DO 10 I=2,LON SMNEW(I)=SM(I,K,L)+(UGRI(I-1,K,L)-UGRI(I,K,L))*DTX 10 CONTINUE SMNEW(1)=SM(1,K,L)+(UGRI(LON,K,L)-UGRI(1,K,L))*DTX C C modifications for extended polar zones C NUMK=NUM(K) LONK=LON/NUMK C IF(NUMK.GT.1) THEN C DO 111 I=1,LON TM(I)=0. 111 CONTINUE DO 112 JV=1,NTRA DO 1120 I=1,LON T0(I,JV)=0. TX(I,JV)=0. TY(I,JV)=0. TZ(I,JV)=0. 1120 CONTINUE 112 CONTINUE C DO 11 I2=1,NUMK C DO 113 I=1,LONK I3=(I-1)*NUMK+I2 TM(I)=TM(I)+SM(I3,K,L) ALF(I)=SM(I3,K,L)/TM(I) ALF1(I)=1.-ALF(I) 113 CONTINUE C DO JV=1,NTRA DO I=1,LONK I3=(I-1)*NUMK+I2 TEMPTM=-ALF(I)*T0(I,JV)+ALF1(I) $ *S0(I3,K,L,JV) T0(I,JV)=T0(I,JV)+S0(I3,K,L,JV) TX(I,JV)=ALF(I) *sx(I3,K,L,JV)+ $ ALF1(I)*TX(I,JV) +3.*TEMPTM TY(I,JV)=TY(I,JV)+sy(I3,K,L,JV) TZ(I,JV)=TZ(I,JV)+sz(I3,K,L,JV) ENDDO ENDDO C 11 CONTINUE C ELSE C DO 115 I=1,LON TM(I)=SM(I,K,L) 115 CONTINUE DO 116 JV=1,NTRA DO 1160 I=1,LON T0(I,JV)=S0(I,K,L,JV) TX(I,JV)=sx(I,K,L,JV) TY(I,JV)=sy(I,K,L,JV) TZ(I,JV)=sz(I,K,L,JV) 1160 CONTINUE 116 CONTINUE C ENDIF C DO 117 I=1,LONK UEXT(I)=UGRI(I*NUMK,K,L) 117 CONTINUE C C place limits on appropriate moments before transport C (if flux-limiting is to be applied) C IF(.NOT.LIMIT) GO TO 13 C DO 12 JV=1,NTRA DO 120 I=1,LONK TX(I,JV)=SIGN(AMIN1(AMAX1(T0(I,JV),0.),ABS(TX(I,JV))),TX(I,JV)) 120 CONTINUE 12 CONTINUE C 13 CONTINUE C C calculate flux and moments between adjacent boxes C 1- create temporary moments/masses for partial boxes in transit C 2- reajusts moments remaining in the box C C flux from IP to I if U(I).lt.0 C DO 140 I=1,LONK-1 IF(UEXT(I).LT.0.) THEN FM(I)=-UEXT(I)*DTX ALF(I)=FM(I)/TM(I+1) TM(I+1)=TM(I+1)-FM(I) ENDIF 140 CONTINUE C I=LONK IF(UEXT(I).LT.0.) THEN FM(I)=-UEXT(I)*DTX ALF(I)=FM(I)/TM(1) TM(1)=TM(1)-FM(I) ENDIF C C flux from I to IP if U(I).gt.0 C DO 141 I=1,LONK IF(UEXT(I).GE.0.) THEN FM(I)=UEXT(I)*DTX ALF(I)=FM(I)/TM(I) TM(I)=TM(I)-FM(I) ENDIF 141 CONTINUE C DO 142 I=1,LONK ALFQ(I)=ALF(I)*ALF(I) ALF1(I)=1.-ALF(I) ALF1Q(I)=ALF1(I)*ALF1(I) 142 CONTINUE C DO 150 JV=1,NTRA DO 1500 I=1,LONK-1 C IF(UEXT(I).LT.0.) THEN C F0(I,JV)=ALF (I)* ( T0(I+1,JV)-ALF1(I)*TX(I+1,JV) ) FX(I,JV)=ALFQ(I)*TX(I+1,JV) FY(I,JV)=ALF (I)*TY(I+1,JV) FZ(I,JV)=ALF (I)*TZ(I+1,JV) C T0(I+1,JV)=T0(I+1,JV)-F0(I,JV) TX(I+1,JV)=ALF1Q(I)*TX(I+1,JV) TY(I+1,JV)=TY(I+1,JV)-FY(I,JV) TZ(I+1,JV)=TZ(I+1,JV)-FZ(I,JV) C ENDIF C 1500 CONTINUE 150 CONTINUE C I=LONK IF(UEXT(I).LT.0.) THEN C DO 151 JV=1,NTRA C F0 (I,JV)=ALF (I)* ( T0(1,JV)-ALF1(I)*TX(1,JV) ) FX (I,JV)=ALFQ(I)*TX(1,JV) FY (I,JV)=ALF (I)*TY(1,JV) FZ (I,JV)=ALF (I)*TZ(1,JV) C T0(1,JV)=T0(1,JV)-F0(I,JV) TX(1,JV)=ALF1Q(I)*TX(1,JV) TY(1,JV)=TY(1,JV)-FY(I,JV) TZ(1,JV)=TZ(1,JV)-FZ(I,JV) C 151 CONTINUE C ENDIF C DO 152 JV=1,NTRA DO 1520 I=1,LONK C IF(UEXT(I).GE.0.) THEN C F0(I,JV)=ALF (I)* ( T0(I,JV)+ALF1(I)*TX(I,JV) ) FX(I,JV)=ALFQ(I)*TX(I,JV) FY(I,JV)=ALF (I)*TY(I,JV) FZ(I,JV)=ALF (I)*TZ(I,JV) C T0(I,JV)=T0(I,JV)-F0(I,JV) TX(I,JV)=ALF1Q(I)*TX(I,JV) TY(I,JV)=TY(I,JV)-FY(I,JV) TZ(I,JV)=TZ(I,JV)-FZ(I,JV) C ENDIF C 1520 CONTINUE 152 CONTINUE C C puts the temporary moments Fi into appropriate neighboring boxes C DO 160 I=1,LONK IF(UEXT(I).LT.0.) THEN TM(I)=TM(I)+FM(I) ALF(I)=FM(I)/TM(I) ENDIF 160 CONTINUE C DO 161 I=1,LONK-1 IF(UEXT(I).GE.0.) THEN TM(I+1)=TM(I+1)+FM(I) ALF(I)=FM(I)/TM(I+1) ENDIF 161 CONTINUE C I=LONK IF(UEXT(I).GE.0.) THEN TM(1)=TM(1)+FM(I) ALF(I)=FM(I)/TM(1) ENDIF C DO 162 I=1,LONK ALF1(I)=1.-ALF(I) 162 CONTINUE C DO 170 JV=1,NTRA DO 1700 I=1,LONK C IF(UEXT(I).LT.0.) THEN C TEMPTM=-ALF(I)*T0(I,JV)+ALF1(I)*F0(I,JV) T0(I,JV)=T0(I,JV)+F0(I,JV) TX(I,JV)=ALF(I)*FX(I,JV)+ALF1(I)*TX(I,JV)+3.*TEMPTM TY(I,JV)=TY(I,JV)+FY(I,JV) TZ(I,JV)=TZ(I,JV)+FZ(I,JV) C ENDIF C 1700 CONTINUE 170 CONTINUE C DO 171 JV=1,NTRA DO 1710 I=1,LONK-1 C IF(UEXT(I).GE.0.) THEN C TEMPTM=ALF(I)*T0(I+1,JV)-ALF1(I)*F0(I,JV) T0(I+1,JV)=T0(I+1,JV)+F0(I,JV) TX(I+1,JV)=ALF(I)*FX(I,JV)+ALF1(I)*TX(I+1,JV)+3.*TEMPTM TY(I+1,JV)=TY(I+1,JV)+FY(I,JV) TZ(I+1,JV)=TZ(I+1,JV)+FZ(I,JV) C ENDIF C 1710 CONTINUE 171 CONTINUE C I=LONK IF(UEXT(I).GE.0.) THEN DO 172 JV=1,NTRA TEMPTM=ALF(I)*T0(1,JV)-ALF1(I)*F0(I,JV) T0(1,JV)=T0(1,JV)+F0(I,JV) TX(1,JV)=ALF(I)*FX(I,JV)+ALF1(I)*TX(1,JV)+3.*TEMPTM TY(1,JV)=TY(1,JV)+FY(I,JV) TZ(1,JV)=TZ(1,JV)+FZ(I,JV) 172 CONTINUE ENDIF C C retour aux mailles d'origine (passage des Tij aux Sij) C IF(NUMK.GT.1) THEN C DO 180 I2=1,NUMK C DO 180 I=1,LONK C I3=I2+(I-1)*NUMK SM(I3,K,L)=SMNEW(I3) ALF(I)=SMNEW(I3)/TM(I) TM(I)=TM(I)-SMNEW(I3) C ALFQ(I)=ALF(I)*ALF(I) ALF1(I)=1.-ALF(I) ALF1Q(I)=ALF1(I)*ALF1(I) C 180 CONTINUE C DO JV=1,NTRA DO I=1,LONK C I3=I2+(I-1)*NUMK S0(I3,K,L,JV)=ALF (I) $ * (T0(I,JV)-ALF1(I)*TX(I,JV)) sx(I3,K,L,JV)=ALFQ(I)*TX(I,JV) sy(I3,K,L,JV)=ALF (I)*TY(I,JV) sz(I3,K,L,JV)=ALF (I)*TZ(I,JV) C C reajusts moments remaining in the box C T0(I,JV)=T0(I,JV)-S0(I3,K,L,JV) TX(I,JV)=ALF1Q(I)*TX(I,JV) TY(I,JV)=TY(I,JV)-sy(I3,K,L,JV) TZ(I,JV)=TZ(I,JV)-sz(I3,K,L,JV) ENDDO ENDDO C C ELSE C DO 190 I=1,LON SM(I,K,L)=TM(I) 190 CONTINUE DO 191 JV=1,NTRA DO 1910 I=1,LON S0(I,K,L,JV)=T0(I,JV) sx(I,K,L,JV)=TX(I,JV) sy(I,K,L,JV)=TY(I,JV) sz(I,K,L,JV)=TZ(I,JV) 1910 CONTINUE 191 CONTINUE C ENDIF C 1 CONTINUE C C ----------- AA Test en fin de ADVX ------ Controle des S* c OK c DO 9998 l = 1, llm c DO 9998 j = 1, jjp1 c DO 9998 i = 1, iip1 c IF (S0(i,j,l,ntra).lt.0..and.LIMIT) THEN c PRINT*, '-------------------' c PRINT*, 'En fin de ADVX' c PRINT*,'SM(',i,j,l,')=',SM(i,j,l) c PRINT*,'S0(',i,j,l,')=',S0(i,j,l,ntra) c print*, 'sx(',i,j,l,')=',sx(i,j,l,ntra) c print*, 'sy(',i,j,l,')=',sy(i,j,l,ntra) c print*, 'sz(',i,j,l,')=',sz(i,j,l,ntra) c WRITE (*,*) 'On arrete !! - pbl en fin de ADVX1' cc STOP c ENDIF c 9998 CONTINUE c C ---------- bouclage cyclique DO itrac=1,ntra DO l = 1,llm DO j = lati,latf SM(iip1,j,l) = SM(1,j,l) S0(iip1,j,l,itrac) = S0(1,j,l,itrac) sx(iip1,j,l,itrac) = sx(1,j,l,itrac) sy(iip1,j,l,itrac) = sy(1,j,l,itrac) sz(iip1,j,l,itrac) = sz(1,j,l,itrac) END DO END DO ENDDO c ----------- qqtite totale de traceur dans tte l'atmosphere DO l = 1, llm DO j = 1, jjp1 DO i = 1, iim cIM 240405 sqf = sqf + S0(i,j,l,9) sqf = sqf + S0(i,j,l,ntra) END DO END DO END DO c PRINT*,'------ DIAG DANS ADVX - SORTIE -----' PRINT*,'sqf=',sqf c------------- RETURN END C_________________________________________________________________ C_________________________________________________________________