Index: trunk/LMDZ.MARS/libf/aeronomars/calchim_mod.F90 =================================================================== --- trunk/LMDZ.MARS/libf/aeronomars/calchim_mod.F90 (revision 3462) +++ trunk/LMDZ.MARS/libf/aeronomars/calchim_mod.F90 (revision 3464) @@ -31,4 +31,6 @@ use conc_mod, only: mmean ! mean molecular mass of the atmosphere use comcstfi_h, only: pi + use chemthermos_mod, only: chemthermos + use photochemistry_mod, only: photochemistry use chemistrydata_mod, only: read_phototable use photolysis_mod, only: init_photolysis, nphot Index: trunk/LMDZ.MARS/libf/aeronomars/chemthermos.F90 =================================================================== --- trunk/LMDZ.MARS/libf/aeronomars/chemthermos.F90 (revision 3462) +++ trunk/LMDZ.MARS/libf/aeronomars/chemthermos.F90 (revision 3464) @@ -1,2 +1,8 @@ + MODULE chemthermos_mod + + IMPLICIT NONE + + CONTAINS + SUBROUTINE chemthermos(ig,nlayer,lswitch,chemthermod,zycol,ztemp, & zdens,zpress,zlocal,zenit,ptimestep,zday,em_no,em_o2) @@ -12,5 +18,7 @@ use param_v4_h, only: Pno, Po2 - USE comcstfi_h + use jthermcalc_e107_mod, only: jthermcalc_e107 + use paramfoto_compact_mod, only: paramfoto_compact + IMPLICIT NONE !======================================================================= @@ -32,5 +40,5 @@ ! ------------------ ! -#include "callkeys.h" + include "callkeys.h" !----------------------------------------------------------------------- ! Input/Output @@ -540,8 +548,8 @@ deallocate(rm) - return - end - - - - + end subroutine chemthermos + + END MODULE chemthermos_mod + + + Index: trunk/LMDZ.MARS/libf/aeronomars/euvheat.F90 =================================================================== --- trunk/LMDZ.MARS/libf/aeronomars/euvheat.F90 (revision 3462) +++ trunk/LMDZ.MARS/libf/aeronomars/euvheat.F90 (revision 3464) @@ -1,2 +1,8 @@ + MODULE euvheat_mod + + IMPLICIT NONE + + CONTAINS + SUBROUTINE euvheat(ngrid,nlayer,nq,pt,pdt,pplev,pplay,zzlay, & mu0,ptimestep,ptime,zday,pq,pdq,pdteuv) @@ -7,4 +13,5 @@ igcm_n, igcm_no, igcm_no2, igcm_n2d, mmol use conc_mod, only: rnew, cpnew + use hrtherm_mod, only: hrtherm IMPLICIT NONE !======================================================================= @@ -31,5 +38,5 @@ ! ------------------ ! -#include "callkeys.h" + include "callkeys.h" !----------------------------------------------------------------------- ! Input/Output @@ -442,4 +449,5 @@ ! deallocate(rm) - return - end + end subroutine euvheat + + END MODULE euvheat_mod Index: trunk/LMDZ.MARS/libf/aeronomars/hrtherm.F =================================================================== --- trunk/LMDZ.MARS/libf/aeronomars/hrtherm.F (revision 3462) +++ trunk/LMDZ.MARS/libf/aeronomars/hrtherm.F (revision 3464) @@ -1,2 +1,8 @@ + MODULE hrtherm_mod + + IMPLICIT NONE + + CONTAINS + c********************************************************************** @@ -11,4 +17,5 @@ use param_v4_h, only: ninter,nabs,jfotsout,fluxtop,freccen + use jthermcalc_e107_mod, only: jthermcalc_e107 implicit none @@ -132,4 +139,6 @@ end do - end + end subroutine hrtherm + + END MODULE hrtherm_mod Index: trunk/LMDZ.MARS/libf/aeronomars/jthermcalc.F =================================================================== --- trunk/LMDZ.MARS/libf/aeronomars/jthermcalc.F (revision 3462) +++ trunk/LMDZ.MARS/libf/aeronomars/jthermcalc.F (revision 3464) @@ -1,2 +1,8 @@ + module jthermcalc_mod + + implicit none + + contains + c********************************************************************** @@ -16,5 +22,6 @@ . co2crsc195,co2crsc295,t0, . jabsifotsintpar,ninter,nz2 - + use jthermcalc_util, only: column, interfast + implicit none @@ -1027,743 +1034,8 @@ - return - - end - - - -c********************************************************************** -c********************************************************************** - - subroutine column(ig,nlayer,chemthermod,rm,nesptherm,tx,iz,zenit, - $ co2colx,o2colx,o3pcolx,h2colx,h2ocolx,h2o2colx,o3colx, - $ n2colx,ncolx,nocolx,cocolx,hcolx,no2colx) - -c nov 2002 fgg first version - -c********************************************************************** - - use tracer_mod, only: igcm_o, igcm_co2, igcm_o2, igcm_h2, - & igcm_h2o_vap, igcm_h2o2, igcm_co, igcm_h, - & igcm_o3, igcm_n2, igcm_n, igcm_no, igcm_no2, - & mmol - use param_v4_h, only: radio,gg,masa,kboltzman,n_avog - - implicit none - - -c common variables and constants - include 'callkeys.h' - - - -c local parameters and variables - - - -c input and output variables - - integer ig,nlayer - integer chemthermod - integer nesptherm !# of species undergoing chemistry, input - real rm(nlayer,nesptherm) !densities (cm-3), input - real tx(nlayer) !temperature profile, input - real iz(nlayer+1) !height profile, input - real zenit !SZA, input - real co2colx(nlayer) !column density of CO2 (cm^-2), output - real o2colx(nlayer) !column density of O2(cm^-2), output - real o3pcolx(nlayer) !column density of O(3P)(cm^-2), output - real h2colx(nlayer) !H2 column density (cm-2), output - real h2ocolx(nlayer) !H2O column density (cm-2), output - real h2o2colx(nlayer) !column density of H2O2(cm^-2), output - real o3colx(nlayer) !O3 column density (cm-2), output - real n2colx(nlayer) !N2 column density (cm-2), output - real ncolx(nlayer) !N column density (cm-2), output - real nocolx(nlayer) !NO column density (cm-2), output - real cocolx(nlayer) !CO column density (cm-2), output - real hcolx(nlayer) !H column density (cm-2), output - real no2colx(nlayer) !NO2 column density (cm-2), output - - -c local variables - - real xx - real grav(nlayer) - real Hco2,Ho3p,Ho2,Hh2,Hh2o,Hh2o2 - real Ho3,Hn2,Hn,Hno,Hco,Hh,Hno2 - - real co2x(nlayer) - real o2x(nlayer) - real o3px(nlayer) - real cox(nlayer) - real hx(nlayer) - real h2x(nlayer) - real h2ox(nlayer) - real h2o2x(nlayer) - real o3x(nlayer) - real n2x(nlayer) - real nx(nlayer) - real nox(nlayer) - real no2x(nlayer) - - integer i,j,k,icol,indexint !indexes - -c variables for optical path calculation - - integer nz3 -! parameter (nz3=nz*2) - - integer jj - real*8 esp(nlayer*2) - real*8 ilayesp(nlayer*2) - real*8 szalayesp(nlayer*2) - integer nlayesp - real*8 zmini - real*8 depth - real*8 espco2, espo2, espo3p, esph2, esph2o, esph2o2,espo3 - real*8 espn2,espn,espno,espco,esph,espno2 - real*8 rcmnz, rcmmini - real*8 szadeg - - - ! Tracer indexes in the thermospheric chemistry: - !!! ATTENTION. These values have to be identical to those in chemthermos.F90 - !!! If the values are changed there, the same has to be done here !!! - integer,parameter :: i_co2 = 1 - integer,parameter :: i_co = 2 - integer,parameter :: i_o = 3 - integer,parameter :: i_o1d = 4 - integer,parameter :: i_o2 = 5 - integer,parameter :: i_o3 = 6 - integer,parameter :: i_h = 7 - integer,parameter :: i_h2 = 8 - integer,parameter :: i_oh = 9 - integer,parameter :: i_ho2 = 10 - integer,parameter :: i_h2o2 = 11 - integer,parameter :: i_h2o = 12 - integer,parameter :: i_n = 13 - integer,parameter :: i_n2d = 14 - integer,parameter :: i_no = 15 - integer,parameter :: i_no2 = 16 - integer,parameter :: i_n2 = 17 -! integer,parameter :: i_co2=1 -! integer,parameter :: i_o2=2 -! integer,parameter :: i_o=3 -! integer,parameter :: i_co=4 -! integer,parameter :: i_h=5 -! integer,parameter :: i_h2=8 -! integer,parameter :: i_h2o=9 -! integer,parameter :: i_h2o2=10 -! integer,parameter :: i_o3=12 -! integer,parameter :: i_n2=13 -! integer,parameter :: i_n=14 -! integer,parameter :: i_no=15 -! integer,parameter :: i_no2=17 - - -c*************************PROGRAM STARTS******************************* - - nz3 = nlayer*2 - do i=1,nlayer - xx = ( radio + iz(i) ) * 1.e5 - grav(i) = gg * masa /(xx**2) - end do - - !Scale heights - xx = kboltzman * tx(nlayer) * n_avog / grav(nlayer) - Ho3p = xx / mmol(igcm_o) - Hco2 = xx / mmol(igcm_co2) - Ho2 = xx / mmol(igcm_o2) - Hh2 = xx / mmol(igcm_h2) - Hh2o = xx / mmol(igcm_h2o_vap) - Hh2o2 = xx / mmol(igcm_h2o2) - Hco = xx / mmol(igcm_co) - Hh = xx / mmol(igcm_h) - !Only if O3 chem. required - if(chemthermod.ge.1) - $ Ho3 = xx / mmol(igcm_o3) - !Only if N or ion chem. - if(chemthermod.ge.2) then - Hn2 = xx / mmol(igcm_n2) - Hn = xx / mmol(igcm_n) - Hno = xx / mmol(igcm_no) - Hno2 = xx / mmol(igcm_no2) - endif - ! first loop in altitude : initialisation - do i=nlayer,1,-1 - !Column initialisation - co2colx(i) = 0. - o2colx(i) = 0. - o3pcolx(i) = 0. - h2colx(i) = 0. - h2ocolx(i) = 0. - h2o2colx(i) = 0. - o3colx(i) = 0. - n2colx(i) = 0. - ncolx(i) = 0. - nocolx(i) = 0. - cocolx(i) = 0. - hcolx(i) = 0. - no2colx(i) = 0. - !Densities - co2x(i) = rm(i,i_co2) - o2x(i) = rm(i,i_o2) - o3px(i) = rm(i,i_o) - h2x(i) = rm(i,i_h2) - h2ox(i) = rm(i,i_h2o) - h2o2x(i) = rm(i,i_h2o2) - cox(i) = rm(i,i_co) - hx(i) = rm(i,i_h) - !Only if O3 chem. required - if(chemthermod.ge.1) - $ o3x(i) = rm(i,i_o3) - !Only if Nitrogen of ion chemistry requested - if(chemthermod.ge.2) then - n2x(i) = rm(i,i_n2) - nx(i) = rm(i,i_n) - nox(i) = rm(i,i_no) - no2x(i) = rm(i,i_no2) - endif - enddo - ! second loop in altitude : column calculations - do i=nlayer,1,-1 - !Routine to calculate the geometrical length of each layer - call espesor_optico_A(ig,i,nlayer,zenit,iz(i),nz3,iz,esp, - $ ilayesp,szalayesp,nlayesp, zmini) - if(ilayesp(nlayesp).eq.-1) then - co2colx(i)=1.e25 - o2colx(i)=1.e25 - o3pcolx(i)=1.e25 - h2colx(i)=1.e25 - h2ocolx(i)=1.e25 - h2o2colx(i)=1.e25 - o3colx(i)=1.e25 - n2colx(i)=1.e25 - ncolx(i)=1.e25 - nocolx(i)=1.e25 - cocolx(i)=1.e25 - hcolx(i)=1.e25 - no2colx(i)=1.e25 - else - rcmnz = ( radio + iz(nlayer) ) * 1.e5 - rcmmini = ( radio + zmini ) * 1.e5 - !Column calculation taking into account the geometrical depth - !calculated before - do j=1,nlayesp - jj=ilayesp(j) - !Top layer - if(jj.eq.nlayer) then - if(zenit.le.60.) then - o3pcolx(i)=o3pcolx(i)+o3px(nlayer)*Ho3p*esp(j) - $ *1.e-5 - co2colx(i)=co2colx(i)+co2x(nlayer)*Hco2*esp(j) - $ *1.e-5 - h2o2colx(i)=h2o2colx(i)+ - $ h2o2x(nlayer)*Hh2o2*esp(j)*1.e-5 - o2colx(i)=o2colx(i)+o2x(nlayer)*Ho2*esp(j) - $ *1.e-5 - h2colx(i)=h2colx(i)+h2x(nlayer)*Hh2*esp(j) - $ *1.e-5 - h2ocolx(i)=h2ocolx(i)+h2ox(nlayer)*Hh2o*esp(j) - $ *1.e-5 - cocolx(i)=cocolx(i)+cox(nlayer)*Hco*esp(j) - $ *1.e-5 - hcolx(i)=hcolx(i)+hx(nlayer)*Hh*esp(j) - $ *1.e-5 - !Only if O3 chemistry required - if(chemthermod.ge.1) o3colx(i)= - $ o3colx(i)+o3x(nlayer)*Ho3*esp(j) - $ *1.e-5 - !Only if N or ion chemistry requested - if(chemthermod.ge.2) then - n2colx(i)=n2colx(i)+n2x(nlayer)*Hn2*esp(j) - $ *1.e-5 - ncolx(i)=ncolx(i)+nx(nlayer)*Hn*esp(j) - $ *1.e-5 - nocolx(i)=nocolx(i)+nox(nlayer)*Hno*esp(j) - $ *1.e-5 - no2colx(i)=no2colx(i)+no2x(nlayer)*Hno2*esp(j) - $ *1.e-5 - endif - else if(zenit.gt.60.) then - espco2 =sqrt((rcmnz+Hco2)**2 -rcmmini**2) - esp(j) - espo2 = sqrt((rcmnz+Ho2)**2 -rcmmini**2) - esp(j) - espo3p = sqrt((rcmnz+Ho3p)**2 -rcmmini**2)- esp(j) - esph2 = sqrt((rcmnz+Hh2)**2 -rcmmini**2) - esp(j) - esph2o = sqrt((rcmnz+Hh2o)**2 -rcmmini**2)- esp(j) - esph2o2= sqrt((rcmnz+Hh2o2)**2-rcmmini**2)- esp(j) - espco = sqrt((rcmnz+Hco)**2 -rcmmini**2) - esp(j) - esph = sqrt((rcmnz+Hh)**2 -rcmmini**2) - esp(j) - !Only if O3 chemistry required - if(chemthermod.ge.1) - $ espo3=sqrt((rcmnz+Ho3)**2-rcmmini**2)-esp(j) - !Only if N or ion chemistry requested - if(chemthermod.ge.2) then - espn2 =sqrt((rcmnz+Hn2)**2-rcmmini**2)-esp(j) - espn =sqrt((rcmnz+Hn)**2-rcmmini**2) - esp(j) - espno =sqrt((rcmnz+Hno)**2-rcmmini**2) - esp(j) - espno2=sqrt((rcmnz+Hno2)**2-rcmmini**2)- esp(j) - endif - co2colx(i) = co2colx(i) + espco2*co2x(nlayer) - o2colx(i) = o2colx(i) + espo2*o2x(nlayer) - o3pcolx(i) = o3pcolx(i) + espo3p*o3px(nlayer) - h2colx(i) = h2colx(i) + esph2*h2x(nlayer) - h2ocolx(i) = h2ocolx(i) + esph2o*h2ox(nlayer) - h2o2colx(i)= h2o2colx(i)+ esph2o2*h2o2x(nlayer) - cocolx(i) = cocolx(i) + espco*cox(nlayer) - hcolx(i) = hcolx(i) + esph*hx(nlayer) - !Only if O3 chemistry required - if(chemthermod.ge.1) - $ o3colx(i) = o3colx(i) + espo3*o3x(nlayer) - !Only if N or ion chemistry requested - if(chemthermod.ge.2) then - n2colx(i) = n2colx(i) + espn2*n2x(nlayer) - ncolx(i) = ncolx(i) + espn*nx(nlayer) - nocolx(i) = nocolx(i) + espno*nox(nlayer) - no2colx(i) = no2colx(i) + espno2*no2x(nlayer) - endif - endif !Of if zenit.lt.60 - !Other layers - else - co2colx(i) = co2colx(i) + - $ esp(j) * (co2x(jj)+co2x(jj+1)) / 2. - o2colx(i) = o2colx(i) + - $ esp(j) * (o2x(jj)+o2x(jj+1)) / 2. - o3pcolx(i) = o3pcolx(i) + - $ esp(j) * (o3px(jj)+o3px(jj+1)) / 2. - h2colx(i) = h2colx(i) + - $ esp(j) * (h2x(jj)+h2x(jj+1)) / 2. - h2ocolx(i) = h2ocolx(i) + - $ esp(j) * (h2ox(jj)+h2ox(jj+1)) / 2. - h2o2colx(i) = h2o2colx(i) + - $ esp(j) * (h2o2x(jj)+h2o2x(jj+1)) / 2. - cocolx(i) = cocolx(i) + - $ esp(j) * (cox(jj)+cox(jj+1)) / 2. - hcolx(i) = hcolx(i) + - $ esp(j) * (hx(jj)+hx(jj+1)) / 2. - !Only if O3 chemistry required - if(chemthermod.ge.1) - $ o3colx(i) = o3colx(i) + - $ esp(j) * (o3x(jj)+o3x(jj+1)) / 2. - !Only if N or ion chemistry requested - if(chemthermod.ge.2) then - n2colx(i) = n2colx(i) + - $ esp(j) * (n2x(jj)+n2x(jj+1)) / 2. - ncolx(i) = ncolx(i) + - $ esp(j) * (nx(jj)+nx(jj+1)) / 2. - nocolx(i) = nocolx(i) + - $ esp(j) * (nox(jj)+nox(jj+1)) / 2. - no2colx(i) = no2colx(i) + - $ esp(j) * (no2x(jj)+no2x(jj+1)) / 2. - endif - endif !Of if jj.eq.nlayer - end do !Of do j=1,nlayesp - endif !Of ilayesp(nlayesp).eq.-1 - enddo !Of do i=nlayer,1,-1 - - return - - - end - - -c********************************************************************** -c********************************************************************** - - subroutine interfast(wm,wp,nm,p,nlayer,pin,nl,limdown,limup) -C -C subroutine to perform linear interpolation in pressure from 1D profile -C escin(nl) sampled on pressure grid pin(nl) to profile -C escout(nlayer) on pressure grid p(nlayer). -C - real*8,intent(out) :: wm(nlayer),wp(nlayer) ! interpolation weights - integer,intent(out) :: nm(nlayer) ! index of nearest point - real*8,intent(in) :: pin(nl),p(nlayer) - real*8,intent(in) :: limup,limdown - integer,intent(in) :: nl,nlayer - integer :: n1,n,np,nini - nini=1 - do n1=1,nlayer - if(p(n1) .gt. limup .or. p(n1) .lt. limdown) then - wm(n1) = 0.d0 - wp(n1) = 0.d0 - else - do n = nini,nl-1 - if (p(n1).ge.pin(n).and.p(n1).le.pin(n+1)) then - nm(n1)=n - np=n+1 - wm(n1)=abs(pin(n)-p(n1))/(pin(np)-pin(n)) - wp(n1)=1.d0 - wm(n1) - nini = n - exit - endif - enddo - endif - enddo - - end - - -c********************************************************************** -c********************************************************************** - - subroutine espesor_optico_A (ig,capa,nlayer, szadeg,z, - @ nz3,iz,esp,ilayesp,szalayesp,nlayesp, zmini) - -c fgg nov 03 Adaptation to Martian model -c malv jul 03 Corrected z grid. Split in alt & frec codes -c fgg feb 03 first version -************************************************************************* - - use param_v4_h, only: radio - implicit none - -c arguments - - real szadeg ! I. SZA [rad] - real z ! I. altitude of interest [km] - integer nz3,ig,nlayer ! I. dimension of esp, ylayesp, etc... - ! (=2*nlayer= max# of layers in ray path) - real iz(nlayer+1) ! I. Altitude of each layer - real*8 esp(nz3) ! O. layer widths after geometrically - ! amplified; in [cm] except at TOA - ! where an auxiliary value is used - real*8 ilayesp(nz3) ! O. Indexes of layers along ray path - real*8 szalayesp(nz3) ! O. SZA [deg] " " " - integer nlayesp -! real*8 nlayesp ! O. # layers along ray path at this z - real*8 zmini ! O. Minimum altitud of ray path [km] - - -c local variables and constants - - integer j,i,capa - integer jmin ! index of min.altitude along ray path - real*8 szarad ! SZA [deg] - real*8 zz - real*8 diz(nlayer+1) - real*8 rkmnz ! distance TOA to center of Planet [km] - real*8 rkmmini ! distance zmini to center of P [km] - real*8 rkmj ! intermediate distance to C of P [km] -c external function - external grid_R8 ! Returns index of layer containing the altitude - ! of interest, z; for example, if - ! zkm(i)=z or zkm(i) grid(z)=i - integer grid_R8 - -************************************************************************* - szarad = dble(szadeg)*3.141592d0/180.d0 - zz=dble(z) - do i=1,nlayer - diz(i)=dble(iz(i)) - enddo - do j=1,nz3 - esp(j) = 0.d0 - szalayesp(j) = 777.d0 - ilayesp(j) = 0 - enddo - nlayesp = 0 - - ! First case: szadeg<60 - ! The optical thickness will be given by 1/cos(sza) - ! We deal with 2 different regions: - ! 1: First, all layers between z and ztop ("upper part of ray") - ! 2: Second, the layer at ztop - if(szadeg.lt.60.d0) then - - zmini = zz - if(abs(zz-diz(nlayer)).lt.1.d-3) goto 1357 - ! 1st Zone: Upper part of ray - ! - do j=grid_R8(zz,diz,nlayer),nlayer-1 - nlayesp = nlayesp + 1 - ilayesp(nlayesp) = j - esp(nlayesp) = (diz(j+1)-diz(j)) / cos(szarad) ! [km] - esp(nlayesp) = esp(nlayesp) * 1.d5 ! [cm] - szalayesp(nlayesp) = szadeg - end do - - ! - ! 2nd Zone: Top layer - 1357 continue - nlayesp = nlayesp + 1 - ilayesp(nlayesp) = nlayer - esp(nlayesp) = 1.d0 / cos(szarad) ! aux. non-dimens. factor - szalayesp(nlayesp) = szadeg - - - ! Second case: 60 < szadeg < 90 - ! The optical thickness is evaluated. - ! (the magnitude of the effect of not using cos(sza) is 3.e-5 - ! for z=60km & sza=30, and 5e-4 for z=60km & sza=60, approximately) - ! We deal with 2 different regions: - ! 1: First, all layers between z and ztop ("upper part of ray") - ! 2: Second, the layer at ztop ("uppermost layer") - else if(szadeg.le.90.d0.and.szadeg.ge.60.d0) then - - zmini=(radio+zz)*sin(szarad)-radio - rkmmini = radio + zmini - - if(abs(zz-diz(nlayer)).lt.1.d-4) goto 1470 - - ! 1st Zone: Upper part of ray - ! - do j=grid_R8(zz,diz,nlayer),nlayer-1 - nlayesp = nlayesp + 1 - ilayesp(nlayesp) = j - esp(nlayesp) = - & sqrt( (radio+diz(j+1))**2 - rkmmini**2 ) - - & sqrt( (radio+diz(j))**2 - rkmmini**2 ) ! [km] - esp(nlayesp) = esp(nlayesp) * 1.d5 ! [cm] - rkmj = radio+diz(j) - szalayesp(nlayesp) = asin( rkmmini/rkmj ) ! [rad] - szalayesp(nlayesp) = szalayesp(nlayesp) * 180.d0/3.141592 ! [deg] - end do - 1470 continue - ! 2nd Zone: Uppermost layer of ray. - ! - nlayesp = nlayesp + 1 - ilayesp(nlayesp) = nlayer - rkmnz = radio+diz(nlayer) - esp(nlayesp) = sqrt( rkmnz**2 - rkmmini**2 ) ! aux.factor[km] - esp(nlayesp) = esp(nlayesp) * 1.d5 ! aux.f. [cm] - szalayesp(nlayesp) = asin( rkmmini/rkmnz ) ! [rad] - szalayesp(nlayesp) = szalayesp(nlayesp) * 180.d0/3.141592! [deg] - - - ! Third case: szadeg > 90 - ! The optical thickness is evaluated. - ! We deal with 5 different regions: - ! 1: all layers between z and ztop ("upper part of ray") - ! 2: the layer at ztop ("uppermost layer") - ! 3: the lowest layer, at zmini - ! 4: the layers increasing from zmini to z (here SZA<90) - ! 5: the layers decreasing from z to zmini (here SZA>90) - else if(szadeg.gt.90.d0) then - - zmini=(radio+zz)*sin(szarad)-radio - !zmini should be lower than zz, as SZA<90. However, in situations - !where SZA is very close to 90, rounding errors can make zmini - !slightly higher than zz, causing problems in the determination - !of the jmin index. A correction is implemented in the determination - !of jmin, some lines below - rkmmini = radio + zmini - - if(zmini.lt.diz(1)) then ! Can see the sun? No => esp(j)=inft - nlayesp = nlayesp + 1 - ilayesp(nlayesp) = - 1 ! Value to mark "no sun on view" -! esp(nlayesp) = 1.e30 - - else - jmin=grid_R8(zmini,diz,nlayer)+1 - !Correction for possible rounding errors when SZA very close - !to 90 degrees - if(jmin.gt.grid_R8(zz,diz,nlayer)) then - write(*,*)'jthermcalc warning: possible rounding error' - write(*,*)'point,sza,layer:',ig,szadeg,capa - jmin=grid_R8(zz,diz,nlayer) - endif - - if(abs(zz-diz(nlayer)).lt.1.d-4) goto 9876 - - ! 1st Zone: Upper part of ray - ! - do j=grid_R8(zz,diz,nlayer),nlayer-1 - nlayesp = nlayesp + 1 - ilayesp(nlayesp) = j - esp(nlayesp) = - $ sqrt( (radio+diz(j+1))**2 - rkmmini**2 ) - - $ sqrt( (radio+diz(j))**2 - rkmmini**2 ) ! [km] - esp(nlayesp) = esp(nlayesp) * 1.d5 ! [cm] - rkmj = radio+diz(j) - szalayesp(nlayesp) = asin( rkmmini/rkmj ) ! [rad] - szalayesp(nlayesp) = szalayesp(nlayesp) *180.d0/3.141592 ! [deg] - end do - - 9876 continue - ! 2nd Zone: Uppermost layer of ray. - ! - nlayesp = nlayesp + 1 - ilayesp(nlayesp) = nlayer - rkmnz = radio+diz(nlayer) - esp(nlayesp) = sqrt( rkmnz**2 - rkmmini**2 ) !aux.factor[km] - esp(nlayesp) = esp(nlayesp) * 1.d5 !aux.f.[cm] - szalayesp(nlayesp) = asin( rkmmini/rkmnz ) ! [rad] - szalayesp(nlayesp) = szalayesp(nlayesp) *180.d0/3.141592 ! [deg] - - ! 3er Zone: Lowestmost layer of ray - ! - if ( jmin .ge. 2 ) then ! If above the planet's surface - j=jmin-1 - nlayesp = nlayesp + 1 - ilayesp(nlayesp) = j - esp(nlayesp) = 2. * - $ sqrt( (radio+diz(j+1))**2 -rkmmini**2 ) ! [km] - esp(nlayesp) = esp(nlayesp) * 1.d5 ! [cm] - rkmj = radio+diz(j+1) - szalayesp(nlayesp) = asin( rkmmini/rkmj ) ! [rad] - szalayesp(nlayesp) = szalayesp(nlayesp) *180.d0/3.141592 ! [deg] - endif - - ! 4th zone: Lower part of ray, increasing from zmin to z - ! ( layers with SZA < 90 deg ) - do j=jmin,grid_R8(zz,diz,nlayer)-1 - nlayesp = nlayesp + 1 - ilayesp(nlayesp) = j - esp(nlayesp) = - $ sqrt( (radio+diz(j+1))**2 - rkmmini**2 ) - $ - sqrt( (radio+diz(j))**2 - rkmmini**2 ) ! [km] - esp(nlayesp) = esp(nlayesp) * 1.d5 ! [cm] - rkmj = radio+diz(j) - szalayesp(nlayesp) = asin( rkmmini/rkmj ) ! [rad] - szalayesp(nlayesp) = szalayesp(nlayesp) *180.d0/3.141592 ! [deg] - end do - - ! 5th zone: Lower part of ray, decreasing from z to zmin - ! ( layers with SZA > 90 deg ) - do j=grid_R8(zz,diz,nlayer)-1, jmin, -1 - nlayesp = nlayesp + 1 - ilayesp(nlayesp) = j - esp(nlayesp) = - $ sqrt( (radio+diz(j+1))**2 - rkmmini**2 ) - $ - sqrt( (radio+diz(j))**2 - rkmmini**2 ) ! [km] - esp(nlayesp) = esp(nlayesp) * 1.d5 ! [cm] - rkmj = radio+diz(j) - szalayesp(nlayesp) = 3.141592 - asin( rkmmini/rkmj ) ! [rad] - szalayesp(nlayesp) = szalayesp(nlayesp)*180.d0/3.141592 ! [deg] - end do - - end if - - end if - - - return - - end - - - -c********************************************************************** -c*********************************************************************** - - function grid_R8 (z, zgrid, nz) - -c Returns the index where z is located within vector zgrid -c The vector zgrid must be monotonously increasing, otherwise program stops. -c If z is outside zgrid limits, or zgrid dimension is nz<2, the program stops. -c -c FGG Aug-2004 Correct z.lt.zgrid(i) to .le. -c MALV Jul-2003 -c*********************************************************************** - - implicit none - -c Arguments - integer nz - real*8 z - real*8 zgrid(nz) - integer grid_R8 - -c Local - integer i, nz1, nznew - -c*** CODE START - - if ( z .lt. zgrid(1) ) then - write (*,*) ' GRID/ z outside bounds of zgrid ' - write (*,*) ' z,zgrid(1),zgrid(nz) =', z,zgrid(1),zgrid(nz) - z = zgrid(1) - write(*,*) 'WARNING: error in grid_r8 (jthermcalc.F)' - write(*,*) 'Please check values of z and zgrid above' - endif - if (z .gt. zgrid(nz) ) then - write (*,*) ' GRID/ z outside bounds of zgrid ' - write (*,*) ' z,zgrid(1),zgrid(nz) =', z,zgrid(1),zgrid(nz) - z = zgrid(nz) - write(*,*) 'WARNING: error in grid_r8 (jthermcalc.F)' - write(*,*) 'Please check values of z and zgrid above' - endif - if ( nz .lt. 2 ) then - write (*,*) ' GRID/ zgrid needs 2 points at least ! ' - stop ' Serious error in GRID.F ' - endif - if ( zgrid(1) .ge. zgrid(nz) ) then - write (*,*) ' GRID/ zgrid must increase with index' - stop ' Serious error in GRID.F ' - endif - - nz1 = 1 - nznew = nz/2 - if ( z .gt. zgrid(nznew) ) then - nz1 = nznew - nznew = nz - endif - do i=nz1+1,nznew - if ( z. eq. zgrid(i) ) then - grid_R8=i - return - elseif ( z .le. zgrid(i) ) then - grid_R8 = i-1 - return - endif - enddo - grid_R8 = nz - return - - - - end - - - -!c*************************************************** -!c*************************************************** - - subroutine flujo(date) - - -!c fgg nov 2002 first version -!c*************************************************** - - use comsaison_h, only: dist_sol - use param_v4_h, only: ninter, - . fluxtop, ct1, ct2, p1, p2 - implicit none - - -! common variables and constants - include "callkeys.h" - - -! Arguments - - real date - - -! Local variable and constants - - integer i - integer inter - real nada - -!c************************************************* - - if(date.lt.1985.) date=1985. - if(date.gt.2001.) date=2001. + end subroutine jthermcalc - do i=1,ninter - fluxtop(i)=1. - !Variation of solar flux with 11 years solar cycle - !For more details, see Gonzalez-Galindo et al. 2005 - !To be improved in next versions - if(i.le.24.and.solvarmod.eq.0) then - fluxtop(i)=(((ct1(i)+p1(i)*date)/2.) - $ *sin(2.*3.1416/11.*(date-1985.-3.1416)) - $ +(ct2(i)+p2(i)*date)+1.)*fluxtop(i) - end if - fluxtop(i)=fluxtop(i)*(1.52/dist_sol)**2 - end do - - return - end + end module jthermcalc_mod + + + Index: trunk/LMDZ.MARS/libf/aeronomars/jthermcalc_e107.F =================================================================== --- trunk/LMDZ.MARS/libf/aeronomars/jthermcalc_e107.F (revision 3462) +++ trunk/LMDZ.MARS/libf/aeronomars/jthermcalc_e107.F (revision 3464) @@ -1,2 +1,8 @@ + module jthermcalc_e107_mod + + implicit none + + contains + c********************************************************************** @@ -19,4 +25,5 @@ . coefit0,coefit1,coefit2,coefit3,coefit4 use comsaison_h, only: dist_sol + use jthermcalc_util, only: column, interfast implicit none @@ -1132,9 +1139,11 @@ jfotsout(:,:,:)=jfotsout(:,:,:)*(1.52/dist_sol)**2 - end - - - - - - + end subroutine jthermcalc_e107 + + end module jthermcalc_e107_mod + + + + + + Index: trunk/LMDZ.MARS/libf/aeronomars/jthermcalc_util.F =================================================================== --- trunk/LMDZ.MARS/libf/aeronomars/jthermcalc_util.F (revision 3464) +++ trunk/LMDZ.MARS/libf/aeronomars/jthermcalc_util.F (revision 3464) @@ -0,0 +1,759 @@ + module jthermcalc_util + + implicit none + + contains + +c********************************************************************** +c********************************************************************** + + subroutine column(ig,nlayer,chemthermod,rm,nesptherm,tx,iz,zenit, + $ co2colx,o2colx,o3pcolx,h2colx,h2ocolx,h2o2colx,o3colx, + $ n2colx,ncolx,nocolx,cocolx,hcolx,no2colx) + +c nov 2002 fgg first version + +c********************************************************************** + + use tracer_mod, only: igcm_o, igcm_co2, igcm_o2, igcm_h2, + & igcm_h2o_vap, igcm_h2o2, igcm_co, igcm_h, + & igcm_o3, igcm_n2, igcm_n, igcm_no, igcm_no2, + & mmol + use param_v4_h, only: radio,gg,masa,kboltzman,n_avog + + implicit none + + +c common variables and constants + include 'callkeys.h' + + + +c local parameters and variables + + + +c input and output variables + + integer ig,nlayer + integer chemthermod + integer nesptherm !# of species undergoing chemistry, input + real rm(nlayer,nesptherm) !densities (cm-3), input + real tx(nlayer) !temperature profile, input + real iz(nlayer+1) !height profile, input + real zenit !SZA, input + real co2colx(nlayer) !column density of CO2 (cm^-2), output + real o2colx(nlayer) !column density of O2(cm^-2), output + real o3pcolx(nlayer) !column density of O(3P)(cm^-2), output + real h2colx(nlayer) !H2 column density (cm-2), output + real h2ocolx(nlayer) !H2O column density (cm-2), output + real h2o2colx(nlayer) !column density of H2O2(cm^-2), output + real o3colx(nlayer) !O3 column density (cm-2), output + real n2colx(nlayer) !N2 column density (cm-2), output + real ncolx(nlayer) !N column density (cm-2), output + real nocolx(nlayer) !NO column density (cm-2), output + real cocolx(nlayer) !CO column density (cm-2), output + real hcolx(nlayer) !H column density (cm-2), output + real no2colx(nlayer) !NO2 column density (cm-2), output + + +c local variables + + real xx + real grav(nlayer) + real Hco2,Ho3p,Ho2,Hh2,Hh2o,Hh2o2 + real Ho3,Hn2,Hn,Hno,Hco,Hh,Hno2 + + real co2x(nlayer) + real o2x(nlayer) + real o3px(nlayer) + real cox(nlayer) + real hx(nlayer) + real h2x(nlayer) + real h2ox(nlayer) + real h2o2x(nlayer) + real o3x(nlayer) + real n2x(nlayer) + real nx(nlayer) + real nox(nlayer) + real no2x(nlayer) + + integer i,j,k,icol,indexint !indexes + +c variables for optical path calculation + + integer nz3 +! parameter (nz3=nz*2) + + integer jj + real*8 esp(nlayer*2) + real*8 ilayesp(nlayer*2) + real*8 szalayesp(nlayer*2) + integer nlayesp + real*8 zmini + real*8 depth + real*8 espco2, espo2, espo3p, esph2, esph2o, esph2o2,espo3 + real*8 espn2,espn,espno,espco,esph,espno2 + real*8 rcmnz, rcmmini + real*8 szadeg + + + ! Tracer indexes in the thermospheric chemistry: + !!! ATTENTION. These values have to be identical to those in chemthermos.F90 + !!! If the values are changed there, the same has to be done here !!! + integer,parameter :: i_co2 = 1 + integer,parameter :: i_co = 2 + integer,parameter :: i_o = 3 + integer,parameter :: i_o1d = 4 + integer,parameter :: i_o2 = 5 + integer,parameter :: i_o3 = 6 + integer,parameter :: i_h = 7 + integer,parameter :: i_h2 = 8 + integer,parameter :: i_oh = 9 + integer,parameter :: i_ho2 = 10 + integer,parameter :: i_h2o2 = 11 + integer,parameter :: i_h2o = 12 + integer,parameter :: i_n = 13 + integer,parameter :: i_n2d = 14 + integer,parameter :: i_no = 15 + integer,parameter :: i_no2 = 16 + integer,parameter :: i_n2 = 17 +! integer,parameter :: i_co2=1 +! integer,parameter :: i_o2=2 +! integer,parameter :: i_o=3 +! integer,parameter :: i_co=4 +! integer,parameter :: i_h=5 +! integer,parameter :: i_h2=8 +! integer,parameter :: i_h2o=9 +! integer,parameter :: i_h2o2=10 +! integer,parameter :: i_o3=12 +! integer,parameter :: i_n2=13 +! integer,parameter :: i_n=14 +! integer,parameter :: i_no=15 +! integer,parameter :: i_no2=17 + + +c*************************PROGRAM STARTS******************************* + + nz3 = nlayer*2 + do i=1,nlayer + xx = ( radio + iz(i) ) * 1.e5 + grav(i) = gg * masa /(xx**2) + end do + + !Scale heights + xx = kboltzman * tx(nlayer) * n_avog / grav(nlayer) + Ho3p = xx / mmol(igcm_o) + Hco2 = xx / mmol(igcm_co2) + Ho2 = xx / mmol(igcm_o2) + Hh2 = xx / mmol(igcm_h2) + Hh2o = xx / mmol(igcm_h2o_vap) + Hh2o2 = xx / mmol(igcm_h2o2) + Hco = xx / mmol(igcm_co) + Hh = xx / mmol(igcm_h) + !Only if O3 chem. required + if(chemthermod.ge.1) + $ Ho3 = xx / mmol(igcm_o3) + !Only if N or ion chem. + if(chemthermod.ge.2) then + Hn2 = xx / mmol(igcm_n2) + Hn = xx / mmol(igcm_n) + Hno = xx / mmol(igcm_no) + Hno2 = xx / mmol(igcm_no2) + endif + ! first loop in altitude : initialisation + do i=nlayer,1,-1 + !Column initialisation + co2colx(i) = 0. + o2colx(i) = 0. + o3pcolx(i) = 0. + h2colx(i) = 0. + h2ocolx(i) = 0. + h2o2colx(i) = 0. + o3colx(i) = 0. + n2colx(i) = 0. + ncolx(i) = 0. + nocolx(i) = 0. + cocolx(i) = 0. + hcolx(i) = 0. + no2colx(i) = 0. + !Densities + co2x(i) = rm(i,i_co2) + o2x(i) = rm(i,i_o2) + o3px(i) = rm(i,i_o) + h2x(i) = rm(i,i_h2) + h2ox(i) = rm(i,i_h2o) + h2o2x(i) = rm(i,i_h2o2) + cox(i) = rm(i,i_co) + hx(i) = rm(i,i_h) + !Only if O3 chem. required + if(chemthermod.ge.1) + $ o3x(i) = rm(i,i_o3) + !Only if Nitrogen of ion chemistry requested + if(chemthermod.ge.2) then + n2x(i) = rm(i,i_n2) + nx(i) = rm(i,i_n) + nox(i) = rm(i,i_no) + no2x(i) = rm(i,i_no2) + endif + enddo + ! second loop in altitude : column calculations + do i=nlayer,1,-1 + !Routine to calculate the geometrical length of each layer + call espesor_optico_A(ig,i,nlayer,zenit,iz(i),nz3,iz,esp, + $ ilayesp,szalayesp,nlayesp, zmini) + if(ilayesp(nlayesp).eq.-1) then + co2colx(i)=1.e25 + o2colx(i)=1.e25 + o3pcolx(i)=1.e25 + h2colx(i)=1.e25 + h2ocolx(i)=1.e25 + h2o2colx(i)=1.e25 + o3colx(i)=1.e25 + n2colx(i)=1.e25 + ncolx(i)=1.e25 + nocolx(i)=1.e25 + cocolx(i)=1.e25 + hcolx(i)=1.e25 + no2colx(i)=1.e25 + else + rcmnz = ( radio + iz(nlayer) ) * 1.e5 + rcmmini = ( radio + zmini ) * 1.e5 + !Column calculation taking into account the geometrical depth + !calculated before + do j=1,nlayesp + jj=ilayesp(j) + !Top layer + if(jj.eq.nlayer) then + if(zenit.le.60.) then + o3pcolx(i)=o3pcolx(i)+o3px(nlayer)*Ho3p*esp(j) + $ *1.e-5 + co2colx(i)=co2colx(i)+co2x(nlayer)*Hco2*esp(j) + $ *1.e-5 + h2o2colx(i)=h2o2colx(i)+ + $ h2o2x(nlayer)*Hh2o2*esp(j)*1.e-5 + o2colx(i)=o2colx(i)+o2x(nlayer)*Ho2*esp(j) + $ *1.e-5 + h2colx(i)=h2colx(i)+h2x(nlayer)*Hh2*esp(j) + $ *1.e-5 + h2ocolx(i)=h2ocolx(i)+h2ox(nlayer)*Hh2o*esp(j) + $ *1.e-5 + cocolx(i)=cocolx(i)+cox(nlayer)*Hco*esp(j) + $ *1.e-5 + hcolx(i)=hcolx(i)+hx(nlayer)*Hh*esp(j) + $ *1.e-5 + !Only if O3 chemistry required + if(chemthermod.ge.1) o3colx(i)= + $ o3colx(i)+o3x(nlayer)*Ho3*esp(j) + $ *1.e-5 + !Only if N or ion chemistry requested + if(chemthermod.ge.2) then + n2colx(i)=n2colx(i)+n2x(nlayer)*Hn2*esp(j) + $ *1.e-5 + ncolx(i)=ncolx(i)+nx(nlayer)*Hn*esp(j) + $ *1.e-5 + nocolx(i)=nocolx(i)+nox(nlayer)*Hno*esp(j) + $ *1.e-5 + no2colx(i)=no2colx(i)+no2x(nlayer)*Hno2*esp(j) + $ *1.e-5 + endif + else if(zenit.gt.60.) then + espco2 =sqrt((rcmnz+Hco2)**2 -rcmmini**2) - esp(j) + espo2 = sqrt((rcmnz+Ho2)**2 -rcmmini**2) - esp(j) + espo3p = sqrt((rcmnz+Ho3p)**2 -rcmmini**2)- esp(j) + esph2 = sqrt((rcmnz+Hh2)**2 -rcmmini**2) - esp(j) + esph2o = sqrt((rcmnz+Hh2o)**2 -rcmmini**2)- esp(j) + esph2o2= sqrt((rcmnz+Hh2o2)**2-rcmmini**2)- esp(j) + espco = sqrt((rcmnz+Hco)**2 -rcmmini**2) - esp(j) + esph = sqrt((rcmnz+Hh)**2 -rcmmini**2) - esp(j) + !Only if O3 chemistry required + if(chemthermod.ge.1) + $ espo3=sqrt((rcmnz+Ho3)**2-rcmmini**2)-esp(j) + !Only if N or ion chemistry requested + if(chemthermod.ge.2) then + espn2 =sqrt((rcmnz+Hn2)**2-rcmmini**2)-esp(j) + espn =sqrt((rcmnz+Hn)**2-rcmmini**2) - esp(j) + espno =sqrt((rcmnz+Hno)**2-rcmmini**2) - esp(j) + espno2=sqrt((rcmnz+Hno2)**2-rcmmini**2)- esp(j) + endif + co2colx(i) = co2colx(i) + espco2*co2x(nlayer) + o2colx(i) = o2colx(i) + espo2*o2x(nlayer) + o3pcolx(i) = o3pcolx(i) + espo3p*o3px(nlayer) + h2colx(i) = h2colx(i) + esph2*h2x(nlayer) + h2ocolx(i) = h2ocolx(i) + esph2o*h2ox(nlayer) + h2o2colx(i)= h2o2colx(i)+ esph2o2*h2o2x(nlayer) + cocolx(i) = cocolx(i) + espco*cox(nlayer) + hcolx(i) = hcolx(i) + esph*hx(nlayer) + !Only if O3 chemistry required + if(chemthermod.ge.1) + $ o3colx(i) = o3colx(i) + espo3*o3x(nlayer) + !Only if N or ion chemistry requested + if(chemthermod.ge.2) then + n2colx(i) = n2colx(i) + espn2*n2x(nlayer) + ncolx(i) = ncolx(i) + espn*nx(nlayer) + nocolx(i) = nocolx(i) + espno*nox(nlayer) + no2colx(i) = no2colx(i) + espno2*no2x(nlayer) + endif + endif !Of if zenit.lt.60 + !Other layers + else + co2colx(i) = co2colx(i) + + $ esp(j) * (co2x(jj)+co2x(jj+1)) / 2. + o2colx(i) = o2colx(i) + + $ esp(j) * (o2x(jj)+o2x(jj+1)) / 2. + o3pcolx(i) = o3pcolx(i) + + $ esp(j) * (o3px(jj)+o3px(jj+1)) / 2. + h2colx(i) = h2colx(i) + + $ esp(j) * (h2x(jj)+h2x(jj+1)) / 2. + h2ocolx(i) = h2ocolx(i) + + $ esp(j) * (h2ox(jj)+h2ox(jj+1)) / 2. + h2o2colx(i) = h2o2colx(i) + + $ esp(j) * (h2o2x(jj)+h2o2x(jj+1)) / 2. + cocolx(i) = cocolx(i) + + $ esp(j) * (cox(jj)+cox(jj+1)) / 2. + hcolx(i) = hcolx(i) + + $ esp(j) * (hx(jj)+hx(jj+1)) / 2. + !Only if O3 chemistry required + if(chemthermod.ge.1) + $ o3colx(i) = o3colx(i) + + $ esp(j) * (o3x(jj)+o3x(jj+1)) / 2. + !Only if N or ion chemistry requested + if(chemthermod.ge.2) then + n2colx(i) = n2colx(i) + + $ esp(j) * (n2x(jj)+n2x(jj+1)) / 2. + ncolx(i) = ncolx(i) + + $ esp(j) * (nx(jj)+nx(jj+1)) / 2. + nocolx(i) = nocolx(i) + + $ esp(j) * (nox(jj)+nox(jj+1)) / 2. + no2colx(i) = no2colx(i) + + $ esp(j) * (no2x(jj)+no2x(jj+1)) / 2. + endif + endif !Of if jj.eq.nlayer + end do !Of do j=1,nlayesp + endif !Of ilayesp(nlayesp).eq.-1 + enddo !Of do i=nlayer,1,-1 + + + end subroutine column + + +c********************************************************************** +c********************************************************************** + + subroutine interfast(wm,wp,nm,p,nlayer,pin,nl,limdown,limup) +C +C subroutine to perform linear interpolation in pressure from 1D profile +C escin(nl) sampled on pressure grid pin(nl) to profile +C escout(nlayer) on pressure grid p(nlayer). +C + real*8,intent(out) :: wm(nlayer),wp(nlayer) ! interpolation weights + integer,intent(out) :: nm(nlayer) ! index of nearest point + real*8,intent(in) :: pin(nl),p(nlayer) + real*8,intent(in) :: limup,limdown + integer,intent(in) :: nl,nlayer + integer :: n1,n,np,nini + logical,parameter :: extra_sanity_checks=.false. + +! Added sanity check: is input p(:) indeed monotonically increasing? + if (extra_sanity_checks) then + do nini=1,nlayer-1 + if (p(nini).gt.p(nini+1)) then + write(*,*) "possible interfast issue, nini=",nini + write(*,*) "p(nini)=",p(nini),"> p(nini+1)=",p(nini+1) + endif + enddo + endif ! of if (extra_sanity_checks) + + nini=1 + do n1=1,nlayer + if(p(n1) .gt. limup .or. p(n1) .lt. limdown) then + wm(n1) = 0.d0 + wp(n1) = 0.d0 + else + do n = nini,nl-1 + if (p(n1).ge.pin(n).and.p(n1).le.pin(n+1)) then + nm(n1)=n + np=n+1 + wm(n1)=abs(pin(n)-p(n1))/(pin(np)-pin(n)) + wp(n1)=1.d0 - wm(n1) + nini = n + exit + endif + enddo + endif + enddo + +! Added sanity check: does nm(:) indeed contain values +! between 1 and nl-1 ? + if (extra_sanity_checks) then + if ((minval(nm)<1).or.(maxval(nm)>nl-1)) then + write(*,*) "interfast issue nm(:) contains incoherent values" + write(*,*) " nm(:) values should be between 1 and ",nl-1 + write(*,*) " but nm(:)=",nm(:) + endif + endif ! of if (extra_sanity_checks) + + end subroutine interfast + + +c********************************************************************** +c********************************************************************** + + subroutine espesor_optico_A (ig,capa,nlayer, szadeg,z, + @ nz3,iz,esp,ilayesp,szalayesp,nlayesp, zmini) + +c fgg nov 03 Adaptation to Martian model +c malv jul 03 Corrected z grid. Split in alt & frec codes +c fgg feb 03 first version +************************************************************************* + + use param_v4_h, only: radio + implicit none + +c arguments + + real szadeg ! I. SZA [rad] + real z ! I. altitude of interest [km] + integer nz3,ig,nlayer ! I. dimension of esp, ylayesp, etc... + ! (=2*nlayer= max# of layers in ray path) + real iz(nlayer+1) ! I. Altitude of each layer + real*8 esp(nz3) ! O. layer widths after geometrically + ! amplified; in [cm] except at TOA + ! where an auxiliary value is used + real*8 ilayesp(nz3) ! O. Indexes of layers along ray path + real*8 szalayesp(nz3) ! O. SZA [deg] " " " + integer nlayesp +! real*8 nlayesp ! O. # layers along ray path at this z + real*8 zmini ! O. Minimum altitud of ray path [km] + + +c local variables and constants + + integer j,i,capa + integer jmin ! index of min.altitude along ray path + real*8 szarad ! SZA [deg] + real*8 zz + real*8 diz(nlayer+1) + real*8 rkmnz ! distance TOA to center of Planet [km] + real*8 rkmmini ! distance zmini to center of P [km] + real*8 rkmj ! intermediate distance to C of P [km] +c external function +c external grid_R8 ! Returns index of layer containing the altitude +c ! of interest, z; for example, if +c ! zkm(i)=z or zkm(i) grid(z)=i +c integer grid_R8 + +************************************************************************* + szarad = dble(szadeg)*3.141592d0/180.d0 + zz=dble(z) + do i=1,nlayer + diz(i)=dble(iz(i)) + enddo + do j=1,nz3 + esp(j) = 0.d0 + szalayesp(j) = 777.d0 + ilayesp(j) = 0 + enddo + nlayesp = 0 + + ! First case: szadeg<60 + ! The optical thickness will be given by 1/cos(sza) + ! We deal with 2 different regions: + ! 1: First, all layers between z and ztop ("upper part of ray") + ! 2: Second, the layer at ztop + if(szadeg.lt.60.d0) then + + zmini = zz + if(abs(zz-diz(nlayer)).lt.1.d-3) goto 1357 + ! 1st Zone: Upper part of ray + ! + do j=grid_R8(zz,diz,nlayer),nlayer-1 + nlayesp = nlayesp + 1 + ilayesp(nlayesp) = j + esp(nlayesp) = (diz(j+1)-diz(j)) / cos(szarad) ! [km] + esp(nlayesp) = esp(nlayesp) * 1.d5 ! [cm] + szalayesp(nlayesp) = szadeg + end do + + ! + ! 2nd Zone: Top layer + 1357 continue + nlayesp = nlayesp + 1 + ilayesp(nlayesp) = nlayer + esp(nlayesp) = 1.d0 / cos(szarad) ! aux. non-dimens. factor + szalayesp(nlayesp) = szadeg + + + ! Second case: 60 < szadeg < 90 + ! The optical thickness is evaluated. + ! (the magnitude of the effect of not using cos(sza) is 3.e-5 + ! for z=60km & sza=30, and 5e-4 for z=60km & sza=60, approximately) + ! We deal with 2 different regions: + ! 1: First, all layers between z and ztop ("upper part of ray") + ! 2: Second, the layer at ztop ("uppermost layer") + else if(szadeg.le.90.d0.and.szadeg.ge.60.d0) then + + zmini=(radio+zz)*sin(szarad)-radio + rkmmini = radio + zmini + + if(abs(zz-diz(nlayer)).lt.1.d-4) goto 1470 + + ! 1st Zone: Upper part of ray + ! + do j=grid_R8(zz,diz,nlayer),nlayer-1 + nlayesp = nlayesp + 1 + ilayesp(nlayesp) = j + esp(nlayesp) = + & sqrt( (radio+diz(j+1))**2 - rkmmini**2 ) - + & sqrt( (radio+diz(j))**2 - rkmmini**2 ) ! [km] + esp(nlayesp) = esp(nlayesp) * 1.d5 ! [cm] + rkmj = radio+diz(j) + szalayesp(nlayesp) = asin( rkmmini/rkmj ) ! [rad] + szalayesp(nlayesp) = szalayesp(nlayesp) * 180.d0/3.141592 ! [deg] + end do + 1470 continue + ! 2nd Zone: Uppermost layer of ray. + ! + nlayesp = nlayesp + 1 + ilayesp(nlayesp) = nlayer + rkmnz = radio+diz(nlayer) + esp(nlayesp) = sqrt( rkmnz**2 - rkmmini**2 ) ! aux.factor[km] + esp(nlayesp) = esp(nlayesp) * 1.d5 ! aux.f. [cm] + szalayesp(nlayesp) = asin( rkmmini/rkmnz ) ! [rad] + szalayesp(nlayesp) = szalayesp(nlayesp) * 180.d0/3.141592! [deg] + + + ! Third case: szadeg > 90 + ! The optical thickness is evaluated. + ! We deal with 5 different regions: + ! 1: all layers between z and ztop ("upper part of ray") + ! 2: the layer at ztop ("uppermost layer") + ! 3: the lowest layer, at zmini + ! 4: the layers increasing from zmini to z (here SZA<90) + ! 5: the layers decreasing from z to zmini (here SZA>90) + else if(szadeg.gt.90.d0) then + + zmini=(radio+zz)*sin(szarad)-radio + !zmini should be lower than zz, as SZA<90. However, in situations + !where SZA is very close to 90, rounding errors can make zmini + !slightly higher than zz, causing problems in the determination + !of the jmin index. A correction is implemented in the determination + !of jmin, some lines below + rkmmini = radio + zmini + + if(zmini.lt.diz(1)) then ! Can see the sun? No => esp(j)=inft + nlayesp = nlayesp + 1 + ilayesp(nlayesp) = - 1 ! Value to mark "no sun on view" +! esp(nlayesp) = 1.e30 + + else + jmin=grid_R8(zmini,diz,nlayer)+1 + !Correction for possible rounding errors when SZA very close + !to 90 degrees + if(jmin.gt.grid_R8(zz,diz,nlayer)) then + write(*,*)'jthermcalc warning: possible rounding error' + write(*,*)'point,sza,layer:',ig,szadeg,capa + jmin=grid_R8(zz,diz,nlayer) + endif + + if(abs(zz-diz(nlayer)).lt.1.d-4) goto 9876 + + ! 1st Zone: Upper part of ray + ! + do j=grid_R8(zz,diz,nlayer),nlayer-1 + nlayesp = nlayesp + 1 + ilayesp(nlayesp) = j + esp(nlayesp) = + $ sqrt( (radio+diz(j+1))**2 - rkmmini**2 ) - + $ sqrt( (radio+diz(j))**2 - rkmmini**2 ) ! [km] + esp(nlayesp) = esp(nlayesp) * 1.d5 ! [cm] + rkmj = radio+diz(j) + szalayesp(nlayesp) = asin( rkmmini/rkmj ) ! [rad] + szalayesp(nlayesp) = szalayesp(nlayesp) *180.d0/3.141592 ! [deg] + end do + + 9876 continue + ! 2nd Zone: Uppermost layer of ray. + ! + nlayesp = nlayesp + 1 + ilayesp(nlayesp) = nlayer + rkmnz = radio+diz(nlayer) + esp(nlayesp) = sqrt( rkmnz**2 - rkmmini**2 ) !aux.factor[km] + esp(nlayesp) = esp(nlayesp) * 1.d5 !aux.f.[cm] + szalayesp(nlayesp) = asin( rkmmini/rkmnz ) ! [rad] + szalayesp(nlayesp) = szalayesp(nlayesp) *180.d0/3.141592 ! [deg] + + ! 3er Zone: Lowestmost layer of ray + ! + if ( jmin .ge. 2 ) then ! If above the planet's surface + j=jmin-1 + nlayesp = nlayesp + 1 + ilayesp(nlayesp) = j + esp(nlayesp) = 2. * + $ sqrt( (radio+diz(j+1))**2 -rkmmini**2 ) ! [km] + esp(nlayesp) = esp(nlayesp) * 1.d5 ! [cm] + rkmj = radio+diz(j+1) + szalayesp(nlayesp) = asin( rkmmini/rkmj ) ! [rad] + szalayesp(nlayesp) = szalayesp(nlayesp) *180.d0/3.141592 ! [deg] + endif + + ! 4th zone: Lower part of ray, increasing from zmin to z + ! ( layers with SZA < 90 deg ) + do j=jmin,grid_R8(zz,diz,nlayer)-1 + nlayesp = nlayesp + 1 + ilayesp(nlayesp) = j + esp(nlayesp) = + $ sqrt( (radio+diz(j+1))**2 - rkmmini**2 ) + $ - sqrt( (radio+diz(j))**2 - rkmmini**2 ) ! [km] + esp(nlayesp) = esp(nlayesp) * 1.d5 ! [cm] + rkmj = radio+diz(j) + szalayesp(nlayesp) = asin( rkmmini/rkmj ) ! [rad] + szalayesp(nlayesp) = szalayesp(nlayesp) *180.d0/3.141592 ! [deg] + end do + + ! 5th zone: Lower part of ray, decreasing from z to zmin + ! ( layers with SZA > 90 deg ) + do j=grid_R8(zz,diz,nlayer)-1, jmin, -1 + nlayesp = nlayesp + 1 + ilayesp(nlayesp) = j + esp(nlayesp) = + $ sqrt( (radio+diz(j+1))**2 - rkmmini**2 ) + $ - sqrt( (radio+diz(j))**2 - rkmmini**2 ) ! [km] + esp(nlayesp) = esp(nlayesp) * 1.d5 ! [cm] + rkmj = radio+diz(j) + szalayesp(nlayesp) = 3.141592 - asin( rkmmini/rkmj ) ! [rad] + szalayesp(nlayesp) = szalayesp(nlayesp)*180.d0/3.141592 ! [deg] + end do + + end if + + end if + + + end subroutine espesor_optico_A + + + +c********************************************************************** +c*********************************************************************** + + function grid_R8 (z, zgrid, nz) + +c Returns the index where z is located within vector zgrid +c The vector zgrid must be monotonously increasing, otherwise program stops. +c If z is outside zgrid limits, or zgrid dimension is nz<2, the program stops. +c +c FGG Aug-2004 Correct z.lt.zgrid(i) to .le. +c MALV Jul-2003 +c*********************************************************************** + + implicit none + +c Arguments + integer nz + real*8 z + real*8 zgrid(nz) + integer grid_R8 + +c Local + integer i, nz1, nznew + +c*** CODE START + + if ( z .lt. zgrid(1) ) then + write (*,*) ' GRID/ z outside bounds of zgrid ' + write (*,*) ' z,zgrid(1),zgrid(nz) =', z,zgrid(1),zgrid(nz) + z = zgrid(1) + write(*,*) 'WARNING: error in grid_r8 (jthermcalc.F)' + write(*,*) 'Please check values of z and zgrid above' + endif + if (z .gt. zgrid(nz) ) then + write (*,*) ' GRID/ z outside bounds of zgrid ' + write (*,*) ' z,zgrid(1),zgrid(nz) =', z,zgrid(1),zgrid(nz) + z = zgrid(nz) + write(*,*) 'WARNING: error in grid_r8 (jthermcalc.F)' + write(*,*) 'Please check values of z and zgrid above' + endif + if ( nz .lt. 2 ) then + write (*,*) ' GRID/ zgrid needs 2 points at least ! ' + stop ' Serious error in GRID.F ' + endif + if ( zgrid(1) .ge. zgrid(nz) ) then + write (*,*) ' GRID/ zgrid must increase with index' + stop ' Serious error in GRID.F ' + endif + + nz1 = 1 + nznew = nz/2 + if ( z .gt. zgrid(nznew) ) then + nz1 = nznew + nznew = nz + endif + do i=nz1+1,nznew + if ( z. eq. zgrid(i) ) then + grid_R8=i + return + elseif ( z .le. zgrid(i) ) then + grid_R8 = i-1 + return + endif + enddo + grid_R8 = nz + + + + end function grid_R8 + + + +!c*************************************************** +!c*************************************************** + + subroutine flujo(date) + + +!c fgg nov 2002 first version +!c*************************************************** + + use comsaison_h, only: dist_sol + use param_v4_h, only: ninter, + . fluxtop, ct1, ct2, p1, p2 + implicit none + + +! common variables and constants + include "callkeys.h" + + +! Arguments + + real date + + +! Local variable and constants + + integer i + integer inter + real nada + +!c************************************************* + + if(date.lt.1985.) date=1985. + if(date.gt.2001.) date=2001. + + do i=1,ninter + fluxtop(i)=1. + !Variation of solar flux with 11 years solar cycle + !For more details, see Gonzalez-Galindo et al. 2005 + !To be improved in next versions + if(i.le.24.and.solvarmod.eq.0) then + fluxtop(i)=(((ct1(i)+p1(i)*date)/2.) + $ *sin(2.*3.1416/11.*(date-1985.-3.1416)) + $ +(ct2(i)+p2(i)*date)+1.)*fluxtop(i) + end if + fluxtop(i)=fluxtop(i)*(1.52/dist_sol)**2 + end do + + end subroutine flujo + + end module jthermcalc_util Index: trunk/LMDZ.MARS/libf/aeronomars/paramfoto_compact.F =================================================================== --- trunk/LMDZ.MARS/libf/aeronomars/paramfoto_compact.F (revision 3462) +++ trunk/LMDZ.MARS/libf/aeronomars/paramfoto_compact.F (revision 3464) @@ -1,2 +1,8 @@ + MODULE paramfoto_compact_mod + + IMPLICIT NONE + + CONTAINS + c********************************************************************** @@ -639,8 +645,5 @@ cccccc End altitude loop - return - - - end + end subroutine paramfoto_compact @@ -693,8 +696,6 @@ if ( c_output.lt.1.d-30) c_output=1.d-30 - - return c END - end + end subroutine implicito @@ -772,5 +773,5 @@ return - end + end function ionsec_nplus @@ -840,5 +841,5 @@ return - end + end function ionsec_n2plus @@ -915,5 +916,5 @@ return - end + end function ionsec_oplus @@ -993,5 +994,5 @@ return - end + end function ionsec_coplus @@ -1069,5 +1070,5 @@ return - end + end function ionsec_co2plus @@ -1139,5 +1140,5 @@ return - end + end function ionsec_o2plus @@ -1346,8 +1347,5 @@ - return - - - end + end subroutine phdisrate @@ -1989,6 +1987,5 @@ endif !Of if(chemthermod.eq.3) - return - end + end subroutine getch @@ -2050,21 +2047,21 @@ - external ionsec_nplus - real*8 ionsec_nplus - - external ionsec_n2plus - real*8 ionsec_n2plus - - external ionsec_oplus - real*8 ionsec_oplus +! external ionsec_nplus +! real*8 ionsec_nplus + +! external ionsec_n2plus +! real*8 ionsec_n2plus + +! external ionsec_oplus +! real*8 ionsec_oplus - external ionsec_coplus - real*8 ionsec_coplus - - external ionsec_co2plus - real*8 ionsec_co2plus - - external ionsec_o2plus - real*8 ionsec_o2plus +! external ionsec_coplus +! real*8 ionsec_coplus + +! external ionsec_co2plus +! real*8 ionsec_co2plus + +! external ionsec_o2plus +! real*8 ionsec_o2plus @@ -2688,7 +2685,6 @@ endif !Of chemthermod.eq.3 - return c END - end + end subroutine lifetimes @@ -2868,7 +2864,6 @@ deltat = tminaux - return c END - end + end subroutine timemarching @@ -2925,28 +2920,6 @@ integer j - external ionsec_nplus - real*8 ionsec_nplus - - external ionsec_n2plus - real*8 ionsec_n2plus - - external ionsec_oplus - real*8 ionsec_oplus - - external ionsec_coplus - real*8 ionsec_coplus - - external ionsec_co2plus - real*8 ionsec_co2plus - - external ionsec_o2plus - real*8 ionsec_o2plus - - logical firstcall - save firstcall - + logical,save :: firstcall=.true. !$OMP THREADPRIVATE(firstcall) - - data firstcall /.true./ ccccccccccccccc CODE STARTS @@ -4100,7 +4073,6 @@ endif !Of chemthermod.eq.3 - return c END - end + end subroutine prodsandlosses @@ -4243,32 +4215,4 @@ parameter (cociopt=1) -c external functions -c - external cociente - real*8 cociente - - external ionsec_nplus - real*8 ionsec_nplus - - external ionsec_n2plus - real*8 ionsec_n2plus - - external ionsec_oplus - real*8 ionsec_oplus - - external ionsec_coplus - real*8 ionsec_coplus - - external ionsec_co2plus - real*8 ionsec_co2plus - - external ionsec_o2plus - real*8 ionsec_o2plus - - external avg - real*8 avg - - external dif - real*8 dif real*8 log1 @@ -5083,6 +5027,5 @@ - return - end + end subroutine EF_oscilacion !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! @@ -5093,5 +5036,5 @@ avg = (log1+log2+log3)*0.333 return - end + end function avg !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! function dif(log1,log2,log3,avg) @@ -5104,5 +5047,5 @@ & abs(log3-avg) ) * 0.333 return - end + end function dif !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! @@ -5177,3 +5120,5 @@ return c END - end + end function cociente + + END MODULE paramfoto_compact_mod Index: trunk/LMDZ.MARS/libf/aeronomars/photochemistry.F90 =================================================================== --- trunk/LMDZ.MARS/libf/aeronomars/photochemistry.F90 (revision 3462) +++ trunk/LMDZ.MARS/libf/aeronomars/photochemistry.F90 (revision 3464) @@ -1,2 +1,8 @@ +module photochemistry_mod + +implicit none + +contains + !**************************************************************** ! @@ -23,4 +29,6 @@ use param_v4_h, only: jion +use jthermcalc_e107_mod, only: jthermcalc_e107 +use paramfoto_compact_mod, only: phdisrate implicit none @@ -4339,2 +4347,4 @@ end subroutine photochemistry + +end module photochemistry_mod Index: trunk/LMDZ.MARS/libf/aeronomars/thermosphere.F =================================================================== --- trunk/LMDZ.MARS/libf/aeronomars/thermosphere.F (revision 3462) +++ trunk/LMDZ.MARS/libf/aeronomars/thermosphere.F (revision 3464) @@ -18,4 +18,5 @@ use moldiff_red_mod, only: moldiff_red ! old molecular diffusion scheme use moldiff_MPF_mod, only: moldiff_MPF ! new molecular diffusion scheme + use euvheat_mod, only: euvheat use conc_mod, only: rnew, cpnew USE comcstfi_h, only: r, cpp