- Timestamp:
- Oct 18, 2024, 10:29:59 AM (2 months ago)
- Location:
- trunk/LMDZ.MARS
- Files:
-
- 1 added
- 10 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/LMDZ.MARS/changelog.txt
r3459 r3464 4709 4709 - Addition in the "util" folder of two python scripts: one is to visualize easily any variable of a NetCDF file; the other is to display useful information about the variables of a NetCDF file to help for debugging. 4710 4710 - Move of the launching scripts for the 1D model from the "deftank" to the "util" folder. 4711 4712 == 18/10/2024 == EM 4713 Some tidying in aeronomars: 4714 - make a jthermcalc_util.F to contain utility routines (used by jthermcal & 4715 jthermcalc_e107). Also add the possibility (turned off by default) in the 4716 interfast routine to do extra sanity checks. 4717 - turn chemthermos, euvheat, hrtherm, jthermcalc, jthermcalc_e107, 4718 paramphoto_compact and photochemistry into modules -
trunk/LMDZ.MARS/libf/aeronomars/calchim_mod.F90
r3461 r3464 31 31 use conc_mod, only: mmean ! mean molecular mass of the atmosphere 32 32 use comcstfi_h, only: pi 33 use chemthermos_mod, only: chemthermos 34 use photochemistry_mod, only: photochemistry 33 35 use chemistrydata_mod, only: read_phototable 34 36 use photolysis_mod, only: init_photolysis, nphot -
trunk/LMDZ.MARS/libf/aeronomars/chemthermos.F90
r2615 r3464 1 MODULE chemthermos_mod 2 3 IMPLICIT NONE 4 5 CONTAINS 6 1 7 SUBROUTINE chemthermos(ig,nlayer,lswitch,chemthermod,zycol,ztemp, & 2 8 zdens,zpress,zlocal,zenit,ptimestep,zday,em_no,em_o2) … … 12 18 13 19 use param_v4_h, only: Pno, Po2 14 USE comcstfi_h 20 use jthermcalc_e107_mod, only: jthermcalc_e107 21 use paramfoto_compact_mod, only: paramfoto_compact 22 15 23 IMPLICIT NONE 16 24 !======================================================================= … … 32 40 ! ------------------ 33 41 ! 34 #include "callkeys.h"42 include "callkeys.h" 35 43 !----------------------------------------------------------------------- 36 44 ! Input/Output … … 540 548 deallocate(rm) 541 549 542 return543 end 544 545 546 547 550 end subroutine chemthermos 551 552 END MODULE chemthermos_mod 553 554 555 -
trunk/LMDZ.MARS/libf/aeronomars/euvheat.F90
r3158 r3464 1 MODULE euvheat_mod 2 3 IMPLICIT NONE 4 5 CONTAINS 6 1 7 SUBROUTINE euvheat(ngrid,nlayer,nq,pt,pdt,pplev,pplay,zzlay, & 2 8 mu0,ptimestep,ptime,zday,pq,pdq,pdteuv) … … 7 13 igcm_n, igcm_no, igcm_no2, igcm_n2d, mmol 8 14 use conc_mod, only: rnew, cpnew 15 use hrtherm_mod, only: hrtherm 9 16 IMPLICIT NONE 10 17 !======================================================================= … … 31 38 ! ------------------ 32 39 ! 33 #include "callkeys.h"40 include "callkeys.h" 34 41 !----------------------------------------------------------------------- 35 42 ! Input/Output … … 442 449 ! deallocate(rm) 443 450 444 return 445 end 451 end subroutine euvheat 452 453 END MODULE euvheat_mod -
trunk/LMDZ.MARS/libf/aeronomars/hrtherm.F
r2615 r3464 1 MODULE hrtherm_mod 2 3 IMPLICIT NONE 4 5 CONTAINS 6 1 7 c********************************************************************** 2 8 … … 11 17 12 18 use param_v4_h, only: ninter,nabs,jfotsout,fluxtop,freccen 19 use jthermcalc_e107_mod, only: jthermcalc_e107 13 20 14 21 implicit none … … 132 139 end do 133 140 134 end 141 end subroutine hrtherm 142 143 END MODULE hrtherm_mod 135 144 -
trunk/LMDZ.MARS/libf/aeronomars/jthermcalc.F
r2674 r3464 1 module jthermcalc_mod 2 3 implicit none 4 5 contains 6 1 7 c********************************************************************** 2 8 … … 16 22 . co2crsc195,co2crsc295,t0, 17 23 . jabsifotsintpar,ninter,nz2 18 24 use jthermcalc_util, only: column, interfast 25 19 26 implicit none 20 27 … … 1027 1034 1028 1035 1029 return 1030 1031 end 1032 1033 1034 1035 c********************************************************************** 1036 c********************************************************************** 1037 1038 subroutine column(ig,nlayer,chemthermod,rm,nesptherm,tx,iz,zenit, 1039 $ co2colx,o2colx,o3pcolx,h2colx,h2ocolx,h2o2colx,o3colx, 1040 $ n2colx,ncolx,nocolx,cocolx,hcolx,no2colx) 1041 1042 c nov 2002 fgg first version 1043 1044 c********************************************************************** 1045 1046 use tracer_mod, only: igcm_o, igcm_co2, igcm_o2, igcm_h2, 1047 & igcm_h2o_vap, igcm_h2o2, igcm_co, igcm_h, 1048 & igcm_o3, igcm_n2, igcm_n, igcm_no, igcm_no2, 1049 & mmol 1050 use param_v4_h, only: radio,gg,masa,kboltzman,n_avog 1051 1052 implicit none 1053 1054 1055 c common variables and constants 1056 include 'callkeys.h' 1057 1058 1059 1060 c local parameters and variables 1061 1062 1063 1064 c input and output variables 1065 1066 integer ig,nlayer 1067 integer chemthermod 1068 integer nesptherm !# of species undergoing chemistry, input 1069 real rm(nlayer,nesptherm) !densities (cm-3), input 1070 real tx(nlayer) !temperature profile, input 1071 real iz(nlayer+1) !height profile, input 1072 real zenit !SZA, input 1073 real co2colx(nlayer) !column density of CO2 (cm^-2), output 1074 real o2colx(nlayer) !column density of O2(cm^-2), output 1075 real o3pcolx(nlayer) !column density of O(3P)(cm^-2), output 1076 real h2colx(nlayer) !H2 column density (cm-2), output 1077 real h2ocolx(nlayer) !H2O column density (cm-2), output 1078 real h2o2colx(nlayer) !column density of H2O2(cm^-2), output 1079 real o3colx(nlayer) !O3 column density (cm-2), output 1080 real n2colx(nlayer) !N2 column density (cm-2), output 1081 real ncolx(nlayer) !N column density (cm-2), output 1082 real nocolx(nlayer) !NO column density (cm-2), output 1083 real cocolx(nlayer) !CO column density (cm-2), output 1084 real hcolx(nlayer) !H column density (cm-2), output 1085 real no2colx(nlayer) !NO2 column density (cm-2), output 1086 1087 1088 c local variables 1089 1090 real xx 1091 real grav(nlayer) 1092 real Hco2,Ho3p,Ho2,Hh2,Hh2o,Hh2o2 1093 real Ho3,Hn2,Hn,Hno,Hco,Hh,Hno2 1094 1095 real co2x(nlayer) 1096 real o2x(nlayer) 1097 real o3px(nlayer) 1098 real cox(nlayer) 1099 real hx(nlayer) 1100 real h2x(nlayer) 1101 real h2ox(nlayer) 1102 real h2o2x(nlayer) 1103 real o3x(nlayer) 1104 real n2x(nlayer) 1105 real nx(nlayer) 1106 real nox(nlayer) 1107 real no2x(nlayer) 1108 1109 integer i,j,k,icol,indexint !indexes 1110 1111 c variables for optical path calculation 1112 1113 integer nz3 1114 ! parameter (nz3=nz*2) 1115 1116 integer jj 1117 real*8 esp(nlayer*2) 1118 real*8 ilayesp(nlayer*2) 1119 real*8 szalayesp(nlayer*2) 1120 integer nlayesp 1121 real*8 zmini 1122 real*8 depth 1123 real*8 espco2, espo2, espo3p, esph2, esph2o, esph2o2,espo3 1124 real*8 espn2,espn,espno,espco,esph,espno2 1125 real*8 rcmnz, rcmmini 1126 real*8 szadeg 1127 1128 1129 ! Tracer indexes in the thermospheric chemistry: 1130 !!! ATTENTION. These values have to be identical to those in chemthermos.F90 1131 !!! If the values are changed there, the same has to be done here !!! 1132 integer,parameter :: i_co2 = 1 1133 integer,parameter :: i_co = 2 1134 integer,parameter :: i_o = 3 1135 integer,parameter :: i_o1d = 4 1136 integer,parameter :: i_o2 = 5 1137 integer,parameter :: i_o3 = 6 1138 integer,parameter :: i_h = 7 1139 integer,parameter :: i_h2 = 8 1140 integer,parameter :: i_oh = 9 1141 integer,parameter :: i_ho2 = 10 1142 integer,parameter :: i_h2o2 = 11 1143 integer,parameter :: i_h2o = 12 1144 integer,parameter :: i_n = 13 1145 integer,parameter :: i_n2d = 14 1146 integer,parameter :: i_no = 15 1147 integer,parameter :: i_no2 = 16 1148 integer,parameter :: i_n2 = 17 1149 ! integer,parameter :: i_co2=1 1150 ! integer,parameter :: i_o2=2 1151 ! integer,parameter :: i_o=3 1152 ! integer,parameter :: i_co=4 1153 ! integer,parameter :: i_h=5 1154 ! integer,parameter :: i_h2=8 1155 ! integer,parameter :: i_h2o=9 1156 ! integer,parameter :: i_h2o2=10 1157 ! integer,parameter :: i_o3=12 1158 ! integer,parameter :: i_n2=13 1159 ! integer,parameter :: i_n=14 1160 ! integer,parameter :: i_no=15 1161 ! integer,parameter :: i_no2=17 1162 1163 1164 c*************************PROGRAM STARTS******************************* 1165 1166 nz3 = nlayer*2 1167 do i=1,nlayer 1168 xx = ( radio + iz(i) ) * 1.e5 1169 grav(i) = gg * masa /(xx**2) 1170 end do 1171 1172 !Scale heights 1173 xx = kboltzman * tx(nlayer) * n_avog / grav(nlayer) 1174 Ho3p = xx / mmol(igcm_o) 1175 Hco2 = xx / mmol(igcm_co2) 1176 Ho2 = xx / mmol(igcm_o2) 1177 Hh2 = xx / mmol(igcm_h2) 1178 Hh2o = xx / mmol(igcm_h2o_vap) 1179 Hh2o2 = xx / mmol(igcm_h2o2) 1180 Hco = xx / mmol(igcm_co) 1181 Hh = xx / mmol(igcm_h) 1182 !Only if O3 chem. required 1183 if(chemthermod.ge.1) 1184 $ Ho3 = xx / mmol(igcm_o3) 1185 !Only if N or ion chem. 1186 if(chemthermod.ge.2) then 1187 Hn2 = xx / mmol(igcm_n2) 1188 Hn = xx / mmol(igcm_n) 1189 Hno = xx / mmol(igcm_no) 1190 Hno2 = xx / mmol(igcm_no2) 1191 endif 1192 ! first loop in altitude : initialisation 1193 do i=nlayer,1,-1 1194 !Column initialisation 1195 co2colx(i) = 0. 1196 o2colx(i) = 0. 1197 o3pcolx(i) = 0. 1198 h2colx(i) = 0. 1199 h2ocolx(i) = 0. 1200 h2o2colx(i) = 0. 1201 o3colx(i) = 0. 1202 n2colx(i) = 0. 1203 ncolx(i) = 0. 1204 nocolx(i) = 0. 1205 cocolx(i) = 0. 1206 hcolx(i) = 0. 1207 no2colx(i) = 0. 1208 !Densities 1209 co2x(i) = rm(i,i_co2) 1210 o2x(i) = rm(i,i_o2) 1211 o3px(i) = rm(i,i_o) 1212 h2x(i) = rm(i,i_h2) 1213 h2ox(i) = rm(i,i_h2o) 1214 h2o2x(i) = rm(i,i_h2o2) 1215 cox(i) = rm(i,i_co) 1216 hx(i) = rm(i,i_h) 1217 !Only if O3 chem. required 1218 if(chemthermod.ge.1) 1219 $ o3x(i) = rm(i,i_o3) 1220 !Only if Nitrogen of ion chemistry requested 1221 if(chemthermod.ge.2) then 1222 n2x(i) = rm(i,i_n2) 1223 nx(i) = rm(i,i_n) 1224 nox(i) = rm(i,i_no) 1225 no2x(i) = rm(i,i_no2) 1226 endif 1227 enddo 1228 ! second loop in altitude : column calculations 1229 do i=nlayer,1,-1 1230 !Routine to calculate the geometrical length of each layer 1231 call espesor_optico_A(ig,i,nlayer,zenit,iz(i),nz3,iz,esp, 1232 $ ilayesp,szalayesp,nlayesp, zmini) 1233 if(ilayesp(nlayesp).eq.-1) then 1234 co2colx(i)=1.e25 1235 o2colx(i)=1.e25 1236 o3pcolx(i)=1.e25 1237 h2colx(i)=1.e25 1238 h2ocolx(i)=1.e25 1239 h2o2colx(i)=1.e25 1240 o3colx(i)=1.e25 1241 n2colx(i)=1.e25 1242 ncolx(i)=1.e25 1243 nocolx(i)=1.e25 1244 cocolx(i)=1.e25 1245 hcolx(i)=1.e25 1246 no2colx(i)=1.e25 1247 else 1248 rcmnz = ( radio + iz(nlayer) ) * 1.e5 1249 rcmmini = ( radio + zmini ) * 1.e5 1250 !Column calculation taking into account the geometrical depth 1251 !calculated before 1252 do j=1,nlayesp 1253 jj=ilayesp(j) 1254 !Top layer 1255 if(jj.eq.nlayer) then 1256 if(zenit.le.60.) then 1257 o3pcolx(i)=o3pcolx(i)+o3px(nlayer)*Ho3p*esp(j) 1258 $ *1.e-5 1259 co2colx(i)=co2colx(i)+co2x(nlayer)*Hco2*esp(j) 1260 $ *1.e-5 1261 h2o2colx(i)=h2o2colx(i)+ 1262 $ h2o2x(nlayer)*Hh2o2*esp(j)*1.e-5 1263 o2colx(i)=o2colx(i)+o2x(nlayer)*Ho2*esp(j) 1264 $ *1.e-5 1265 h2colx(i)=h2colx(i)+h2x(nlayer)*Hh2*esp(j) 1266 $ *1.e-5 1267 h2ocolx(i)=h2ocolx(i)+h2ox(nlayer)*Hh2o*esp(j) 1268 $ *1.e-5 1269 cocolx(i)=cocolx(i)+cox(nlayer)*Hco*esp(j) 1270 $ *1.e-5 1271 hcolx(i)=hcolx(i)+hx(nlayer)*Hh*esp(j) 1272 $ *1.e-5 1273 !Only if O3 chemistry required 1274 if(chemthermod.ge.1) o3colx(i)= 1275 $ o3colx(i)+o3x(nlayer)*Ho3*esp(j) 1276 $ *1.e-5 1277 !Only if N or ion chemistry requested 1278 if(chemthermod.ge.2) then 1279 n2colx(i)=n2colx(i)+n2x(nlayer)*Hn2*esp(j) 1280 $ *1.e-5 1281 ncolx(i)=ncolx(i)+nx(nlayer)*Hn*esp(j) 1282 $ *1.e-5 1283 nocolx(i)=nocolx(i)+nox(nlayer)*Hno*esp(j) 1284 $ *1.e-5 1285 no2colx(i)=no2colx(i)+no2x(nlayer)*Hno2*esp(j) 1286 $ *1.e-5 1287 endif 1288 else if(zenit.gt.60.) then 1289 espco2 =sqrt((rcmnz+Hco2)**2 -rcmmini**2) - esp(j) 1290 espo2 = sqrt((rcmnz+Ho2)**2 -rcmmini**2) - esp(j) 1291 espo3p = sqrt((rcmnz+Ho3p)**2 -rcmmini**2)- esp(j) 1292 esph2 = sqrt((rcmnz+Hh2)**2 -rcmmini**2) - esp(j) 1293 esph2o = sqrt((rcmnz+Hh2o)**2 -rcmmini**2)- esp(j) 1294 esph2o2= sqrt((rcmnz+Hh2o2)**2-rcmmini**2)- esp(j) 1295 espco = sqrt((rcmnz+Hco)**2 -rcmmini**2) - esp(j) 1296 esph = sqrt((rcmnz+Hh)**2 -rcmmini**2) - esp(j) 1297 !Only if O3 chemistry required 1298 if(chemthermod.ge.1) 1299 $ espo3=sqrt((rcmnz+Ho3)**2-rcmmini**2)-esp(j) 1300 !Only if N or ion chemistry requested 1301 if(chemthermod.ge.2) then 1302 espn2 =sqrt((rcmnz+Hn2)**2-rcmmini**2)-esp(j) 1303 espn =sqrt((rcmnz+Hn)**2-rcmmini**2) - esp(j) 1304 espno =sqrt((rcmnz+Hno)**2-rcmmini**2) - esp(j) 1305 espno2=sqrt((rcmnz+Hno2)**2-rcmmini**2)- esp(j) 1306 endif 1307 co2colx(i) = co2colx(i) + espco2*co2x(nlayer) 1308 o2colx(i) = o2colx(i) + espo2*o2x(nlayer) 1309 o3pcolx(i) = o3pcolx(i) + espo3p*o3px(nlayer) 1310 h2colx(i) = h2colx(i) + esph2*h2x(nlayer) 1311 h2ocolx(i) = h2ocolx(i) + esph2o*h2ox(nlayer) 1312 h2o2colx(i)= h2o2colx(i)+ esph2o2*h2o2x(nlayer) 1313 cocolx(i) = cocolx(i) + espco*cox(nlayer) 1314 hcolx(i) = hcolx(i) + esph*hx(nlayer) 1315 !Only if O3 chemistry required 1316 if(chemthermod.ge.1) 1317 $ o3colx(i) = o3colx(i) + espo3*o3x(nlayer) 1318 !Only if N or ion chemistry requested 1319 if(chemthermod.ge.2) then 1320 n2colx(i) = n2colx(i) + espn2*n2x(nlayer) 1321 ncolx(i) = ncolx(i) + espn*nx(nlayer) 1322 nocolx(i) = nocolx(i) + espno*nox(nlayer) 1323 no2colx(i) = no2colx(i) + espno2*no2x(nlayer) 1324 endif 1325 endif !Of if zenit.lt.60 1326 !Other layers 1327 else 1328 co2colx(i) = co2colx(i) + 1329 $ esp(j) * (co2x(jj)+co2x(jj+1)) / 2. 1330 o2colx(i) = o2colx(i) + 1331 $ esp(j) * (o2x(jj)+o2x(jj+1)) / 2. 1332 o3pcolx(i) = o3pcolx(i) + 1333 $ esp(j) * (o3px(jj)+o3px(jj+1)) / 2. 1334 h2colx(i) = h2colx(i) + 1335 $ esp(j) * (h2x(jj)+h2x(jj+1)) / 2. 1336 h2ocolx(i) = h2ocolx(i) + 1337 $ esp(j) * (h2ox(jj)+h2ox(jj+1)) / 2. 1338 h2o2colx(i) = h2o2colx(i) + 1339 $ esp(j) * (h2o2x(jj)+h2o2x(jj+1)) / 2. 1340 cocolx(i) = cocolx(i) + 1341 $ esp(j) * (cox(jj)+cox(jj+1)) / 2. 1342 hcolx(i) = hcolx(i) + 1343 $ esp(j) * (hx(jj)+hx(jj+1)) / 2. 1344 !Only if O3 chemistry required 1345 if(chemthermod.ge.1) 1346 $ o3colx(i) = o3colx(i) + 1347 $ esp(j) * (o3x(jj)+o3x(jj+1)) / 2. 1348 !Only if N or ion chemistry requested 1349 if(chemthermod.ge.2) then 1350 n2colx(i) = n2colx(i) + 1351 $ esp(j) * (n2x(jj)+n2x(jj+1)) / 2. 1352 ncolx(i) = ncolx(i) + 1353 $ esp(j) * (nx(jj)+nx(jj+1)) / 2. 1354 nocolx(i) = nocolx(i) + 1355 $ esp(j) * (nox(jj)+nox(jj+1)) / 2. 1356 no2colx(i) = no2colx(i) + 1357 $ esp(j) * (no2x(jj)+no2x(jj+1)) / 2. 1358 endif 1359 endif !Of if jj.eq.nlayer 1360 end do !Of do j=1,nlayesp 1361 endif !Of ilayesp(nlayesp).eq.-1 1362 enddo !Of do i=nlayer,1,-1 1363 1364 return 1365 1366 1367 end 1368 1369 1370 c********************************************************************** 1371 c********************************************************************** 1372 1373 subroutine interfast(wm,wp,nm,p,nlayer,pin,nl,limdown,limup) 1374 C 1375 C subroutine to perform linear interpolation in pressure from 1D profile 1376 C escin(nl) sampled on pressure grid pin(nl) to profile 1377 C escout(nlayer) on pressure grid p(nlayer). 1378 C 1379 real*8,intent(out) :: wm(nlayer),wp(nlayer) ! interpolation weights 1380 integer,intent(out) :: nm(nlayer) ! index of nearest point 1381 real*8,intent(in) :: pin(nl),p(nlayer) 1382 real*8,intent(in) :: limup,limdown 1383 integer,intent(in) :: nl,nlayer 1384 integer :: n1,n,np,nini 1385 nini=1 1386 do n1=1,nlayer 1387 if(p(n1) .gt. limup .or. p(n1) .lt. limdown) then 1388 wm(n1) = 0.d0 1389 wp(n1) = 0.d0 1390 else 1391 do n = nini,nl-1 1392 if (p(n1).ge.pin(n).and.p(n1).le.pin(n+1)) then 1393 nm(n1)=n 1394 np=n+1 1395 wm(n1)=abs(pin(n)-p(n1))/(pin(np)-pin(n)) 1396 wp(n1)=1.d0 - wm(n1) 1397 nini = n 1398 exit 1399 endif 1400 enddo 1401 endif 1402 enddo 1403 1404 end 1405 1406 1407 c********************************************************************** 1408 c********************************************************************** 1409 1410 subroutine espesor_optico_A (ig,capa,nlayer, szadeg,z, 1411 @ nz3,iz,esp,ilayesp,szalayesp,nlayesp, zmini) 1412 1413 c fgg nov 03 Adaptation to Martian model 1414 c malv jul 03 Corrected z grid. Split in alt & frec codes 1415 c fgg feb 03 first version 1416 ************************************************************************* 1417 1418 use param_v4_h, only: radio 1419 implicit none 1420 1421 c arguments 1422 1423 real szadeg ! I. SZA [rad] 1424 real z ! I. altitude of interest [km] 1425 integer nz3,ig,nlayer ! I. dimension of esp, ylayesp, etc... 1426 ! (=2*nlayer= max# of layers in ray path) 1427 real iz(nlayer+1) ! I. Altitude of each layer 1428 real*8 esp(nz3) ! O. layer widths after geometrically 1429 ! amplified; in [cm] except at TOA 1430 ! where an auxiliary value is used 1431 real*8 ilayesp(nz3) ! O. Indexes of layers along ray path 1432 real*8 szalayesp(nz3) ! O. SZA [deg] " " " 1433 integer nlayesp 1434 ! real*8 nlayesp ! O. # layers along ray path at this z 1435 real*8 zmini ! O. Minimum altitud of ray path [km] 1436 1437 1438 c local variables and constants 1439 1440 integer j,i,capa 1441 integer jmin ! index of min.altitude along ray path 1442 real*8 szarad ! SZA [deg] 1443 real*8 zz 1444 real*8 diz(nlayer+1) 1445 real*8 rkmnz ! distance TOA to center of Planet [km] 1446 real*8 rkmmini ! distance zmini to center of P [km] 1447 real*8 rkmj ! intermediate distance to C of P [km] 1448 c external function 1449 external grid_R8 ! Returns index of layer containing the altitude 1450 ! of interest, z; for example, if 1451 ! zkm(i)=z or zkm(i)<z<zkm(i+1) => grid(z)=i 1452 integer grid_R8 1453 1454 ************************************************************************* 1455 szarad = dble(szadeg)*3.141592d0/180.d0 1456 zz=dble(z) 1457 do i=1,nlayer 1458 diz(i)=dble(iz(i)) 1459 enddo 1460 do j=1,nz3 1461 esp(j) = 0.d0 1462 szalayesp(j) = 777.d0 1463 ilayesp(j) = 0 1464 enddo 1465 nlayesp = 0 1466 1467 ! First case: szadeg<60 1468 ! The optical thickness will be given by 1/cos(sza) 1469 ! We deal with 2 different regions: 1470 ! 1: First, all layers between z and ztop ("upper part of ray") 1471 ! 2: Second, the layer at ztop 1472 if(szadeg.lt.60.d0) then 1473 1474 zmini = zz 1475 if(abs(zz-diz(nlayer)).lt.1.d-3) goto 1357 1476 ! 1st Zone: Upper part of ray 1477 ! 1478 do j=grid_R8(zz,diz,nlayer),nlayer-1 1479 nlayesp = nlayesp + 1 1480 ilayesp(nlayesp) = j 1481 esp(nlayesp) = (diz(j+1)-diz(j)) / cos(szarad) ! [km] 1482 esp(nlayesp) = esp(nlayesp) * 1.d5 ! [cm] 1483 szalayesp(nlayesp) = szadeg 1484 end do 1485 1486 ! 1487 ! 2nd Zone: Top layer 1488 1357 continue 1489 nlayesp = nlayesp + 1 1490 ilayesp(nlayesp) = nlayer 1491 esp(nlayesp) = 1.d0 / cos(szarad) ! aux. non-dimens. factor 1492 szalayesp(nlayesp) = szadeg 1493 1494 1495 ! Second case: 60 < szadeg < 90 1496 ! The optical thickness is evaluated. 1497 ! (the magnitude of the effect of not using cos(sza) is 3.e-5 1498 ! for z=60km & sza=30, and 5e-4 for z=60km & sza=60, approximately) 1499 ! We deal with 2 different regions: 1500 ! 1: First, all layers between z and ztop ("upper part of ray") 1501 ! 2: Second, the layer at ztop ("uppermost layer") 1502 else if(szadeg.le.90.d0.and.szadeg.ge.60.d0) then 1503 1504 zmini=(radio+zz)*sin(szarad)-radio 1505 rkmmini = radio + zmini 1506 1507 if(abs(zz-diz(nlayer)).lt.1.d-4) goto 1470 1508 1509 ! 1st Zone: Upper part of ray 1510 ! 1511 do j=grid_R8(zz,diz,nlayer),nlayer-1 1512 nlayesp = nlayesp + 1 1513 ilayesp(nlayesp) = j 1514 esp(nlayesp) = 1515 & sqrt( (radio+diz(j+1))**2 - rkmmini**2 ) - 1516 & sqrt( (radio+diz(j))**2 - rkmmini**2 ) ! [km] 1517 esp(nlayesp) = esp(nlayesp) * 1.d5 ! [cm] 1518 rkmj = radio+diz(j) 1519 szalayesp(nlayesp) = asin( rkmmini/rkmj ) ! [rad] 1520 szalayesp(nlayesp) = szalayesp(nlayesp) * 180.d0/3.141592 ! [deg] 1521 end do 1522 1470 continue 1523 ! 2nd Zone: Uppermost layer of ray. 1524 ! 1525 nlayesp = nlayesp + 1 1526 ilayesp(nlayesp) = nlayer 1527 rkmnz = radio+diz(nlayer) 1528 esp(nlayesp) = sqrt( rkmnz**2 - rkmmini**2 ) ! aux.factor[km] 1529 esp(nlayesp) = esp(nlayesp) * 1.d5 ! aux.f. [cm] 1530 szalayesp(nlayesp) = asin( rkmmini/rkmnz ) ! [rad] 1531 szalayesp(nlayesp) = szalayesp(nlayesp) * 180.d0/3.141592! [deg] 1532 1533 1534 ! Third case: szadeg > 90 1535 ! The optical thickness is evaluated. 1536 ! We deal with 5 different regions: 1537 ! 1: all layers between z and ztop ("upper part of ray") 1538 ! 2: the layer at ztop ("uppermost layer") 1539 ! 3: the lowest layer, at zmini 1540 ! 4: the layers increasing from zmini to z (here SZA<90) 1541 ! 5: the layers decreasing from z to zmini (here SZA>90) 1542 else if(szadeg.gt.90.d0) then 1543 1544 zmini=(radio+zz)*sin(szarad)-radio 1545 !zmini should be lower than zz, as SZA<90. However, in situations 1546 !where SZA is very close to 90, rounding errors can make zmini 1547 !slightly higher than zz, causing problems in the determination 1548 !of the jmin index. A correction is implemented in the determination 1549 !of jmin, some lines below 1550 rkmmini = radio + zmini 1551 1552 if(zmini.lt.diz(1)) then ! Can see the sun? No => esp(j)=inft 1553 nlayesp = nlayesp + 1 1554 ilayesp(nlayesp) = - 1 ! Value to mark "no sun on view" 1555 ! esp(nlayesp) = 1.e30 1556 1557 else 1558 jmin=grid_R8(zmini,diz,nlayer)+1 1559 !Correction for possible rounding errors when SZA very close 1560 !to 90 degrees 1561 if(jmin.gt.grid_R8(zz,diz,nlayer)) then 1562 write(*,*)'jthermcalc warning: possible rounding error' 1563 write(*,*)'point,sza,layer:',ig,szadeg,capa 1564 jmin=grid_R8(zz,diz,nlayer) 1565 endif 1566 1567 if(abs(zz-diz(nlayer)).lt.1.d-4) goto 9876 1568 1569 ! 1st Zone: Upper part of ray 1570 ! 1571 do j=grid_R8(zz,diz,nlayer),nlayer-1 1572 nlayesp = nlayesp + 1 1573 ilayesp(nlayesp) = j 1574 esp(nlayesp) = 1575 $ sqrt( (radio+diz(j+1))**2 - rkmmini**2 ) - 1576 $ sqrt( (radio+diz(j))**2 - rkmmini**2 ) ! [km] 1577 esp(nlayesp) = esp(nlayesp) * 1.d5 ! [cm] 1578 rkmj = radio+diz(j) 1579 szalayesp(nlayesp) = asin( rkmmini/rkmj ) ! [rad] 1580 szalayesp(nlayesp) = szalayesp(nlayesp) *180.d0/3.141592 ! [deg] 1581 end do 1582 1583 9876 continue 1584 ! 2nd Zone: Uppermost layer of ray. 1585 ! 1586 nlayesp = nlayesp + 1 1587 ilayesp(nlayesp) = nlayer 1588 rkmnz = radio+diz(nlayer) 1589 esp(nlayesp) = sqrt( rkmnz**2 - rkmmini**2 ) !aux.factor[km] 1590 esp(nlayesp) = esp(nlayesp) * 1.d5 !aux.f.[cm] 1591 szalayesp(nlayesp) = asin( rkmmini/rkmnz ) ! [rad] 1592 szalayesp(nlayesp) = szalayesp(nlayesp) *180.d0/3.141592 ! [deg] 1593 1594 ! 3er Zone: Lowestmost layer of ray 1595 ! 1596 if ( jmin .ge. 2 ) then ! If above the planet's surface 1597 j=jmin-1 1598 nlayesp = nlayesp + 1 1599 ilayesp(nlayesp) = j 1600 esp(nlayesp) = 2. * 1601 $ sqrt( (radio+diz(j+1))**2 -rkmmini**2 ) ! [km] 1602 esp(nlayesp) = esp(nlayesp) * 1.d5 ! [cm] 1603 rkmj = radio+diz(j+1) 1604 szalayesp(nlayesp) = asin( rkmmini/rkmj ) ! [rad] 1605 szalayesp(nlayesp) = szalayesp(nlayesp) *180.d0/3.141592 ! [deg] 1606 endif 1607 1608 ! 4th zone: Lower part of ray, increasing from zmin to z 1609 ! ( layers with SZA < 90 deg ) 1610 do j=jmin,grid_R8(zz,diz,nlayer)-1 1611 nlayesp = nlayesp + 1 1612 ilayesp(nlayesp) = j 1613 esp(nlayesp) = 1614 $ sqrt( (radio+diz(j+1))**2 - rkmmini**2 ) 1615 $ - sqrt( (radio+diz(j))**2 - rkmmini**2 ) ! [km] 1616 esp(nlayesp) = esp(nlayesp) * 1.d5 ! [cm] 1617 rkmj = radio+diz(j) 1618 szalayesp(nlayesp) = asin( rkmmini/rkmj ) ! [rad] 1619 szalayesp(nlayesp) = szalayesp(nlayesp) *180.d0/3.141592 ! [deg] 1620 end do 1621 1622 ! 5th zone: Lower part of ray, decreasing from z to zmin 1623 ! ( layers with SZA > 90 deg ) 1624 do j=grid_R8(zz,diz,nlayer)-1, jmin, -1 1625 nlayesp = nlayesp + 1 1626 ilayesp(nlayesp) = j 1627 esp(nlayesp) = 1628 $ sqrt( (radio+diz(j+1))**2 - rkmmini**2 ) 1629 $ - sqrt( (radio+diz(j))**2 - rkmmini**2 ) ! [km] 1630 esp(nlayesp) = esp(nlayesp) * 1.d5 ! [cm] 1631 rkmj = radio+diz(j) 1632 szalayesp(nlayesp) = 3.141592 - asin( rkmmini/rkmj ) ! [rad] 1633 szalayesp(nlayesp) = szalayesp(nlayesp)*180.d0/3.141592 ! [deg] 1634 end do 1635 1636 end if 1637 1638 end if 1639 1640 1641 return 1642 1643 end 1644 1645 1646 1647 c********************************************************************** 1648 c*********************************************************************** 1649 1650 function grid_R8 (z, zgrid, nz) 1651 1652 c Returns the index where z is located within vector zgrid 1653 c The vector zgrid must be monotonously increasing, otherwise program stops. 1654 c If z is outside zgrid limits, or zgrid dimension is nz<2, the program stops. 1655 c 1656 c FGG Aug-2004 Correct z.lt.zgrid(i) to .le. 1657 c MALV Jul-2003 1658 c*********************************************************************** 1659 1660 implicit none 1661 1662 c Arguments 1663 integer nz 1664 real*8 z 1665 real*8 zgrid(nz) 1666 integer grid_R8 1667 1668 c Local 1669 integer i, nz1, nznew 1670 1671 c*** CODE START 1672 1673 if ( z .lt. zgrid(1) ) then 1674 write (*,*) ' GRID/ z outside bounds of zgrid ' 1675 write (*,*) ' z,zgrid(1),zgrid(nz) =', z,zgrid(1),zgrid(nz) 1676 z = zgrid(1) 1677 write(*,*) 'WARNING: error in grid_r8 (jthermcalc.F)' 1678 write(*,*) 'Please check values of z and zgrid above' 1679 endif 1680 if (z .gt. zgrid(nz) ) then 1681 write (*,*) ' GRID/ z outside bounds of zgrid ' 1682 write (*,*) ' z,zgrid(1),zgrid(nz) =', z,zgrid(1),zgrid(nz) 1683 z = zgrid(nz) 1684 write(*,*) 'WARNING: error in grid_r8 (jthermcalc.F)' 1685 write(*,*) 'Please check values of z and zgrid above' 1686 endif 1687 if ( nz .lt. 2 ) then 1688 write (*,*) ' GRID/ zgrid needs 2 points at least ! ' 1689 stop ' Serious error in GRID.F ' 1690 endif 1691 if ( zgrid(1) .ge. zgrid(nz) ) then 1692 write (*,*) ' GRID/ zgrid must increase with index' 1693 stop ' Serious error in GRID.F ' 1694 endif 1695 1696 nz1 = 1 1697 nznew = nz/2 1698 if ( z .gt. zgrid(nznew) ) then 1699 nz1 = nznew 1700 nznew = nz 1701 endif 1702 do i=nz1+1,nznew 1703 if ( z. eq. zgrid(i) ) then 1704 grid_R8=i 1705 return 1706 elseif ( z .le. zgrid(i) ) then 1707 grid_R8 = i-1 1708 return 1709 endif 1710 enddo 1711 grid_R8 = nz 1712 return 1713 1714 1715 1716 end 1717 1718 1719 1720 !c*************************************************** 1721 !c*************************************************** 1722 1723 subroutine flujo(date) 1724 1725 1726 !c fgg nov 2002 first version 1727 !c*************************************************** 1728 1729 use comsaison_h, only: dist_sol 1730 use param_v4_h, only: ninter, 1731 . fluxtop, ct1, ct2, p1, p2 1732 implicit none 1733 1734 1735 ! common variables and constants 1736 include "callkeys.h" 1737 1738 1739 ! Arguments 1740 1741 real date 1742 1743 1744 ! Local variable and constants 1745 1746 integer i 1747 integer inter 1748 real nada 1749 1750 !c************************************************* 1751 1752 if(date.lt.1985.) date=1985. 1753 if(date.gt.2001.) date=2001. 1036 end subroutine jthermcalc 1754 1037 1755 do i=1,ninter 1756 fluxtop(i)=1. 1757 !Variation of solar flux with 11 years solar cycle 1758 !For more details, see Gonzalez-Galindo et al. 2005 1759 !To be improved in next versions 1760 if(i.le.24.and.solvarmod.eq.0) then 1761 fluxtop(i)=(((ct1(i)+p1(i)*date)/2.) 1762 $ *sin(2.*3.1416/11.*(date-1985.-3.1416)) 1763 $ +(ct2(i)+p2(i)*date)+1.)*fluxtop(i) 1764 end if 1765 fluxtop(i)=fluxtop(i)*(1.52/dist_sol)**2 1766 end do 1767 1768 return 1769 end 1038 end module jthermcalc_mod 1039 1040 1041 -
trunk/LMDZ.MARS/libf/aeronomars/jthermcalc_e107.F
r2808 r3464 1 module jthermcalc_e107_mod 2 3 implicit none 4 5 contains 6 1 7 c********************************************************************** 2 8 … … 19 25 . coefit0,coefit1,coefit2,coefit3,coefit4 20 26 use comsaison_h, only: dist_sol 27 use jthermcalc_util, only: column, interfast 21 28 22 29 implicit none … … 1132 1139 jfotsout(:,:,:)=jfotsout(:,:,:)*(1.52/dist_sol)**2 1133 1140 1134 end 1135 1136 1137 1138 1139 1140 1141 end subroutine jthermcalc_e107 1142 1143 end module jthermcalc_e107_mod 1144 1145 1146 1147 1148 1149 -
trunk/LMDZ.MARS/libf/aeronomars/paramfoto_compact.F
r2615 r3464 1 MODULE paramfoto_compact_mod 2 3 IMPLICIT NONE 4 5 CONTAINS 6 1 7 c********************************************************************** 2 8 … … 639 645 cccccc End altitude loop 640 646 641 return 642 643 644 end 647 end subroutine paramfoto_compact 645 648 646 649 … … 693 696 if ( c_output.lt.1.d-30) c_output=1.d-30 694 697 695 696 return697 698 c END 698 end 699 end subroutine implicito 699 700 700 701 … … 772 773 return 773 774 774 end 775 end function ionsec_nplus 775 776 776 777 … … 840 841 return 841 842 842 end 843 end function ionsec_n2plus 843 844 844 845 … … 915 916 return 916 917 917 end 918 end function ionsec_oplus 918 919 919 920 … … 993 994 return 994 995 995 end 996 end function ionsec_coplus 996 997 997 998 … … 1069 1070 return 1070 1071 1071 end 1072 end function ionsec_co2plus 1072 1073 1073 1074 … … 1139 1140 return 1140 1141 1141 end 1142 end function ionsec_o2plus 1142 1143 1143 1144 … … 1346 1347 1347 1348 1348 return 1349 1350 1351 end 1349 end subroutine phdisrate 1352 1350 1353 1351 … … 1989 1987 endif !Of if(chemthermod.eq.3) 1990 1988 1991 return 1992 end 1989 end subroutine getch 1993 1990 1994 1991 … … 2050 2047 2051 2048 2052 external ionsec_nplus2053 real*8 ionsec_nplus2054 2055 external ionsec_n2plus2056 real*8 ionsec_n2plus2057 2058 external ionsec_oplus2059 real*8 ionsec_oplus2049 ! external ionsec_nplus 2050 ! real*8 ionsec_nplus 2051 2052 ! external ionsec_n2plus 2053 ! real*8 ionsec_n2plus 2054 2055 ! external ionsec_oplus 2056 ! real*8 ionsec_oplus 2060 2057 2061 external ionsec_coplus2062 real*8 ionsec_coplus2063 2064 external ionsec_co2plus2065 real*8 ionsec_co2plus2066 2067 external ionsec_o2plus2068 real*8 ionsec_o2plus2058 ! external ionsec_coplus 2059 ! real*8 ionsec_coplus 2060 2061 ! external ionsec_co2plus 2062 ! real*8 ionsec_co2plus 2063 2064 ! external ionsec_o2plus 2065 ! real*8 ionsec_o2plus 2069 2066 2070 2067 … … 2688 2685 endif !Of chemthermod.eq.3 2689 2686 2690 return2691 2687 c END 2692 end 2688 end subroutine lifetimes 2693 2689 2694 2690 … … 2868 2864 deltat = tminaux 2869 2865 2870 return2871 2866 c END 2872 end 2867 end subroutine timemarching 2873 2868 2874 2869 … … 2925 2920 integer j 2926 2921 2927 external ionsec_nplus 2928 real*8 ionsec_nplus 2929 2930 external ionsec_n2plus 2931 real*8 ionsec_n2plus 2932 2933 external ionsec_oplus 2934 real*8 ionsec_oplus 2935 2936 external ionsec_coplus 2937 real*8 ionsec_coplus 2938 2939 external ionsec_co2plus 2940 real*8 ionsec_co2plus 2941 2942 external ionsec_o2plus 2943 real*8 ionsec_o2plus 2944 2945 logical firstcall 2946 save firstcall 2947 2922 logical,save :: firstcall=.true. 2948 2923 !$OMP THREADPRIVATE(firstcall) 2949 2950 data firstcall /.true./2951 2924 2952 2925 ccccccccccccccc CODE STARTS … … 4100 4073 endif !Of chemthermod.eq.3 4101 4074 4102 return4103 4075 c END 4104 end 4076 end subroutine prodsandlosses 4105 4077 4106 4078 … … 4243 4215 parameter (cociopt=1) 4244 4216 4245 c external functions4246 c4247 external cociente4248 real*8 cociente4249 4250 external ionsec_nplus4251 real*8 ionsec_nplus4252 4253 external ionsec_n2plus4254 real*8 ionsec_n2plus4255 4256 external ionsec_oplus4257 real*8 ionsec_oplus4258 4259 external ionsec_coplus4260 real*8 ionsec_coplus4261 4262 external ionsec_co2plus4263 real*8 ionsec_co2plus4264 4265 external ionsec_o2plus4266 real*8 ionsec_o2plus4267 4268 external avg4269 real*8 avg4270 4271 external dif4272 real*8 dif4273 4217 4274 4218 real*8 log1 … … 5083 5027 5084 5028 5085 return 5086 end 5029 end subroutine EF_oscilacion 5087 5030 5088 5031 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! … … 5093 5036 avg = (log1+log2+log3)*0.333 5094 5037 return 5095 end 5038 end function avg 5096 5039 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 5097 5040 function dif(log1,log2,log3,avg) … … 5104 5047 & abs(log3-avg) ) * 0.333 5105 5048 return 5106 end 5049 end function dif 5107 5050 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 5108 5051 … … 5177 5120 return 5178 5121 c END 5179 end 5122 end function cociente 5123 5124 END MODULE paramfoto_compact_mod -
trunk/LMDZ.MARS/libf/aeronomars/photochemistry.F90
r3462 r3464 1 module photochemistry_mod 2 3 implicit none 4 5 contains 6 1 7 !**************************************************************** 2 8 ! … … 23 29 24 30 use param_v4_h, only: jion 31 use jthermcalc_e107_mod, only: jthermcalc_e107 32 use paramfoto_compact_mod, only: phdisrate 25 33 26 34 implicit none … … 4339 4347 4340 4348 end subroutine photochemistry 4349 4350 end module photochemistry_mod -
trunk/LMDZ.MARS/libf/aeronomars/thermosphere.F
r3443 r3464 18 18 use moldiff_red_mod, only: moldiff_red ! old molecular diffusion scheme 19 19 use moldiff_MPF_mod, only: moldiff_MPF ! new molecular diffusion scheme 20 use euvheat_mod, only: euvheat 20 21 use conc_mod, only: rnew, cpnew 21 22 USE comcstfi_h, only: r, cpp
Note: See TracChangeset
for help on using the changeset viewer.