Changeset 1310 for trunk/LMDZ.VENUS/libf/phyvenus/physiq.F
- Timestamp:
- Jul 10, 2014, 5:59:16 PM (10 years ago)
- File:
-
- 1 edited
-
trunk/LMDZ.VENUS/libf/phyvenus/physiq.F (modified) (10 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/LMDZ.VENUS/libf/phyvenus/physiq.F
r1306 r1310 56 56 USE ioipsl 57 57 ! USE histcom ! not needed; histcom is included in ioipsl 58 USE chemparam_mod59 58 USE infotrac 60 59 USE control_mod … … 66 65 USE iophy 67 66 use cpdet_mod, only: cpdet, t2tpot 67 USE chemparam_mod 68 USE conc 69 USE compo_hedin83_mod2 70 use moyzon_mod, only: tmoy 68 71 use ieee_arithmetic 69 72 IMPLICIT none … … 86 89 #include "logic.h" 87 90 #include "tabcontrol.h" 91 #include "nirdata.h" 92 #include "hedin.h" 88 93 c====================================================================== 89 94 LOGICAL ok_journe ! sortir le fichier journalier … … 223 228 EXTERNAL mucorr 224 229 EXTERNAL phytrac 230 EXTERNAL nirco2abs 231 EXTERNAL nir_leedat 232 EXTERNAL nltecool 233 EXTERNAL nlte_tcool 234 EXTERNAL nlte_setup 235 EXTERNAL blendrad 236 EXTERNAL nlthermeq 237 EXTERNAL euvheat 238 EXTERNAL param_read 239 EXTERNAL param_read_e107 240 EXTERNAL conduction 241 EXTERNAL molvis 225 242 c 226 243 c Variables locales … … 247 264 REAL zphi(klon,klev) 248 265 REAL zzlev(klon,klev+1),zzlay(klon,klev),z1,z2 266 real tsurf(klon) 267 268 c va avec nlte_tcool 269 INTEGER ierr_nlte 270 REAL varerr 249 271 250 272 c Variables du changement … … 268 290 REAL d_u_hin(klon,klev), d_v_hin(klon,klev) 269 291 REAL d_t_hin(klon,klev) 292 293 c Tendencies due to radiative scheme [K/s] 294 c d_t_rad,dtsw,dtlw,d_t_nirco2,d_t_nlte,d_t_euv 295 c are not computed at each physical timestep 296 c therefore, they are defined and saved in phys_state_var_mod 297 298 c Tendencies due to molecular viscosity and conduction 299 real d_t_conduc(klon,klev) ! [K/s] 300 real d_u_molvis(klon,klev) ! (m/s) /s 301 real d_v_molvis(klon,klev) ! (m/s) /s 270 302 271 303 c … … 386 418 call phys_state_var_init 387 419 c 420 c Initialising Hedin model for upper atm 421 c (to be revised when coupled to chemistry) : 422 call conc_init 423 c 388 424 c Initialiser les compteurs: 389 425 c … … 458 494 C source dans couche limite 459 495 source = 0.0 ! pas de source, pour l'instant 460 C461 496 c--------- 497 498 c--------- 499 c INITIALIZE THERMOSPHERIC PARAMETERS 500 c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 501 502 if (callthermos) then 503 if(solvarmod.eq.0) call param_read 504 if(solvarmod.eq.1) call param_read_e107 505 endif 506 507 c Initialisation (recomputed in concentration2) 508 do ig=1,klon 509 do j=1,klev 510 rnew(ig,j)=R 511 cpnew(ig,j)=cpdet(tmoy(j)) 512 c print*, ' physique l503' 513 c print*, j, cpdet(tmoy(j)) 514 mmean(ig,j)=RMD 515 akknew(ig,j)=1.e-4 516 enddo 517 c stop 518 519 enddo 520 521 IF(callthermos.or.callnlte.or.callnirco2) THEN 522 call compo_hedin83_init2 523 ENDIF 524 if (callnlte) call nlte_setup 525 if(callnirco2.and.(nircorr.eq.1)) call nir_leedat 526 c--------- 527 462 528 c 463 529 c Verifications: … … 684 750 DO k=1,klev 685 751 DO i=1,klon 686 zzlay(i,k)=zphi(i,k)/RG 752 zzlay(i,k)=zphi(i,k)/RG ! [m] 687 753 ENDDO 688 754 ENDDO 689 755 DO i=1,klon 690 zzlev(i,1)=pphis(i)/RG 756 zzlev(i,1)=pphis(i)/RG ! [m] 691 757 ENDDO 692 758 DO k=2,klev … … 1042 1108 c==================================================================== 1043 1109 1110 c----------------------------------------------------------------------- 1111 c . Compute radiative tendencies : 1112 c------------------------------------ 1113 c==================================================================== 1044 1114 IF (MOD(itaprad,radpas).EQ.0) THEN 1045 c print*,'RAYONNEMENT ', 1046 c $ ' (itaprad=',itaprad,'/radpas=',radpas,')' 1115 c==================================================================== 1047 1116 1048 1117 dtimerad = dtime*REAL(radpas) ! pas de temps du rayonnement (s) 1049 1118 c PRINT*,'dtimerad,dtime,radpas',dtimerad,dtime,radpas 1050 1119 1120 c Calcul pour Cp rnew et mmean avec traceurs (Cp independant de T !! ) 1121 IF(callnlte.or.callthermos) THEN 1122 call compo_hedin83_mod(pplay,rmu0, 1123 & co2vmr_gcm,covmr_gcm,ovmr_gcm,n2vmr_gcm,nvmr_gcm) 1124 ENDIF 1125 1126 if(callthermos) then 1127 call concentrations2(pplay,t_seri,d_t,co2vmr_gcm, n2vmr_gcm, 1128 & covmr_gcm, ovmr_gcm,nvmr_gcm,pdtphys) 1129 endif 1130 1131 1132 c NLTE cooling from CO2 emission 1133 c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 1134 1135 IF(callnlte) THEN 1136 if(nltemodel.eq.0.or.nltemodel.eq.1) then 1137 CALL nltecool(klon, klev, pplay*9.869e-6, t_seri, 1138 $ co2vmr_gcm,n2vmr_gcm, covmr_gcm, ovmr_gcm, 1139 $ d_t_nlte) 1140 else if(nltemodel.eq.2) then 1141 CALL nlte_tcool(klon,klev,pplay*9.869e-6, 1142 $ t_seri,zzlay,co2vmr_gcm, n2vmr_gcm, covmr_gcm, 1143 $ ovmr_gcm,d_t_nlte,ierr_nlte,varerr ) 1144 if(ierr_nlte.gt.0) then 1145 write(*,*) 1146 $ 'WARNING: nlte_tcool output with error message', 1147 $ 'ierr_nlte=',ierr_nlte,'varerr=',varerr 1148 write(*,*)'I will continue anyway' 1149 endif 1150 1151 endif 1152 1153 ELSE 1154 1155 d_t_nlte(:,:)=0. 1156 1157 ENDIF 1158 1159 c Find number of layers for LTE radiation calculations 1160 1161 IF(callnlte .or. callnirco2) 1162 $ CALL nlthermeq(klon, klev, paprs, pplay) 1163 1164 1165 c LTE radiative transfert / solar / IR matrix 1166 c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 1167 1051 1168 CALL radlwsw 1052 1169 e (dist, rmu0, fract, zzlev, 1053 e paprs, pplay,ftsol, t_seri, 1054 s heat,cool,radsol, 1055 s topsw,toplw,solsw,sollw, 1056 s sollwdown, 1057 s lwnet, swnet) 1058 1059 itaprad = 0 1170 e paprs, pplay,ftsol, t_seri) 1171 1172 1173 c CO2 near infrared absorption 1174 c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 1175 1176 d_t_nirco2(:,:)=0. 1177 if (callnirco2) then 1178 call nirco2abs (klon, klev, pplay, dist, nqmax, qx, 1179 . rmu0, fract, d_t_nirco2, 1180 . co2vmr_gcm, ovmr_gcm) 1181 endif 1182 1183 1184 c Net atmospheric radiative heating rate (K.s-1) 1185 c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 1186 1187 IF(callnlte.or.callnirco2) THEN 1188 CALL blendrad(klon, klev, pplay,heat, 1189 & cool, d_t_nirco2,d_t_nlte, dtsw, dtlw) 1190 ELSE 1191 dtsw(:,:)=heat(:,:) 1192 dtlw(:,:)=-1*cool(:,:) 1193 ENDIF 1194 1195 DO k=1,klev 1196 DO i=1,klon 1197 d_t_rad(i,k) = dtsw(i,k) + dtlw(i,k) ! K/s 1198 ENDDO 1199 ENDDO 1200 1201 1202 cc--------------------------------------------- 1203 1204 c EUV heating rate (K.s-1) 1205 c ~~~~~~~~~~~~~~~~~~~~~~~~ 1206 1207 d_t_euv(:,:)=0. 1208 1209 IF (callthermos) THEN 1210 1211 c call thermosphere(zplev,zplay,dist_sol, 1212 c $ mu0,ptimestep,ptime,zday,tsurf,zzlev,zzlay, 1213 c & pt,pq,pu,pv,pdt,pdq, 1214 c $ zdteuv,zdtconduc,zdumolvis,zdvmolvis,zdqmoldiff) 1215 1216 call euvheat(klon, klev, t_seri,paprs,pplay,zzlay, 1217 $ rmu0,pdtphys,gmtime,rjourvrai, 1218 C $ pq,pdq,zdteuv) 1219 $ co2vmr_gcm, n2vmr_gcm, covmr_gcm, 1220 $ ovmr_gcm,nvmr_gcm,d_t_euv ) 1221 1222 DO k=1,klev 1223 DO ig=1,klon 1224 d_t_rad(ig,k)=d_t_rad(ig,k)+d_t_euv(ig,k) 1225 1226 ENDDO 1227 ENDDO 1228 1229 ENDIF ! callthermos 1230 1231 c==================================================================== 1232 itaprad = 0 1233 ENDIF ! radpas 1234 c==================================================================== 1235 c 1236 c Ajouter la tendance des rayonnements (tous les pas) 1237 c 1060 1238 DO k = 1, klev 1061 1239 DO i = 1, klon 1062 dtrad(i,k) = heat(i,k)-cool(i,k) !K/s 1063 ENDDO 1064 ENDDO 1065 1066 ENDIF 1067 itaprad = itaprad + 1 1068 c==================================================================== 1069 c 1070 c Ajouter la tendance des rayonnements (tous les pas) 1071 c 1072 DO k = 1, klev 1073 DO i = 1, klon 1074 t_seri(i,k) = t_seri(i,k) + dtrad(i,k) * dtime 1075 ENDDO 1076 ENDDO 1077 1240 t_seri(i,k) = t_seri(i,k) + d_t_rad(i,k) * dtime 1241 ENDDO 1242 ENDDO 1243 1244 ! CONDUCTION and MOLECULAR VISCOSITY 1245 c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 1246 1247 d_t_conduc(:,:)=0. 1248 d_u_molvis(:,:)=0. 1249 d_v_molvis(:,:)=0. 1250 1251 IF (callthermos) THEN 1252 1253 tsurf(:)=t_seri(:,1) 1254 call conduction(klon, klev,pdtphys, 1255 $ pplay,paprs,t_seri, 1256 $ tsurf,zzlev,zzlay,d_t_conduc) 1257 1258 call molvis(klon, klev,pdtphys, 1259 $ pplay,paprs,t_seri, 1260 $ u,tsurf,zzlev,zzlay,d_u_molvis) 1261 1262 call molvis(klon, klev, pdtphys, 1263 $ pplay,paprs,t_seri, 1264 $ v,tsurf,zzlev,zzlay,d_u_molvis) 1265 1266 DO k=1,klev 1267 DO ig=1,klon 1268 t_seri(ig,k)= t_seri(ig,k)+ d_t_conduc(ig,k)*dtime ! [K] 1269 u_seri(ig,k)= u_seri(ig,k)+ d_u_molvis(ig,k)*dtime ! m/s 1270 v_seri(ig,k)= v_seri(ig,k)+ d_v_molvis(ig,k)*dtime ! m/s 1271 ENDDO 1272 ENDDO 1273 1274 ENDIF ! callthermos 1275 1276 c==================================================================== 1078 1277 ! tests: output tendencies 1079 1278 ! call writefield_phy('physiq_dtrad',dtrad,klev)
Note: See TracChangeset
for help on using the changeset viewer.
