Changeset 414 for trunk/LMDZ.MARS/libf
- Timestamp:
- Nov 23, 2011, 10:14:21 AM (14 years ago)
- Location:
- trunk/LMDZ.MARS/libf
- Files:
-
- 54 added
- 7 edited
Legend:
- Unmodified
- Added
- Removed
-
TabularUnified trunk/LMDZ.MARS/libf/aeronomars/moldiff.F ¶
r38 r414 32 32 c Local 33 33 c 34 real hco2(ncomptot),ho 34 35 integer,parameter :: ncompmoldiff = 14 36 real hco2(ncompmoldiff),ho 35 37 36 38 integer ig,nz,l,n,nn,iq … … 39 41 real hp(nlayermx) 40 42 real tt(nlayermx) 41 real qq(nlayermx,ncomp tot)43 real qq(nlayermx,ncompmoldiff) 42 44 real dmmeandz(nlayermx) 43 real qnew(nlayermx,ncomp tot)45 real qnew(nlayermx,ncompmoldiff) 44 46 real zlocal(nlayermx) 45 real alf(ncomp tot-1,ncomptot-1)46 real alfinv(nlayermx,ncomp tot-1,ncomptot-1)47 real indx(ncomp tot-1)48 real b(nlayermx,ncomp tot-1)49 real y(ncomp tot-1,ncomptot-1)50 real aa(nlayermx,ncomp tot-1,ncomptot-1)51 real bb(nlayermx,ncomp tot-1,ncomptot-1)52 real cc(nlayermx,ncomp tot-1,ncomptot-1)47 real alf(ncompmoldiff-1,ncompmoldiff-1) 48 real alfinv(nlayermx,ncompmoldiff-1,ncompmoldiff-1) 49 real indx(ncompmoldiff-1) 50 real b(nlayermx,ncompmoldiff-1) 51 real y(ncompmoldiff-1,ncompmoldiff-1) 52 real aa(nlayermx,ncompmoldiff-1,ncompmoldiff-1) 53 real bb(nlayermx,ncompmoldiff-1,ncompmoldiff-1) 54 real cc(nlayermx,ncompmoldiff-1,ncompmoldiff-1) 53 55 real atri(nlayermx-2) 54 56 real btri(nlayermx-2) … … 56 58 real rtri(nlayermx-2) 57 59 real qtri(nlayermx-2) 58 real alfdiag(ncomp tot-1)59 real wi(ncomp tot), flux(ncomptot), pote60 real alfdiag(ncompmoldiff-1) 61 real wi(ncompmoldiff), flux(ncompmoldiff), pote 60 62 61 63 cccccccccccccccccccccccccccccccccccccccccccccccccccccccc … … 95 97 integer,save :: g_h2o=0 96 98 97 integer,save :: gcmind(ncomp tot) ! array of GCM indexes99 integer,save :: gcmind(ncompmoldiff) ! array of GCM indexes 98 100 integer ierr 99 101 100 102 logical,save :: firstcall=.true. 101 real abfac(ncomp tot)102 real,save :: dij(ncomp tot,ncomptot)103 real abfac(ncompmoldiff) 104 real,save :: dij(ncompmoldiff,ncompmoldiff) 103 105 104 106 ! Initializations at first call … … 215 217 & +pdtconduc(ig,l)*ptimestep 216 218 217 do nn=1,ncomp tot219 do nn=1,ncompmoldiff 218 220 qq(l,nn)=pq(ig,l,gcmind(nn))+pdq(ig,l,gcmind(nn))*ptimestep 219 221 qq(l,nn)=max(qq(l,nn),1.e-30) … … 230 232 & +pdtconduc(ig,nz)*ptimestep 231 233 232 do nn=1,ncomp tot234 do nn=1,ncompmoldiff 233 235 qq(1,nn)=pq(ig,1,gcmind(nn))+pdq(ig,1,gcmind(nn))*ptimestep 234 236 qq(nz,nn)=pq(ig,nz,gcmind(nn))+pdq(ig,nz,gcmind(nn))*ptimestep … … 251 253 ntot=pplay(ig,l)/(1.381e-23*tt(l)) ! in #/m3 252 254 253 do nn=1,ncomp tot-1 ! rows255 do nn=1,ncompmoldiff-1 ! rows 254 256 alfdiag(nn)=0. 255 257 dcoef1=dij(nn,i_o)*ptfac 256 do n=1,ncomp tot-1 ! columns258 do n=1,ncompmoldiff-1 ! columns 257 259 y(nn,n)=0. 258 260 dcoef=dij(nn,n)*ptfac … … 277 279 c Inverting the alfa matrix 278 280 c 279 call ludcmp(alf,ncomp tot-1,ncomptot-1,indx,d,ierr)281 call ludcmp(alf,ncompmoldiff-1,ncompmoldiff-1,indx,d,ierr) 280 282 281 283 c TEMPORAIRE ***************************** … … 291 293 end if 292 294 c ******************************************* 293 do n=1,ncomp tot-1294 call lubksb(alf,ncomp tot-1,ncomptot-1,indx,y(1,n))295 do nn=1,ncomp tot-1295 do n=1,ncompmoldiff-1 296 call lubksb(alf,ncompmoldiff-1,ncompmoldiff-1,indx,y(1,n)) 297 do nn=1,ncompmoldiff-1 296 298 alfinv(l,nn,n)=y(nn,n)/hh 297 299 enddo … … 309 311 del1=hp(l)*pplay(ig,l)/(g*ptimestep) 310 312 del2=(hp(l)/2)**2*pplay(ig,l)/(g*ptimestep) 311 do nn=1,ncomp tot-1312 do n=1,ncomp tot-1313 do nn=1,ncompmoldiff-1 314 do n=1,ncompmoldiff-1 313 315 dalfinvdz=(alfinv(l+1,nn,n)-alfinv(l-1,nn,n))/hp(l) 314 316 aa(l,nn,n)=-dalfinvdz/del1+alfinv(l,nn,n)/del2 … … 336 338 c (Hunten 1973, eq. 5) 337 339 338 do n=1,ncomp tot340 do n=1,ncompmoldiff 339 341 wi(n)=1. 340 342 flux(n)=0. … … 379 381 c 380 382 381 do nn=1,ncomp tot-1383 do nn=1,ncompmoldiff-1 382 384 do l=2,nz-1 383 385 atri(l-1)=aa(l,nn,nn) … … 385 387 ctri(l-1)=cc(l,nn,nn) 386 388 rtri(l-1)=qq(l,nn) 387 do n=1,ncomp tot-1389 do n=1,ncompmoldiff-1 388 390 rtri(l-1)=rtri(l-1)-(aa(l,nn,n)*qq(l-1,n) 389 391 & +bb(l,nn,n)*qq(l,n) … … 434 436 if(zlocal(l).gt.65000.)then 435 437 pdqdiff(ig,l,g_o)=0. 436 do n=1,ncomp tot-1438 do n=1,ncompmoldiff-1 437 439 pdqdiff(ig,l,gcmind(n))=(qnew(l,n)-qq(l,n)) 438 440 pdqdiff(ig,l,g_o)=pdqdiff(ig,l,g_o)-(qnew(l,n)-qq(l,n)) -
TabularUnified trunk/LMDZ.MARS/libf/aeronomars/moldiffcoeff.F ¶
r38 r414 22 22 c Input/Output 23 23 c ------------ 24 real dij(ncomptot,ncomptot) 24 integer,parameter :: ncompmoldiff = 14 25 real dij(ncompmoldiff,ncompmoldiff) 25 26 26 27 c Local variables: … … 63 64 integer,save :: g_h2o=0 64 65 65 integer,save :: gcmind(ncomp tot)66 integer,save :: gcmind(ncompmoldiff) 66 67 67 68 real dnh … … 187 188 print*,'moldiffcoef: COEFF CALC' 188 189 open(56,file='coeffs.dat',status='unknown') 189 do n=1,ncomp tot190 do n=1,ncompmoldiff 190 191 if (dij(i_h2,n).gt.0.0) then 191 do nn=n,ncomp tot192 do nn=n,ncompmoldiff 192 193 dij(nn,n)=dij(i_h2,n) 193 194 & *sqrt(mmol(g_h2)/mmol(gcmind(nn))) … … 198 199 if (dij(i_h2,n).eq.0.0) then 199 200 dnh=dij(i_h,i_o)*sqrt(mmol(g_o)/mmol(gcmind(n))) 200 do nn=n,ncomp tot201 do nn=n,ncompmoldiff 201 202 dij(nn,n)=dnh*sqrt(mmol(g_h)/mmol(gcmind(nn))) 202 203 if(n.eq.nn) dij(nn,n)=1.0 … … 206 207 enddo 207 208 208 do n=1,ncomp tot209 do nn=n,ncomp tot209 do n=1,ncompmoldiff 210 do nn=n,ncompmoldiff 210 211 write(56,*) n,nn,dij(n,nn) !*1.e5/1.381e-23/(273**1.75) 211 212 enddo -
TabularUnified trunk/LMDZ.MARS/libf/phymars/callkeys.h ¶
r358 r414 16 16 17 17 COMMON/callkeys_i/iradia,iaervar,iddist,ilwd,ilwb,ilwn,ncouche & 18 & ,dustbin,nqchem_min 18 & ,dustbin,nqchem_min,nltemodel,nircorr 19 19 20 20 COMMON/callkeys_r/topdustref,solarcondate,semi,alphan … … 52 52 logical photochem 53 53 integer nqchem_min 54 integer nltemodel 55 integer nircorr 54 56 55 57 integer swrtype ! type of short wave (solar wavelength) radiative -
TabularUnified trunk/LMDZ.MARS/libf/phymars/inifis.F ¶
r358 r414 134 134 135 135 write(*,*) "Directory where external input files are:" 136 datafile="/ u/forget/WWW/datagcm/datafile"136 datafile="/data/fgglmd/datagcm/datafile" 137 137 call getin("datadir",datafile) ! default path 138 138 write(*,*) " datafile = ",trim(datafile) … … 223 223 write(*,*) " callnlte = ",callnlte 224 224 225 nltemodel=0 !default value 226 write(*,*) "NLTE model?" 227 write(*,*) "0 -> old model, static O" 228 write(*,*) "1 -> old model, dynamic O" 229 write(*,*) "2 -> new model" 230 write(*,*) "(matters only if callnlte=T)" 231 call getin("nltemodel",nltemodel) 232 write(*,*) " nltemodel = ",nltemodel 233 225 234 write(*,*) "call CO2 NIR absorption ?", 226 235 & "(matters only if callrad=T)" … … 228 237 call getin("callnirco2",callnirco2) 229 238 write(*,*) " callnirco2 = ",callnirco2 230 239 240 write(*,*) "New NIR NLTE correction ?", 241 $ "0-> old model (no correction)", 242 $ "1-> new correction", 243 $ "(matters only if callnirco2=T)" 244 call getin("nircorr",nircorr) 245 write(*,*) " nircorr = ",nircorr 246 231 247 write(*,*) "call turbulent vertical diffusion ?" 232 248 calldifv=.true. ! default value -
TabularUnified trunk/LMDZ.MARS/libf/phymars/nirco2abs.F ¶
r38 r414 1 SUBROUTINE nirco2abs(ngrid,nlayer,pplay,dist_sol, 1 SUBROUTINE nirco2abs(ngrid,nlayer,pplay,dist_sol,nq,pq, 2 2 $ mu0,fract,declin,pdtnirco2) 3 3 … … 16 16 c see NLTE correction-factor of Lopez-Valverde et al (1998) 17 17 c Stephen Lewis 2000 18 c 18 c 19 c jul 2011 malv+fgg New corrections for NLTE implemented 19 20 c 08/2002 : correction for bug when running with diurnal=F 20 21 c … … 48 49 #include "callkeys.h" 49 50 #include "comdiurn.h" 50 51 #include "NIRdata.h" 51 52 52 53 c----------------------------------------------------------------------- … … 62 63 c Local variables : 63 64 c ----------------- 64 INTEGER l,ig, n, nstep 65 INTEGER l,ig, n, nstep,i 65 66 REAL co2heat0, zmu(ngridmx) 66 67 … … 84 85 parameter (n_p0=0.0015888279) 85 86 87 c Variables added to implement NLTE correction factor (feb 2011) 88 real pyy(nlayer) 89 real cor1(nlayer),oldoco2(nlayer),alfa2(nlayer) 90 real p2011,cociente1,merge 91 real cor0,oco2gcm 92 integer nq 93 real pq(ngrid,nlayer,nq) 94 86 95 c---------------------------------------------------------------------- 87 96 … … 97 106 do ig=1,ngrid 98 107 zmu(ig)=sqrt(1224.*mu0(ig)*mu0(ig)+1.)/35. 99 enddo 100 do l=1,nlayer 101 do ig=1,ngrid 102 if(fract(ig).gt.0.) pdtnirco2(ig,l)= 108 109 if(nircorr.eq.1) then 110 do l=1,nlayer 111 pyy(l)=pplay(ig,l) 112 enddo 113 114 call interpnir(cor1,pyy,nlayer,corgcm,pres1d,npres) 115 call interpnir(oldoco2,pyy,nlayer,oco21d,pres1d,npres) 116 call interpnir(alfa2,pyy,nlayer,alfa,pres1d,npres) 117 endif 118 119 do l=1,nlayer 120 ! Calculations for the O/CO2 correction 121 if(nircorr.eq.1) then 122 cor0=1./(1.+n_p0/pplay(ig,l))**n_b 123 if(pq(ig,l,1).gt.1.e-6) then 124 oco2gcm=pq(ig,l,3)/pq(ig,l,1) 125 else 126 oco2gcm=1.e6 127 endif 128 cociente1=oco2gcm/oldoco2(l) 129 merge=alog10(cociente1)*alfa2(l)+alog10(cor0)* 130 $ (1.-alfa2(l)) 131 merge=10**merge 132 p2011=sqrt(merge)*cor0 133 else if (nircorr.eq.0) then 134 p2011=1. 135 cor1(l)=1. 136 endif 137 138 if(fract(ig).gt.0.) pdtnirco2(ig,l)= 103 139 & co2heat0*sqrt((700.*zmu(ig))/pplay(ig,l)) 104 140 & /(1.+n_p0/pplay(ig,l))**n_b 105 141 ! Corrections from tabulation 142 $ * cor1(l) * p2011 106 143 c OLD SCHEME (forget et al. 1999) 107 144 c s co2heat0*sqrt((700.*zmu(ig))/pplay(ig,l)) … … 109 146 enddo 110 147 enddo 148 111 149 112 150 c Averaging over diurnal cycle (if diurnal=F) … … 128 166 do ig=1,ngrid 129 167 zmu(ig)=sqrt(1224.*mu0_int(ig)*mu0_int(ig)+1.)/35. 130 enddo 131 do l=1,nlayer 132 do ig=1,ngrid 168 169 if(nircorr.eq.1) then 170 do l=1,nlayer 171 pyy(l)=pplay(ig,l) 172 enddo 173 call interpnir(cor1,pyy,nlayer,corgcm,pres1d,npres) 174 call interpnir(oldoco2,pyy,nlayer,oco21d,pres1d,npres) 175 call interpnir(alfa2,pyy,nlayer,alfa,pres1d,npres) 176 endif 177 178 do l=1,nlayer 179 if(nircorr.eq.1) then 180 cor0=1./(1.+n_p0/pplay(ig,l))**n_b 181 oco2gcm=pq(ig,l,3)/pq(ig,l,1) 182 cociente1=oco2gcm/oldoco2(l) 183 merge=alog10(cociente1)*alfa2(l)+alog10(cor0)* 184 $ (1.-alfa2(l)) 185 merge=10**merge 186 p2011=sqrt(merge)*cor0 187 else if (nircorr.eq.0) then 188 p2011=1. 189 cor1(l)=1. 190 endif 191 133 192 if(fract_int(ig).gt.0.) pdtnirco2(ig,l)= 134 193 & pdtnirco2(ig,l) + (1/float(nstep))* 135 194 & co2heat0*sqrt((700.*zmu(ig))/pplay(ig,l)) 136 195 & /(1.+n_p0/pplay(ig,l))**n_b 196 ! Corrections from tabulation 197 $ * cor1(l) * p2011 137 198 enddo 138 199 enddo … … 143 204 end 144 205 206 207 208 subroutine interpnir(escout,p,nlayer,escin,pin,nl) 209 C 210 C subroutine to perform linear interpolation in pressure from 1D profile 211 C escin(nl) sampled on pressure grid pin(nl) to profile 212 C escout(nlayer) on pressure grid p(nlayer). 213 C 214 real escout(nlayer),p(nlayer) 215 real escin(nl),pin(nl),wm,wp 216 integer nl,nlayer,n1,n,nm,np 217 do n1=1,nlayer 218 if(p(n1) .gt. 1500. .or. p(n1) .lt. 1.0e-13) then 219 escout(n1) = 0.0 220 else 221 do n = 1,nl-1 222 if (p(n1).le.pin(n).and.p(n1).ge.pin(n+1)) then 223 nm=n 224 np=n+1 225 wm=abs(pin(np)-p(n1))/(pin(nm)-pin(np)) 226 wp=1.0 - wm 227 endif 228 enddo 229 escout(n1) = escin(nm)*wm + escin(np)*wp 230 endif 231 enddo 232 return 233 end -
TabularUnified trunk/LMDZ.MARS/libf/phymars/nltecool.F ¶
r38 r414 1 1 c************************************************************************** 2 2 c 3 subroutine nltecool(ngrid,nlayer, pplay,pt,dtnlte)3 subroutine nltecool(ngrid,nlayer,nq,pplay,pt,pq,dtnlte) 4 4 c 5 5 c This code was designed as a delivery for the "Martian Environment Models" … … 25 25 c Version 1d data contained in original input file "nlte_escape.dat" 26 26 c now stored in include file "nltedata.h" Y.Wanherdrick 09/2000 27 28 c jul 2011 fgg Modified to allow variable O 27 29 c 28 30 c*************************************************************************** … … 31 33 32 34 #include "nltedata.h" ! (Equivalent to the reading of the "nlte_escape.dat" file) 33 35 #include "dimensions.h" 36 #include "dimphys.h" 37 #include "chimiedata.h" 38 #include "conc.h" !Added to have "dynamic composition" in the scheme 39 #include "tracer.h" !" 40 #include "callkeys.h" 41 34 42 c Input and output variables 35 43 c 36 44 integer ngrid ! no. of horiz. gridpoints 37 45 integer nlayer ! no. of atmospheric layers 46 integer nq ! no. of tracers 38 47 real pplay(ngrid,nlayer) ! input pressure grid 39 48 real pt(ngrid,nlayer) ! input temperatures 49 real pq(ngrid,nlayer,nq) ! input mmrs 40 50 real dtnlte(ngrid,nlayer) ! output temp. tendencies 41 51 … … 150 160 call interp1(escf2,pyy,nlayer,ef2,pnb,np) 151 161 call interp1(escf1,pyy,nlayer,ef1,pnb,np) 152 call interp3(co2,o3p,n2co,pyy,nlayer, 153 & co2vmr,o3pvmr,n2covmr,pnb,np) 162 if(nltemodel.eq.0) then 163 call interp3(co2,o3p,n2co,pyy,nlayer, 164 & co2vmr,o3pvmr,n2covmr,pnb,np) 165 endif 154 166 155 167 do i=1,nlayer ! loop over layers … … 165 177 c if(pt(j,i).lt.1.0)print*,pt(j,i) 166 178 nt = pyy(i)/(1.381e-17*pt(j,i)) ! nt in cm-3 179 !Dynamic composition 180 if(nltemodel.eq.1) then 181 co2(i)=pq(j,i,1)*mmean(j,i)/mmol(1) 182 o3p(i)=pq(j,i,3)*mmean(j,i)/mmol(3) 183 n2co(i)=pq(j,i,2)*mmean(j,i)/mmol(2) + 184 $ pq(j,i,12)*mmean(j,i)/mmol(12) 185 endif 186 187 !Mixing ratio to density 167 188 co2(i)=co2(i)*nt ! CO2 density in cm-3 168 189 o3p(i)=o3p(i)*nt ! O3p density in cm-3 -
TabularUnified trunk/LMDZ.MARS/libf/phymars/physiq.F ¶
r411 r414 65 65 c Nb: See callradite.F for more information. 66 66 c Mesoscale lines: Aymeric Spiga (2007 - 2011) -- check MESOSCALE flags 67 c jul 2011 malv+fgg: Modified calls to NIR heating routine and 15 um cooling parameterization 67 68 c 68 69 c arguments: … … 323 324 REAL time_phys 324 325 326 ! Added for new NLTE scheme 327 328 real co2vmr_gcm(ngridmx,nlayermx) 329 real n2vmr_gcm(ngridmx,nlayermx) 330 real ovmr_gcm(ngridmx,nlayermx) 331 real covmr_gcm(ngridmx,nlayermx) 332 333 325 334 c Variables for PBL 326 335 … … 422 431 endif 423 432 433 if(callnlte.and.nltemodel.eq.2) call NLTE_leedat 434 if(callnirco2.and.nircorr.eq.1) call NIR_leedat 435 424 436 IF (tracer.AND.water.AND.(ngridmx.NE.1)) THEN 425 437 write(*,*)"physiq: water_param Surface water ice albedo:", … … 428 440 429 441 ENDIF ! (end of "if firstcall") 442 430 443 431 444 c --------------------------------------------------- … … 532 545 c NLTE cooling from CO2 emission 533 546 c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 534 535 IF(callnlte) CALL nltecool(ngrid,nlayer,pplay,pt,zdtnlte) 547 IF(callnlte) then 548 if(nltemodel.eq.0.or.nltemodel.eq.1) then 549 CALL nltecool(ngrid,nlayer,nq,pplay,pt,pq,zdtnlte) 550 else if(nltemodel.eq.2) then 551 do ig=1,ngrid 552 do l=1,nlayer 553 co2vmr_gcm(ig,l)=pq(ig,l,igcm_co2)* 554 $ mmean(ig,l)/mmol(igcm_co2) 555 n2vmr_gcm(ig,l)=pq(ig,l,igcm_n2)* 556 $ mmean(ig,l)/mmol(igcm_n2) 557 covmr_gcm(ig,l)=pq(ig,l,igcm_co)* 558 $ mmean(ig,l)/mmol(igcm_co) 559 ovmr_gcm(ig,l)=pq(ig,l,igcm_o)* 560 $ mmean(ig,l)/mmol(igcm_o) 561 enddo 562 enddo 563 564 CALL NLTEdlvr09_TCOOL(ngrid,nlayer,pplay*9.869e-6, 565 $ pt,zzlay,co2vmr_gcm, n2vmr_gcm, covmr_gcm, 566 $ ovmr_gcm, zdtnlte ) 567 568 do ig=1,ngrid 569 do l=1,nlayer 570 zdtnlte(ig,l)=zdtnlte(ig,l)/86400. 571 enddo 572 enddo 573 endif 574 endif 536 575 537 576 c Find number of layers for LTE radiation calculations … … 605 644 zdtnirco2(:,:)=0 606 645 if (callnirco2) then 607 call nirco2abs (ngrid,nlayer,pplay,dist_sol, 646 call nirco2abs (ngrid,nlayer,pplay,dist_sol,nq,pq, 608 647 . mu0,fract,declin, zdtnirco2) 609 648 endif … … 995 1034 c -------------- 996 1035 IF (photochem .or. thermochem) then 1036 !NB: Photochemistry includes condensation of H2O2 997 1037 PRINT*, 'SURFDUST,SURFICE TO BE IMPLEMENTED. YAAAAAARG.' 998 1038 call calchim(ptimestep,pplay,pplev,pt,pdt,dist_sol,mu0,
Note: See TracChangeset
for help on using the changeset viewer.