Changeset 5246 for LMDZ6/trunk/libf/dyn3d_common/diagedyn.f90
- Timestamp:
- Oct 21, 2024, 2:58:45 PM (22 hours ago)
- File:
-
- 1 moved
-
LMDZ6/trunk/libf/dyn3d_common/diagedyn.f90 (moved) (moved from LMDZ6/trunk/libf/dyn3d_common/diagedyn.F) (1 diff)
Legend:
- Unmodified
- Added
- Removed
-
LMDZ6/trunk/libf/dyn3d_common/diagedyn.f90
r5245 r5246 3 3 ! 4 4 5 C====================================================================== 6 SUBROUTINE diagedyn(tit,iprt,idiag,idiag2,dtime 7 e , ucov , vcov , ps, p ,pk , teta , q, ql) 8 C====================================================================== 9 C 10 C Purpose: 11 C Calcul la difference d'enthalpie et de masse d'eau entre 2 appels, 12 C et calcul le flux de chaleur et le flux d'eau necessaire a ces 13 C changements. Ces valeurs sont moyennees sur la surface de tout 14 C le globe et sont exprime en W/2 et kg/s/m2 15 C Outil pour diagnostiquer la conservation de l'energie 16 C et de la masse dans la dynamique. 17 C 18 C 19 c====================================================================== 20 C Arguments: 21 C tit-----imput-A15- Comment added in PRINT (CHARACTER*15) 22 C iprt----input-I- PRINT level ( <=1 : no PRINT) 23 C idiag---input-I- indice dans lequel sera range les nouveaux 24 C bilans d' entalpie et de masse 25 C idiag2--input-I-les nouveaux bilans d'entalpie et de masse 26 C sont compare au bilan de d'enthalpie de masse de 27 C l'indice numero idiag2 28 C Cas parriculier : si idiag2=0, pas de comparaison, on 29 c sort directement les bilans d'enthalpie et de masse 30 C dtime----input-R- time step (s) 31 C uconv, vconv-input-R- vents covariants (m/s) 32 C ps-------input-R- Surface pressure (Pa) 33 C p--------input-R- pressure at the interfaces 34 C pk-------input-R- pk= (p/Pref)**kappa 35 c teta-----input-R- potential temperature (K) 36 c q--------input-R- vapeur d'eau (kg/kg) 37 c ql-------input-R- liquid watter (kg/kg) 38 c aire-----input-R- mesh surafce (m2) 39 c 40 C the following total value are computed by UNIT of earth surface 41 C 42 C d_h_vcol--output-R- Heat flux (W/m2) define as the Enthalpy 43 c change (J/m2) during one time step (dtime) for the whole 44 C atmosphere (air, watter vapour, liquid and solid) 45 C d_qt------output-R- total water mass flux (kg/m2/s) defined as the 46 C total watter (kg/m2) change during one time step (dtime), 47 C d_qw------output-R- same, for the watter vapour only (kg/m2/s) 48 C d_ql------output-R- same, for the liquid watter only (kg/m2/s) 49 C d_ec------output-R- Cinetic Energy Budget (W/m2) for vertical air column 50 C 51 C 52 C J.L. Dufresne, July 2002 53 c====================================================================== 54 55 USE control_mod, ONLY : planet_type 56 57 IMPLICIT NONE 58 C 59 INCLUDE "dimensions.h" 60 INCLUDE "paramet.h" 61 INCLUDE "comgeom.h" 62 INCLUDE "iniprint.h" 63 64 !#ifdef CPP_EARTH 65 ! INCLUDE "../phylmd/YOMCST.h" 66 ! INCLUDE "../phylmd/YOETHF.h" 67 !#endif 68 ! Ehouarn: for now set these parameters to what is in Earth physics... 69 ! (cf ../phylmd/suphel.h) 70 ! this should be generalized... 71 REAL,PARAMETER :: RCPD= 72 & 3.5*(1000.*(6.0221367E+23*1.380658E-23)/28.9644) 73 REAL,PARAMETER :: RCPV= 74 & 4.*(1000.*(6.0221367E+23*1.380658E-23)/18.0153) 75 REAL,PARAMETER :: RCS=RCPV 76 REAL,PARAMETER :: RCW=RCPV 77 REAL,PARAMETER :: RLSTT=2.8345E+6 78 REAL,PARAMETER :: RLVTT=2.5008E+6 79 ! 80 C 81 INTEGER imjmp1 82 PARAMETER( imjmp1=iim*jjp1) 83 c Input variables 84 CHARACTER*15 tit 85 INTEGER iprt,idiag, idiag2 86 REAL dtime 87 REAL vcov(ip1jm,llm),ucov(ip1jmp1,llm) ! vents covariants 88 REAL ps(ip1jmp1) ! pression au sol 89 REAL p (ip1jmp1,llmp1 ) ! pression aux interfac.des couches 90 REAL pk (ip1jmp1,llm ) ! = (p/Pref)**kappa 91 REAL teta(ip1jmp1,llm) ! temperature potentielle 92 REAL q(ip1jmp1,llm) ! champs eau vapeur 93 REAL ql(ip1jmp1,llm) ! champs eau liquide 94 95 96 c Output variables 97 REAL d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec 98 C 99 C Local variables 100 c 101 REAL h_vcol_tot, h_dair_tot, h_qw_tot, h_ql_tot 102 . , h_qs_tot, qw_tot, ql_tot, qs_tot , ec_tot 103 c h_vcol_tot-- total enthalpy of vertical air column 104 C (air with watter vapour, liquid and solid) (J/m2) 105 c h_dair_tot-- total enthalpy of dry air (J/m2) 106 c h_qw_tot---- total enthalpy of watter vapour (J/m2) 107 c h_ql_tot---- total enthalpy of liquid watter (J/m2) 108 c h_qs_tot---- total enthalpy of solid watter (J/m2) 109 c qw_tot------ total mass of watter vapour (kg/m2) 110 c ql_tot------ total mass of liquid watter (kg/m2) 111 c qs_tot------ total mass of solid watter (kg/m2) 112 c ec_tot------ total cinetic energy (kg/m2) 113 C 114 REAL masse(ip1jmp1,llm) ! masse d'air 115 REAL vcont(ip1jm,llm),ucont(ip1jmp1,llm) 116 REAL ecin(ip1jmp1,llm) 117 118 REAL zaire(imjmp1) 119 REAL zps(imjmp1) 120 REAL zairm(imjmp1,llm) 121 REAL zecin(imjmp1,llm) 122 REAL zpaprs(imjmp1,llm) 123 REAL zpk(imjmp1,llm) 124 REAL zt(imjmp1,llm) 125 REAL zh(imjmp1,llm) 126 REAL zqw(imjmp1,llm) 127 REAL zql(imjmp1,llm) 128 REAL zqs(imjmp1,llm) 129 130 REAL zqw_col(imjmp1) 131 REAL zql_col(imjmp1) 132 REAL zqs_col(imjmp1) 133 REAL zec_col(imjmp1) 134 REAL zh_dair_col(imjmp1) 135 REAL zh_qw_col(imjmp1), zh_ql_col(imjmp1), zh_qs_col(imjmp1) 136 C 137 REAL d_h_dair, d_h_qw, d_h_ql, d_h_qs 138 C 139 REAL airetot, zcpvap, zcwat, zcice 140 C 141 INTEGER i, k, jj, ij , l ,ip1jjm1 142 C 143 INTEGER ndiag ! max number of diagnostic in parallel 144 PARAMETER (ndiag=10) 145 integer pas(ndiag) 146 save pas 147 data pas/ndiag*0/ 148 C 149 REAL h_vcol_pre(ndiag), h_dair_pre(ndiag), h_qw_pre(ndiag) 150 $ , h_ql_pre(ndiag), h_qs_pre(ndiag), qw_pre(ndiag) 151 $ , ql_pre(ndiag), qs_pre(ndiag) , ec_pre(ndiag) 152 SAVE h_vcol_pre, h_dair_pre, h_qw_pre, h_ql_pre 153 $ , h_qs_pre, qw_pre, ql_pre, qs_pre , ec_pre 154 155 156 !#ifdef CPP_EARTH 157 IF (planet_type=="earth") THEN 158 159 c====================================================================== 160 C Compute Kinetic enrgy 161 CALL covcont ( llm , ucov , vcov , ucont, vcont ) 162 CALL enercin ( vcov , ucov , vcont , ucont , ecin ) 163 CALL massdair( p, masse ) 164 c====================================================================== 165 C 166 C 167 print*,'MAIS POURQUOI DONC DIAGEDYN NE MARCHE PAS ?' 168 return 169 C On ne garde les donnees que dans les colonnes i=1,iim 170 DO jj = 1,jjp1 171 ip1jjm1=iip1*(jj-1) 172 DO ij = 1,iim 173 i=iim*(jj-1)+ij 174 zaire(i)=aire(ij+ip1jjm1) 175 zps(i)=ps(ij+ip1jjm1) 176 ENDDO 177 ENDDO 178 C 3D arrays 179 DO l = 1, llm 180 DO jj = 1,jjp1 181 ip1jjm1=iip1*(jj-1) 182 DO ij = 1,iim 183 i=iim*(jj-1)+ij 184 zairm(i,l) = masse(ij+ip1jjm1,l) 185 zecin(i,l) = ecin(ij+ip1jjm1,l) 186 zpaprs(i,l) = p(ij+ip1jjm1,l) 187 zpk(i,l) = pk(ij+ip1jjm1,l) 188 zh(i,l) = teta(ij+ip1jjm1,l) 189 zqw(i,l) = q(ij+ip1jjm1,l) 190 zql(i,l) = ql(ij+ip1jjm1,l) 191 zqs(i,l) = 0. 192 ENDDO 193 ENDDO 194 ENDDO 195 C 196 C Reset variables 197 DO i = 1, imjmp1 198 zqw_col(i)=0. 199 zql_col(i)=0. 200 zqs_col(i)=0. 201 zec_col(i) = 0. 202 zh_dair_col(i) = 0. 203 zh_qw_col(i) = 0. 204 zh_ql_col(i) = 0. 205 zh_qs_col(i) = 0. 5 !====================================================================== 6 SUBROUTINE diagedyn(tit,iprt,idiag,idiag2,dtime & 7 , ucov , vcov , ps, p ,pk , teta , q, ql) 8 !====================================================================== 9 ! 10 ! Purpose: 11 ! Calcul la difference d'enthalpie et de masse d'eau entre 2 appels, 12 ! et calcul le flux de chaleur et le flux d'eau necessaire a ces 13 ! changements. Ces valeurs sont moyennees sur la surface de tout 14 ! le globe et sont exprime en W/2 et kg/s/m2 15 ! Outil pour diagnostiquer la conservation de l'energie 16 ! et de la masse dans la dynamique. 17 ! 18 ! 19 !====================================================================== 20 ! Arguments: 21 ! tit-----imput-A15- Comment added in PRINT (CHARACTER*15) 22 ! iprt----input-I- PRINT level ( <=1 : no PRINT) 23 ! idiag---input-I- indice dans lequel sera range les nouveaux 24 ! bilans d' entalpie et de masse 25 ! idiag2--input-I-les nouveaux bilans d'entalpie et de masse 26 ! sont compare au bilan de d'enthalpie de masse de 27 ! l'indice numero idiag2 28 ! Cas parriculier : si idiag2=0, pas de comparaison, on 29 ! sort directement les bilans d'enthalpie et de masse 30 ! dtime----input-R- time step (s) 31 ! uconv, vconv-input-R- vents covariants (m/s) 32 ! ps-------input-R- Surface pressure (Pa) 33 ! p--------input-R- pressure at the interfaces 34 ! pk-------input-R- pk= (p/Pref)**kappa 35 ! teta-----input-R- potential temperature (K) 36 ! q--------input-R- vapeur d'eau (kg/kg) 37 ! ql-------input-R- liquid watter (kg/kg) 38 ! aire-----input-R- mesh surafce (m2) 39 ! 40 ! the following total value are computed by UNIT of earth surface 41 ! 42 ! d_h_vcol--output-R- Heat flux (W/m2) define as the Enthalpy 43 ! change (J/m2) during one time step (dtime) for the whole 44 ! atmosphere (air, watter vapour, liquid and solid) 45 ! d_qt------output-R- total water mass flux (kg/m2/s) defined as the 46 ! total watter (kg/m2) change during one time step (dtime), 47 ! d_qw------output-R- same, for the watter vapour only (kg/m2/s) 48 ! d_ql------output-R- same, for the liquid watter only (kg/m2/s) 49 ! d_ec------output-R- Cinetic Energy Budget (W/m2) for vertical air column 50 ! 51 ! 52 ! J.L. Dufresne, July 2002 53 !====================================================================== 54 55 USE control_mod, ONLY : planet_type 56 57 IMPLICIT NONE 58 ! 59 INCLUDE "dimensions.h" 60 INCLUDE "paramet.h" 61 INCLUDE "comgeom.h" 62 INCLUDE "iniprint.h" 63 64 !#ifdef CPP_EARTH 65 ! INCLUDE "../phylmd/YOMCST.h" 66 ! INCLUDE "../phylmd/YOETHF.h" 67 !#endif 68 ! Ehouarn: for now set these parameters to what is in Earth physics... 69 ! (cf ../phylmd/suphel.h) 70 ! this should be generalized... 71 REAL,PARAMETER :: RCPD= & 72 3.5*(1000.*(6.0221367E+23*1.380658E-23)/28.9644) 73 REAL,PARAMETER :: RCPV= & 74 4.*(1000.*(6.0221367E+23*1.380658E-23)/18.0153) 75 REAL,PARAMETER :: RCS=RCPV 76 REAL,PARAMETER :: RCW=RCPV 77 REAL,PARAMETER :: RLSTT=2.8345E+6 78 REAL,PARAMETER :: RLVTT=2.5008E+6 79 ! 80 ! 81 INTEGER :: imjmp1 82 PARAMETER( imjmp1=iim*jjp1) 83 ! Input variables 84 CHARACTER(len=15) :: tit 85 INTEGER :: iprt,idiag, idiag2 86 REAL :: dtime 87 REAL :: vcov(ip1jm,llm),ucov(ip1jmp1,llm) ! vents covariants 88 REAL :: ps(ip1jmp1) ! pression au sol 89 REAL :: p (ip1jmp1,llmp1 ) ! pression aux interfac.des couches 90 REAL :: pk (ip1jmp1,llm ) ! = (p/Pref)**kappa 91 REAL :: teta(ip1jmp1,llm) ! temperature potentielle 92 REAL :: q(ip1jmp1,llm) ! champs eau vapeur 93 REAL :: ql(ip1jmp1,llm) ! champs eau liquide 94 95 96 ! Output variables 97 REAL :: d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec 98 ! 99 ! Local variables 100 ! 101 REAL :: h_vcol_tot, h_dair_tot, h_qw_tot, h_ql_tot & 102 , h_qs_tot, qw_tot, ql_tot, qs_tot , ec_tot 103 ! h_vcol_tot-- total enthalpy of vertical air column 104 ! (air with watter vapour, liquid and solid) (J/m2) 105 ! h_dair_tot-- total enthalpy of dry air (J/m2) 106 ! h_qw_tot---- total enthalpy of watter vapour (J/m2) 107 ! h_ql_tot---- total enthalpy of liquid watter (J/m2) 108 ! h_qs_tot---- total enthalpy of solid watter (J/m2) 109 ! qw_tot------ total mass of watter vapour (kg/m2) 110 ! ql_tot------ total mass of liquid watter (kg/m2) 111 ! qs_tot------ total mass of solid watter (kg/m2) 112 ! ec_tot------ total cinetic energy (kg/m2) 113 ! 114 REAL :: masse(ip1jmp1,llm) ! masse d'air 115 REAL :: vcont(ip1jm,llm),ucont(ip1jmp1,llm) 116 REAL :: ecin(ip1jmp1,llm) 117 118 REAL :: zaire(imjmp1) 119 REAL :: zps(imjmp1) 120 REAL :: zairm(imjmp1,llm) 121 REAL :: zecin(imjmp1,llm) 122 REAL :: zpaprs(imjmp1,llm) 123 REAL :: zpk(imjmp1,llm) 124 REAL :: zt(imjmp1,llm) 125 REAL :: zh(imjmp1,llm) 126 REAL :: zqw(imjmp1,llm) 127 REAL :: zql(imjmp1,llm) 128 REAL :: zqs(imjmp1,llm) 129 130 REAL :: zqw_col(imjmp1) 131 REAL :: zql_col(imjmp1) 132 REAL :: zqs_col(imjmp1) 133 REAL :: zec_col(imjmp1) 134 REAL :: zh_dair_col(imjmp1) 135 REAL :: zh_qw_col(imjmp1), zh_ql_col(imjmp1), zh_qs_col(imjmp1) 136 ! 137 REAL :: d_h_dair, d_h_qw, d_h_ql, d_h_qs 138 ! 139 REAL :: airetot, zcpvap, zcwat, zcice 140 ! 141 INTEGER :: i, k, jj, ij , l ,ip1jjm1 142 ! 143 INTEGER :: ndiag ! max number of diagnostic in parallel 144 PARAMETER (ndiag=10) 145 integer :: pas(ndiag) 146 save pas 147 data pas/ndiag*0/ 148 ! 149 REAL :: h_vcol_pre(ndiag), h_dair_pre(ndiag), h_qw_pre(ndiag) & 150 , h_ql_pre(ndiag), h_qs_pre(ndiag), qw_pre(ndiag) & 151 , ql_pre(ndiag), qs_pre(ndiag) , ec_pre(ndiag) 152 SAVE h_vcol_pre, h_dair_pre, h_qw_pre, h_ql_pre & 153 , h_qs_pre, qw_pre, ql_pre, qs_pre , ec_pre 154 155 156 !#ifdef CPP_EARTH 157 IF (planet_type=="earth") THEN 158 159 !====================================================================== 160 ! Compute Kinetic enrgy 161 CALL covcont ( llm , ucov , vcov , ucont, vcont ) 162 CALL enercin ( vcov , ucov , vcont , ucont , ecin ) 163 CALL massdair( p, masse ) 164 !====================================================================== 165 ! 166 ! 167 print*,'MAIS POURQUOI DONC DIAGEDYN NE MARCHE PAS ?' 168 return 169 ! On ne garde les donnees que dans les colonnes i=1,iim 170 DO jj = 1,jjp1 171 ip1jjm1=iip1*(jj-1) 172 DO ij = 1,iim 173 i=iim*(jj-1)+ij 174 zaire(i)=aire(ij+ip1jjm1) 175 zps(i)=ps(ij+ip1jjm1) 176 ENDDO 177 ENDDO 178 ! 3D arrays 179 DO l = 1, llm 180 DO jj = 1,jjp1 181 ip1jjm1=iip1*(jj-1) 182 DO ij = 1,iim 183 i=iim*(jj-1)+ij 184 zairm(i,l) = masse(ij+ip1jjm1,l) 185 zecin(i,l) = ecin(ij+ip1jjm1,l) 186 zpaprs(i,l) = p(ij+ip1jjm1,l) 187 zpk(i,l) = pk(ij+ip1jjm1,l) 188 zh(i,l) = teta(ij+ip1jjm1,l) 189 zqw(i,l) = q(ij+ip1jjm1,l) 190 zql(i,l) = ql(ij+ip1jjm1,l) 191 zqs(i,l) = 0. 206 192 ENDDO 207 C 208 zcpvap=RCPV 209 zcwat=RCW 210 zcice=RCS 211 C 212 C Compute vertical sum for each atmospheric column 213 C ================================================ 214 DO k = 1, llm 215 DO i = 1, imjmp1 216 C Watter mass 217 zqw_col(i) = zqw_col(i) + zqw(i,k)*zairm(i,k) 218 zql_col(i) = zql_col(i) + zql(i,k)*zairm(i,k) 219 zqs_col(i) = zqs_col(i) + zqs(i,k)*zairm(i,k) 220 C Cinetic Energy 221 zec_col(i) = zec_col(i) 222 $ +zecin(i,k)*zairm(i,k) 223 C Air enthalpy 224 zt(i,k)= zh(i,k) * zpk(i,k) / RCPD 225 zh_dair_col(i) = zh_dair_col(i) 226 $ + RCPD*(1.-zqw(i,k)-zql(i,k)-zqs(i,k))*zairm(i,k)*zt(i,k) 227 zh_qw_col(i) = zh_qw_col(i) 228 $ + zcpvap*zqw(i,k)*zairm(i,k)*zt(i,k) 229 zh_ql_col(i) = zh_ql_col(i) 230 $ + zcwat*zql(i,k)*zairm(i,k)*zt(i,k) 231 $ - RLVTT*zql(i,k)*zairm(i,k) 232 zh_qs_col(i) = zh_qs_col(i) 233 $ + zcice*zqs(i,k)*zairm(i,k)*zt(i,k) 234 $ - RLSTT*zqs(i,k)*zairm(i,k) 235 236 END DO 237 ENDDO 238 C 239 C Mean over the planete surface 240 C ============================= 241 qw_tot = 0. 242 ql_tot = 0. 243 qs_tot = 0. 244 ec_tot = 0. 245 h_vcol_tot = 0. 246 h_dair_tot = 0. 247 h_qw_tot = 0. 248 h_ql_tot = 0. 249 h_qs_tot = 0. 250 airetot=0. 251 C 252 do i=1,imjmp1 253 qw_tot = qw_tot + zqw_col(i) 254 ql_tot = ql_tot + zql_col(i) 255 qs_tot = qs_tot + zqs_col(i) 256 ec_tot = ec_tot + zec_col(i) 257 h_dair_tot = h_dair_tot + zh_dair_col(i) 258 h_qw_tot = h_qw_tot + zh_qw_col(i) 259 h_ql_tot = h_ql_tot + zh_ql_col(i) 260 h_qs_tot = h_qs_tot + zh_qs_col(i) 261 airetot=airetot+zaire(i) 262 END DO 263 C 264 qw_tot = qw_tot/airetot 265 ql_tot = ql_tot/airetot 266 qs_tot = qs_tot/airetot 267 ec_tot = ec_tot/airetot 268 h_dair_tot = h_dair_tot/airetot 269 h_qw_tot = h_qw_tot/airetot 270 h_ql_tot = h_ql_tot/airetot 271 h_qs_tot = h_qs_tot/airetot 272 C 273 h_vcol_tot = h_dair_tot+h_qw_tot+h_ql_tot+h_qs_tot 274 C 275 C Compute the change of the atmospheric state compare to the one 276 C stored in "idiag2", and convert it in flux. THis computation 277 C is performed IF idiag2 /= 0 and IF it is not the first CALL 278 c for "idiag" 279 C =================================== 280 C 281 IF ( (idiag2.gt.0) .and. (pas(idiag2) .ne. 0) ) THEN 282 d_h_vcol = (h_vcol_tot - h_vcol_pre(idiag2) )/dtime 283 d_h_dair = (h_dair_tot- h_dair_pre(idiag2))/dtime 284 d_h_qw = (h_qw_tot - h_qw_pre(idiag2) )/dtime 285 d_h_ql = (h_ql_tot - h_ql_pre(idiag2) )/dtime 286 d_h_qs = (h_qs_tot - h_qs_pre(idiag2) )/dtime 287 d_qw = (qw_tot - qw_pre(idiag2) )/dtime 288 d_ql = (ql_tot - ql_pre(idiag2) )/dtime 289 d_qs = (qs_tot - qs_pre(idiag2) )/dtime 290 d_ec = (ec_tot - ec_pre(idiag2) )/dtime 291 d_qt = d_qw + d_ql + d_qs 292 ELSE 293 d_h_vcol = 0. 294 d_h_dair = 0. 295 d_h_qw = 0. 296 d_h_ql = 0. 297 d_h_qs = 0. 298 d_qw = 0. 299 d_ql = 0. 300 d_qs = 0. 301 d_ec = 0. 302 d_qt = 0. 303 ENDIF 304 C 305 IF (iprt.ge.2) THEN 306 WRITE(6,9000) tit,pas(idiag),d_qt,d_qw,d_ql,d_qs 307 9000 format('Dyn3d. Watter Mass Budget (kg/m2/s)',A15 308 $ ,1i6,10(1pE14.6)) 309 WRITE(6,9001) tit,pas(idiag), d_h_vcol 193 ENDDO 194 ENDDO 195 ! 196 ! Reset variables 197 DO i = 1, imjmp1 198 zqw_col(i)=0. 199 zql_col(i)=0. 200 zqs_col(i)=0. 201 zec_col(i) = 0. 202 zh_dair_col(i) = 0. 203 zh_qw_col(i) = 0. 204 zh_ql_col(i) = 0. 205 zh_qs_col(i) = 0. 206 ENDDO 207 ! 208 zcpvap=RCPV 209 zcwat=RCW 210 zcice=RCS 211 ! 212 ! Compute vertical sum for each atmospheric column 213 ! ================================================ 214 DO k = 1, llm 215 DO i = 1, imjmp1 216 ! Watter mass 217 zqw_col(i) = zqw_col(i) + zqw(i,k)*zairm(i,k) 218 zql_col(i) = zql_col(i) + zql(i,k)*zairm(i,k) 219 zqs_col(i) = zqs_col(i) + zqs(i,k)*zairm(i,k) 220 ! Cinetic Energy 221 zec_col(i) = zec_col(i) & 222 +zecin(i,k)*zairm(i,k) 223 ! Air enthalpy 224 zt(i,k)= zh(i,k) * zpk(i,k) / RCPD 225 zh_dair_col(i) = zh_dair_col(i) & 226 + RCPD*(1.-zqw(i,k)-zql(i,k)-zqs(i,k))*zairm(i,k)*zt(i,k) 227 zh_qw_col(i) = zh_qw_col(i) & 228 + zcpvap*zqw(i,k)*zairm(i,k)*zt(i,k) 229 zh_ql_col(i) = zh_ql_col(i) & 230 + zcwat*zql(i,k)*zairm(i,k)*zt(i,k) & 231 - RLVTT*zql(i,k)*zairm(i,k) 232 zh_qs_col(i) = zh_qs_col(i) & 233 + zcice*zqs(i,k)*zairm(i,k)*zt(i,k) & 234 - RLSTT*zqs(i,k)*zairm(i,k) 235 236 END DO 237 ENDDO 238 ! 239 ! Mean over the planete surface 240 ! ============================= 241 qw_tot = 0. 242 ql_tot = 0. 243 qs_tot = 0. 244 ec_tot = 0. 245 h_vcol_tot = 0. 246 h_dair_tot = 0. 247 h_qw_tot = 0. 248 h_ql_tot = 0. 249 h_qs_tot = 0. 250 airetot=0. 251 ! 252 do i=1,imjmp1 253 qw_tot = qw_tot + zqw_col(i) 254 ql_tot = ql_tot + zql_col(i) 255 qs_tot = qs_tot + zqs_col(i) 256 ec_tot = ec_tot + zec_col(i) 257 h_dair_tot = h_dair_tot + zh_dair_col(i) 258 h_qw_tot = h_qw_tot + zh_qw_col(i) 259 h_ql_tot = h_ql_tot + zh_ql_col(i) 260 h_qs_tot = h_qs_tot + zh_qs_col(i) 261 airetot=airetot+zaire(i) 262 END DO 263 ! 264 qw_tot = qw_tot/airetot 265 ql_tot = ql_tot/airetot 266 qs_tot = qs_tot/airetot 267 ec_tot = ec_tot/airetot 268 h_dair_tot = h_dair_tot/airetot 269 h_qw_tot = h_qw_tot/airetot 270 h_ql_tot = h_ql_tot/airetot 271 h_qs_tot = h_qs_tot/airetot 272 ! 273 h_vcol_tot = h_dair_tot+h_qw_tot+h_ql_tot+h_qs_tot 274 ! 275 ! Compute the change of the atmospheric state compare to the one 276 ! stored in "idiag2", and convert it in flux. THis computation 277 ! is performed IF idiag2 /= 0 and IF it is not the first CALL 278 ! for "idiag" 279 ! =================================== 280 ! 281 IF ( (idiag2.gt.0) .and. (pas(idiag2) .ne. 0) ) THEN 282 d_h_vcol = (h_vcol_tot - h_vcol_pre(idiag2) )/dtime 283 d_h_dair = (h_dair_tot- h_dair_pre(idiag2))/dtime 284 d_h_qw = (h_qw_tot - h_qw_pre(idiag2) )/dtime 285 d_h_ql = (h_ql_tot - h_ql_pre(idiag2) )/dtime 286 d_h_qs = (h_qs_tot - h_qs_pre(idiag2) )/dtime 287 d_qw = (qw_tot - qw_pre(idiag2) )/dtime 288 d_ql = (ql_tot - ql_pre(idiag2) )/dtime 289 d_qs = (qs_tot - qs_pre(idiag2) )/dtime 290 d_ec = (ec_tot - ec_pre(idiag2) )/dtime 291 d_qt = d_qw + d_ql + d_qs 292 ELSE 293 d_h_vcol = 0. 294 d_h_dair = 0. 295 d_h_qw = 0. 296 d_h_ql = 0. 297 d_h_qs = 0. 298 d_qw = 0. 299 d_ql = 0. 300 d_qs = 0. 301 d_ec = 0. 302 d_qt = 0. 303 ENDIF 304 ! 305 IF (iprt.ge.2) THEN 306 WRITE(6,9000) tit,pas(idiag),d_qt,d_qw,d_ql,d_qs 307 9000 format('Dyn3d. Watter Mass Budget (kg/m2/s)',A15 & 308 ,1i6,10(1pE14.6)) 309 WRITE(6,9001) tit,pas(idiag), d_h_vcol 310 310 9001 format('Dyn3d. Enthalpy Budget (W/m2) ',A15,1i6,10(F8.2)) 311 311 WRITE(6,9002) tit,pas(idiag), d_ec 312 312 9002 format('Dyn3d. Cinetic Energy Budget (W/m2) ',A15,1i6,10(F8.2)) 313 CWRITE(6,9003) tit,pas(idiag), ec_tot313 ! WRITE(6,9003) tit,pas(idiag), ec_tot 314 314 9003 format('Dyn3d. Cinetic Energy (W/m2) ',A15,1i6,10(E15.6)) 315 315 WRITE(6,9004) tit,pas(idiag), d_h_vcol+d_ec 316 316 9004 format('Dyn3d. Total Energy Budget (W/m2) ',A15,1i6,10(F8.2)) 317 END IF318 C 319 CStore the new atmospheric state in "idiag"320 C 321 322 323 324 325 326 327 328 329 330 331 C 332 !#else333 334 335 336 !#endif337 ! #endif of #ifdef CPP_EARTH 338 RETURN339 END 317 END IF 318 ! 319 ! Store the new atmospheric state in "idiag" 320 ! 321 pas(idiag)=pas(idiag)+1 322 h_vcol_pre(idiag) = h_vcol_tot 323 h_dair_pre(idiag) = h_dair_tot 324 h_qw_pre(idiag) = h_qw_tot 325 h_ql_pre(idiag) = h_ql_tot 326 h_qs_pre(idiag) = h_qs_tot 327 qw_pre(idiag) = qw_tot 328 ql_pre(idiag) = ql_tot 329 qs_pre(idiag) = qs_tot 330 ec_pre (idiag) = ec_tot 331 ! 332 !#else 333 ELSE 334 write(lunout,*)'diagedyn: set to function with Earth parameters' 335 ENDIF ! of if (planet_type=="earth") 336 !#endif 337 ! #endif of #ifdef CPP_EARTH 338 RETURN 339 END SUBROUTINE diagedyn
Note: See TracChangeset
for help on using the changeset viewer.
