| 1 | MODULE newsedim_mod |
|---|
| 2 | |
|---|
| 3 | IMPLICIT NONE |
|---|
| 4 | |
|---|
| 5 | CONTAINS |
|---|
| 6 | |
|---|
| 7 | SUBROUTINE newsedim(ngrid,nlay,naersize,nrhosize,ptimestep, |
|---|
| 8 | & pplev,masse,epaisseur,pt,rd,rho,pqi,wq,beta) |
|---|
| 9 | USE comcstfi_h, ONLY: r,g |
|---|
| 10 | IMPLICIT NONE |
|---|
| 11 | |
|---|
| 12 | c======================================================================= |
|---|
| 13 | c |
|---|
| 14 | c Compute sedimentation of 1 tracer |
|---|
| 15 | c of radius rd (m) and density rho (kg.m-3) |
|---|
| 16 | c |
|---|
| 17 | c======================================================================= |
|---|
| 18 | |
|---|
| 19 | c----------------------------------------------------------------------- |
|---|
| 20 | c declarations: |
|---|
| 21 | c ------------- |
|---|
| 22 | |
|---|
| 23 | c |
|---|
| 24 | c arguments: |
|---|
| 25 | c ---------- |
|---|
| 26 | |
|---|
| 27 | INTEGER,INTENT(IN) :: ngrid,nlay,naersize,nrhosize |
|---|
| 28 | REAL,INTENT(IN) :: ptimestep ! pas de temps physique (s) |
|---|
| 29 | REAL,INTENT(IN) :: pplev(ngrid,nlay+1) ! pression aux inter-couches (Pa) |
|---|
| 30 | REAL,INTENT(IN) :: pt(ngrid,nlay) ! temperature au centre des couches (K) |
|---|
| 31 | real,intent(in) :: masse (ngrid,nlay) ! masse d'une couche (kg) |
|---|
| 32 | real,intent(in) :: epaisseur (ngrid,nlay) ! epaisseur d'une couche (m) |
|---|
| 33 | real,intent(in) :: rd(naersize) ! particle radius (m) |
|---|
| 34 | real,intent(in) :: rho(nrhosize) ! particle density (kg.m-3) |
|---|
| 35 | |
|---|
| 36 | |
|---|
| 37 | c Traceurs : |
|---|
| 38 | real,intent(inout) :: pqi(ngrid,nlay) ! traceur (e.g. ?/kg) |
|---|
| 39 | real,intent(out) :: wq(ngrid,nlay+1) ! flux de traceur durant timestep (?/m-2) |
|---|
| 40 | real,intent(in) :: beta ! correction for the shape of the particles |
|---|
| 41 | ! (see Murphy et al. JGR 1990 vol.95) |
|---|
| 42 | ! beta=1 for spheres |
|---|
| 43 | ! beta=0.85 for irregular particles |
|---|
| 44 | ! beta=0.5 for disk shaped particles |
|---|
| 45 | |
|---|
| 46 | c local: |
|---|
| 47 | c ------ |
|---|
| 48 | |
|---|
| 49 | INTEGER l,ig, k, i |
|---|
| 50 | REAL rfall,rhofall |
|---|
| 51 | |
|---|
| 52 | LOGICAL,SAVE :: firstcall=.true. |
|---|
| 53 | |
|---|
| 54 | c Traceurs : |
|---|
| 55 | c ~~~~~~~~ |
|---|
| 56 | real traversee (ngrid,nlay) |
|---|
| 57 | real vstokes(ngrid,nlay) |
|---|
| 58 | real w(ngrid,nlay) |
|---|
| 59 | real ptop, dztop, Ep, Stra |
|---|
| 60 | |
|---|
| 61 | |
|---|
| 62 | c Physical constant |
|---|
| 63 | c ~~~~~~~~~~~~~~~~~ |
|---|
| 64 | c Gas molecular viscosity (N.s.m-2) |
|---|
| 65 | real,parameter :: visc=1.e-5 ! CO2 |
|---|
| 66 | c Effective gas molecular radius (m) |
|---|
| 67 | real,parameter :: molrad=2.2e-10 ! CO2 |
|---|
| 68 | |
|---|
| 69 | c local and saved variable |
|---|
| 70 | real,save :: a,b |
|---|
| 71 | |
|---|
| 72 | |
|---|
| 73 | c ** un petit test de coherence |
|---|
| 74 | c -------------------------- |
|---|
| 75 | |
|---|
| 76 | ! AS: OK firstcall absolute |
|---|
| 77 | IF (firstcall) THEN |
|---|
| 78 | |
|---|
| 79 | firstcall=.false. |
|---|
| 80 | |
|---|
| 81 | |
|---|
| 82 | c Preliminary calculations for sedimenation velocity : |
|---|
| 83 | c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
|---|
| 84 | |
|---|
| 85 | c - Constant to compute stokes speed simple formulae |
|---|
| 86 | c (Vstokes = b * rho* r**2 avec b= (2/9) * rho * g / visc |
|---|
| 87 | b = 2./9. * g / visc |
|---|
| 88 | ENDIF ! of IF(firstcall) |
|---|
| 89 | |
|---|
| 90 | c - Constant to compute gas mean free path |
|---|
| 91 | c l= (T/P)*a, avec a = ( 0.707*8.31/(4*pi*molrad**2 * avogadro)) |
|---|
| 92 | a = 0.707*8.31/(4*3.1416* molrad*molrad * 6.023e23) |
|---|
| 93 | |
|---|
| 94 | c - Correction to account for non-spherical shape (Murphy et al. 1990) |
|---|
| 95 | a = a * beta |
|---|
| 96 | |
|---|
| 97 | |
|---|
| 98 | |
|---|
| 99 | c----------------------------------------------------------------------- |
|---|
| 100 | c 1. initialisation |
|---|
| 101 | c ----------------- |
|---|
| 102 | |
|---|
| 103 | c Sedimentation velocity (m/s) |
|---|
| 104 | c ~~~~~~~~~~~~~~~~~~~~~~ |
|---|
| 105 | c (stokes law corrected for low pressure by the Cunningham |
|---|
| 106 | c slip-flow correction according to Rossow (Icarus 36, 1-50, 1978) |
|---|
| 107 | |
|---|
| 108 | do l=1,nlay |
|---|
| 109 | do ig=1, ngrid |
|---|
| 110 | c radius |
|---|
| 111 | if (naersize.eq.1) then |
|---|
| 112 | rfall=rd(1) |
|---|
| 113 | else |
|---|
| 114 | i=ngrid*(l-1)+ig |
|---|
| 115 | rfall=rd(i) |
|---|
| 116 | endif |
|---|
| 117 | c density |
|---|
| 118 | if (nrhosize.eq.1) then |
|---|
| 119 | rhofall=rho(1) |
|---|
| 120 | else |
|---|
| 121 | i=ngrid*(l-1)+ig |
|---|
| 122 | rhofall=rho(i) |
|---|
| 123 | endif |
|---|
| 124 | c vstokes |
|---|
| 125 | vstokes(ig,l) = b * rhofall * rfall*rfall * |
|---|
| 126 | & (1 + 1.333* ( a*pt(ig,l)/pplev(ig,l) )/rfall) |
|---|
| 127 | |
|---|
| 128 | c Layer crossing time (s) : |
|---|
| 129 | traversee(ig,l)= epaisseur(ig,l)/vstokes(ig,l) |
|---|
| 130 | end do |
|---|
| 131 | end do |
|---|
| 132 | |
|---|
| 133 | |
|---|
| 134 | c Calcul de la masse d'atmosphere correspondant a q transferee |
|---|
| 135 | c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
|---|
| 136 | c (e.g. on recherche le niveau en dessous de laquelle le traceur |
|---|
| 137 | c va traverser le niveau intercouche l : "dztop" est sa hauteur |
|---|
| 138 | c au dessus de l (m), "ptop" est sa pression (Pa)) |
|---|
| 139 | |
|---|
| 140 | do l=1,nlay |
|---|
| 141 | do ig=1, ngrid |
|---|
| 142 | |
|---|
| 143 | dztop = vstokes(ig,l)* ptimestep |
|---|
| 144 | Ep=0 |
|---|
| 145 | k=0 |
|---|
| 146 | |
|---|
| 147 | w(ig,l) = 0. !! JF+AS ajout initialisation |
|---|
| 148 | c ************************************************************** |
|---|
| 149 | c Simple Method |
|---|
| 150 | |
|---|
| 151 | cc w(ig,l) = |
|---|
| 152 | cc & (1.- exp(-dztop*g/(r*pt(ig,l))))*pplev(ig,l) / g |
|---|
| 153 | cccc write(*,*) 'OK simple method l,w =', l, w(ig,l) |
|---|
| 154 | cccc write(*,*) 'OK simple method dztop =', dztop |
|---|
| 155 | |
|---|
| 156 | w(ig,l) = 1. - exp(-dztop*g/(r*pt(ig,l))) |
|---|
| 157 | !!! Diagnostic: JF. Fix: AS. Date: 05/11 |
|---|
| 158 | !!! Probleme arrondi avec la quantite ci-dessus |
|---|
| 159 | !!! ---> vaut 0 pour -dztop*g/(r*pt(ig,l)) trop petit |
|---|
| 160 | !!! ---> dans ce cas on utilise le developpement limite ! |
|---|
| 161 | !!! ---> exp(-x) = 1 - x lorsque x --> 0 avec une erreur de x^2 / 2 |
|---|
| 162 | IF ( w(ig,l) .eq. 0. ) THEN |
|---|
| 163 | w(ig,l) = ( dztop*g/(r*pt(ig,l)) ) * pplev(ig,l) / g |
|---|
| 164 | ELSE |
|---|
| 165 | w(ig,l) = w(ig,l) * pplev(ig,l) / g |
|---|
| 166 | ENDIF |
|---|
| 167 | |
|---|
| 168 | |
|---|
| 169 | c ************************************************************** |
|---|
| 170 | cccc Complex method : |
|---|
| 171 | if (dztop.gt.epaisseur(ig,l)) then !!!if on traverse plus d'une couche |
|---|
| 172 | cccc Cas ou on "epuise" la couche l : On calcule le flux |
|---|
| 173 | cccc Venant de dessus en tenant compte de la variation de Vstokes |
|---|
| 174 | c ************************************************************** |
|---|
| 175 | Ep= epaisseur(ig,l) |
|---|
| 176 | Stra= traversee(ig,l) |
|---|
| 177 | do while(dztop.gt.Ep.and.l+k+1.le.nlay) |
|---|
| 178 | k=k+1 |
|---|
| 179 | dztop= Ep + vstokes(ig,l+k)*(ptimestep -Stra) |
|---|
| 180 | Ep = Ep + epaisseur(ig,l+k) |
|---|
| 181 | Stra = Stra + traversee(ig,l+k) |
|---|
| 182 | enddo |
|---|
| 183 | Ep = Ep - epaisseur(ig,l+k) |
|---|
| 184 | !ptop=pplev(ig,l+k)*exp(-(dztop-Ep)*g/(r*pt(ig,l+k))) |
|---|
| 185 | |
|---|
| 186 | !!! JF+AS 05/11 Probleme arrondi potentiel, meme solution que ci-dessus |
|---|
| 187 | ptop=exp(-(dztop-Ep)*g/(r*pt(ig,l+k))) |
|---|
| 188 | IF ( ptop .eq. 1. ) THEN |
|---|
| 189 | !PRINT*, 'newsedim: exposant trop petit ', ig, l |
|---|
| 190 | ptop=pplev(ig,l+k) * ( 1. - (dztop-Ep)*g/(r*pt(ig,l+k))) |
|---|
| 191 | ELSE |
|---|
| 192 | ptop=pplev(ig,l+k) * ptop |
|---|
| 193 | ENDIF |
|---|
| 194 | w(ig,l) = (pplev(ig,l) - Ptop)/g |
|---|
| 195 | |
|---|
| 196 | endif !!!!!if complex method |
|---|
| 197 | |
|---|
| 198 | |
|---|
| 199 | cc write(*,*) 'OK new method l,w =', l, w(ig,l) |
|---|
| 200 | cc write(*,*) 'OK new method dztop =', dztop |
|---|
| 201 | cc if(l.eq.7)write(*,*)'l=7,k,pplev,Ptop',pplev(ig,l),Ptop |
|---|
| 202 | cc if(l.eq.7)write(*,*)'l=7,dztop,Ep',dztop,Ep |
|---|
| 203 | cc if(l.eq.6)write(*,*)'l=6,k, w',k, w(1,l) |
|---|
| 204 | cc if(l.eq.7)write(*,*)'l=7,k, w',k, w(1,l) |
|---|
| 205 | cc if(l.eq.8)write(*,*)'l=8,k, w',k, w(1,l) |
|---|
| 206 | c ************************************************************** |
|---|
| 207 | |
|---|
| 208 | |
|---|
| 209 | end do |
|---|
| 210 | end do |
|---|
| 211 | call vlz_fi(ngrid,nlay,pqi,2.,masse,w,wq) |
|---|
| 212 | c write(*,*) ' newsed: wq(6), wq(7), q(6)', |
|---|
| 213 | c & wq(1,6),wq(1,7),pqi(1,6) |
|---|
| 214 | |
|---|
| 215 | |
|---|
| 216 | END SUBROUTINE newsedim |
|---|
| 217 | |
|---|
| 218 | END MODULE newsedim_mod |
|---|
| 219 | |
|---|
