| 1 | SUBROUTINE newsedim(ngrid,nlay,naersize,ptimestep, |
|---|
| 2 | & pplev,masse,epaisseur,pt,rd,rho,pqi,wq,iq,pphi) |
|---|
| 3 | |
|---|
| 4 | use ioipsl_getin_p_mod, only: getin_p |
|---|
| 5 | use comcstfi_mod, only: r, g, rad |
|---|
| 6 | use gases_h |
|---|
| 7 | ! use tracer_h, only : igcm_h2o_ice |
|---|
| 8 | ! use watercommon_h, only: T_h2O_ice_liq,T_h2O_ice_clouds |
|---|
| 9 | ! use radii_mod, only: h2o_cloudrad |
|---|
| 10 | |
|---|
| 11 | IMPLICIT NONE |
|---|
| 12 | |
|---|
| 13 | !================================================================== |
|---|
| 14 | ! |
|---|
| 15 | ! Purpose |
|---|
| 16 | ! ------- |
|---|
| 17 | ! Calculates sedimentation of 1 tracer |
|---|
| 18 | ! of radius rd (m) and density rho (kg.m-3) |
|---|
| 19 | ! |
|---|
| 20 | !================================================================== |
|---|
| 21 | |
|---|
| 22 | !----------------------------------------------------------------------- |
|---|
| 23 | ! declarations |
|---|
| 24 | ! ------------ |
|---|
| 25 | |
|---|
| 26 | ! arguments |
|---|
| 27 | ! --------- |
|---|
| 28 | |
|---|
| 29 | integer,intent(in) :: ngrid ! number of atmospheric columns |
|---|
| 30 | integer,intent(in) :: nlay ! number of atmospheric layers |
|---|
| 31 | integer,intent(in) :: naersize ! number of particle sizes (1 or number |
|---|
| 32 | ! of grid boxes) |
|---|
| 33 | real,intent(in) :: ptimestep ! physics time step (s) |
|---|
| 34 | real,intent(in) :: pplev(ngrid,nlay+1) ! inter-layer pressures (Pa) |
|---|
| 35 | real,intent(in) :: pt(ngrid,nlay) ! mid-layer temperatures (K) |
|---|
| 36 | real,intent(in) :: masse (ngrid,nlay) ! atmospheric mass (kg) |
|---|
| 37 | real,intent(in) :: epaisseur (ngrid,nlay) ! thickness of atm. layers (m) |
|---|
| 38 | real,intent(in) :: rd(naersize) ! particle radius (m) |
|---|
| 39 | real,intent(in) :: rho ! particle density (kg.m-3) |
|---|
| 40 | real,intent(inout) :: pqi(ngrid,nlay) ! tracer (e.g. ?/kg) |
|---|
| 41 | real,intent(out) :: wq(ngrid,nlay+1) ! flux of tracer during timestep (?/m-2) |
|---|
| 42 | integer,intent(in) :: iq ! tracer index |
|---|
| 43 | real,intent(in):: pphi(ngrid,nlay) ! geopotential |
|---|
| 44 | |
|---|
| 45 | c local: |
|---|
| 46 | c ------ |
|---|
| 47 | |
|---|
| 48 | INTEGER l,ig, k, i, igas |
|---|
| 49 | REAL rfall, rsurf, Reynolds, Cd, zfrac |
|---|
| 50 | REAL reffh2oliq(ngrid,nlay), reffh2oice(ngrid,nlay) |
|---|
| 51 | |
|---|
| 52 | LOGICAL,SAVE :: firstcall=.true. |
|---|
| 53 | !$OMP THREADPRIVATE(firstcall) |
|---|
| 54 | LOGICAL,SAVE :: crystal_shape |
|---|
| 55 | !$OMP THREADPRIVATE(crystal_shape) |
|---|
| 56 | |
|---|
| 57 | c Traceurs : |
|---|
| 58 | c ~~~~~~~~ |
|---|
| 59 | real traversee (ngrid,nlay) |
|---|
| 60 | real vstokes(ngrid,nlay) |
|---|
| 61 | real w(ngrid,nlay) |
|---|
| 62 | real ptop, dztop, Ep, Stra |
|---|
| 63 | |
|---|
| 64 | |
|---|
| 65 | c Physical constant |
|---|
| 66 | c ~~~~~~~~~~~~~~~~~ |
|---|
| 67 | c Gas molecular viscosity (N.s.m-2) |
|---|
| 68 | real, allocatable, save :: visc(:,:) |
|---|
| 69 | !$OMP THREADPRIVATE(visc) |
|---|
| 70 | c Effective gas molecular radius (m) |
|---|
| 71 | real,save :: molrad |
|---|
| 72 | |
|---|
| 73 | c local and saved variable |
|---|
| 74 | real,save :: a,b |
|---|
| 75 | !$OMP THREADPRIVATE(a,b) |
|---|
| 76 | |
|---|
| 77 | c ** un petit test de coherence |
|---|
| 78 | c -------------------------- |
|---|
| 79 | |
|---|
| 80 | !print*,'temporary fixed particle rad in newsedim!!' |
|---|
| 81 | |
|---|
| 82 | IF (firstcall) THEN |
|---|
| 83 | firstcall=.false. |
|---|
| 84 | |
|---|
| 85 | c Determination of the viscosity a(N.s.m-2) and the mean molecular radius (m) |
|---|
| 86 | write(*,*) "Calculation of the viscosity and the mean molecular" |
|---|
| 87 | & ," radius from gases.def" |
|---|
| 88 | allocate(visc(ngrid,nlay)) |
|---|
| 89 | visc(:,:)=0.0 |
|---|
| 90 | molrad=0. |
|---|
| 91 | do igas=1, ngasmx |
|---|
| 92 | if(gfrac(igas).ge.0.0) then |
|---|
| 93 | if(igas.eq.igas_CO2) then |
|---|
| 94 | molrad = molrad + gfrac(igas)*2.2e-10 ! CO2 |
|---|
| 95 | visc(:,:) = visc(:,:) + gfrac(igas)*1.0e-5 ! CO2 |
|---|
| 96 | elseif(igas.eq.igas_N2) then |
|---|
| 97 | molrad = molrad + gfrac(igas)*1.8e-10 ! N2 (Kunze et al. 2022) |
|---|
| 98 | visc(:,:) = visc(:,:) + gfrac(igas)*1.0e-5 ! N2 |
|---|
| 99 | elseif(igas.eq.igas_H2) then |
|---|
| 100 | molrad = molrad + gfrac(igas)*1.41e-10 ! H2 (Ackerman & Marley 2001) |
|---|
| 101 | visc(:,:) = visc(:,:) + gfrac(igas)*2.0d-07*pt(:,:)**0.66 ! H2 (from Rosner 2000) |
|---|
| 102 | elseif(igas.eq.igas_H2O) then |
|---|
| 103 | molrad = molrad + gfrac(igas)*2.3e-10 ! H2O (Crifo 1989 at 300K) |
|---|
| 104 | visc(:,:) = visc(:,:) + gfrac(igas)*8e-6 ! H2O (Sengers & Kamgar-Parsi 1984) |
|---|
| 105 | elseif(igas.eq.igas_He) then |
|---|
| 106 | molrad = molrad + gfrac(igas)*1.1e-10 ! He (Kunze et al. 2022) |
|---|
| 107 | visc(:,:) = visc(:,:) + ! He |
|---|
| 108 | & gfrac(igas)*1.9e-5*(pt(:,:)/273.15)**0.7 ! He (Petersen 1970) |
|---|
| 109 | elseif(igas.eq.igas_CH4) then |
|---|
| 110 | molrad = molrad + gfrac(igas)*1.9e-10 ! CH4 (Ismail et al. 2015) |
|---|
| 111 | visc(:,:) = visc(:,:) + gfrac(igas)*1.0e-5 ! CH4 |
|---|
| 112 | else |
|---|
| 113 | molrad = molrad + gfrac(igas)*1.93e-10 ! N2 by default |
|---|
| 114 | visc(:,:) = visc(:,:) + 6.67e-6 ! N2 by default |
|---|
| 115 | write(*,*) trim(gnom(igas))," is not included in" |
|---|
| 116 | & ," newsedim, N2 is used by default" |
|---|
| 117 | endif |
|---|
| 118 | endif |
|---|
| 119 | enddo |
|---|
| 120 | write(*,*) "visc(1,1)=",visc(1,1),"N.s.m-2; ", |
|---|
| 121 | & "molrad=",molrad,"m" |
|---|
| 122 | |
|---|
| 123 | c Correction for non-spherical water ice particles |
|---|
| 124 | write(*,*) "Use non-spherical water ice particles", |
|---|
| 125 | & " for the sedimentation ?" |
|---|
| 126 | crystal_shape=.false. !default |
|---|
| 127 | call getin_p("crystal_shape",crystal_shape) |
|---|
| 128 | write(*,*) " crystal_shape = ",crystal_shape |
|---|
| 129 | |
|---|
| 130 | |
|---|
| 131 | !======================================================================= |
|---|
| 132 | ! Preliminary calculations for sedimenation velocity |
|---|
| 133 | |
|---|
| 134 | ! - Constant to compute stokes speed simple formulae |
|---|
| 135 | ! (Vstokes = b / visc * rho* r**2 avec b= (2/9) * rho * g |
|---|
| 136 | b = 2./9. * g |
|---|
| 137 | |
|---|
| 138 | ! - Constant to compute gas mean free path |
|---|
| 139 | ! l= (T/P)*a, avec a = ( 0.707*8.31/(4*pi*molrad**2 * avogadro)) |
|---|
| 140 | a = 0.707*8.31/(4*3.1416* molrad**2 * 6.023e23) |
|---|
| 141 | |
|---|
| 142 | c - Correction to account for non-spherical shape (Murphy et al. 1990) |
|---|
| 143 | c (correction = 0.85 for irregular particles, 0.5 for disk shaped particles) |
|---|
| 144 | c a = a * 0.85 |
|---|
| 145 | |
|---|
| 146 | |
|---|
| 147 | ENDIF |
|---|
| 148 | |
|---|
| 149 | c----------------------------------------------------------------------- |
|---|
| 150 | c 1. initialisation |
|---|
| 151 | c ----------------- |
|---|
| 152 | |
|---|
| 153 | c Sedimentation velocity (m/s) |
|---|
| 154 | c ~~~~~~~~~~~~~~~~~~~~~~ |
|---|
| 155 | c (stokes law corrected for low pressure by the Cunningham |
|---|
| 156 | c slip-flow correction according to Rossow (Icarus 36, 1-50, 1978) |
|---|
| 157 | |
|---|
| 158 | c Compute liquid and ice particle radii |
|---|
| 159 | ! if((iq.eq.igcm_h2o_ice).and.crystal_shape) then |
|---|
| 160 | ! call h2o_cloudrad(ngrid,nlay,pqi,reffh2oliq,reffh2oice) |
|---|
| 161 | ! endif |
|---|
| 162 | |
|---|
| 163 | |
|---|
| 164 | do l=1,nlay |
|---|
| 165 | do ig=1, ngrid |
|---|
| 166 | if (naersize.eq.1) then |
|---|
| 167 | rfall=rd(1) |
|---|
| 168 | else |
|---|
| 169 | i=ngrid*(l-1)+ig |
|---|
| 170 | rfall=rd(i) ! how can this be correct!!? |
|---|
| 171 | endif |
|---|
| 172 | |
|---|
| 173 | c Correction for non-spherical water ice particles |
|---|
| 174 | ! if((iq.eq.igcm_h2o_ice).and.crystal_shape) then |
|---|
| 175 | ! zfrac= (pt(ig,l)-T_h2O_ice_clouds) / |
|---|
| 176 | ! & (T_h2O_ice_liq-T_h2O_ice_clouds) |
|---|
| 177 | ! zfrac= MAX(zfrac, 0.0) |
|---|
| 178 | ! zfrac= MIN(zfrac, 1.0) |
|---|
| 179 | ! rsurf=max(reffh2oice(ig,l),45.6*reffh2oice(ig,l)**1.3) ! surface radius (formula for rimed dendrites from Hemsfield 1977, transition at around 30 microns) |
|---|
| 180 | ! rsurf=1/(zfrac/reffh2oice(ig,l)+(1-zfrac)/rsurf) ! radius giving the mean velocity between liquid and ice particles |
|---|
| 181 | ! else |
|---|
| 182 | rsurf=rfall |
|---|
| 183 | ! endif |
|---|
| 184 | |
|---|
| 185 | ! b = 2./9. * g |
|---|
| 186 | b = 2./9. * ((g*rad-pphi(ig,l))**2/(g*(rad**2))) ! AF24: from Pluto.old |
|---|
| 187 | |
|---|
| 188 | vstokes(ig,l) = b / visc(ig,l) * rho * rfall**3 / rsurf * |
|---|
| 189 | & (1 + 1.333* ( a*pt(ig,l)/pplev(ig,l) )/rsurf) |
|---|
| 190 | |
|---|
| 191 | c Correction for high Reynolds number |
|---|
| 192 | Reynolds=2. * pplev(ig,l) / r / pt(ig,l) * |
|---|
| 193 | & rsurf * vstokes(ig,l) / visc(ig,l) |
|---|
| 194 | if(Reynolds.ge.1.0) then |
|---|
| 195 | do i=1,5 |
|---|
| 196 | Cd=24. / Reynolds * (1. + 0.15 * Reynolds**0.687) + |
|---|
| 197 | & 0.42 / (1. + 42500 / Reynolds**1.16) ! (Formula from Bagheri 2018) |
|---|
| 198 | vstokes(ig,l) =(8./3.*pplev(ig,l)/r/pt(ig,l)*g*rfall**3 / |
|---|
| 199 | & rsurf**2/rho/Cd * |
|---|
| 200 | & (1.+1.333*(a*pt(ig,l)/pplev(ig,l))/rsurf))**0.5 |
|---|
| 201 | Reynolds=2. * pplev(ig,l) / r / pt(ig,l) * |
|---|
| 202 | & rsurf * vstokes(ig,l) / visc(ig,l) |
|---|
| 203 | enddo |
|---|
| 204 | endif |
|---|
| 205 | |
|---|
| 206 | c Layer crossing time (s) : |
|---|
| 207 | traversee(ig,l)= epaisseur(ig,l)/vstokes(ig,l) |
|---|
| 208 | end do |
|---|
| 209 | end do |
|---|
| 210 | |
|---|
| 211 | |
|---|
| 212 | c Calcul de la masse d'atmosphere correspondant a q transferee |
|---|
| 213 | c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
|---|
| 214 | c (e.g. on recherche le niveau en dessous de laquelle le traceur |
|---|
| 215 | c va traverser le niveau intercouche l : "dztop" est sa hauteur |
|---|
| 216 | c au dessus de l (m), "ptop" est sa pression (Pa)) |
|---|
| 217 | do l=1,nlay |
|---|
| 218 | do ig=1, ngrid |
|---|
| 219 | |
|---|
| 220 | dztop = vstokes(ig,l)* ptimestep |
|---|
| 221 | Ep=0 |
|---|
| 222 | k=0 |
|---|
| 223 | w(ig,l) = 0. !! JF+AS ajout initialisation (LK MARS) |
|---|
| 224 | c ************************************************************** |
|---|
| 225 | c Simple Method |
|---|
| 226 | cc w(ig,l) = |
|---|
| 227 | cc & (1- exp(-dztop*g/(r*pt(ig,l))))*pplev(ig,l) / g |
|---|
| 228 | cc write(*,*) 'OK simple method l,w =', l, w(ig,l) |
|---|
| 229 | cc write(*,*) 'OK simple method dztop =', dztop |
|---|
| 230 | w(ig,l) = 1. - exp(-dztop*g/(r*pt(ig,l))) |
|---|
| 231 | !!! Diagnostic: JF. Fix: AS. Date: 05/11 |
|---|
| 232 | !!! Probleme arrondi avec la quantite ci-dessus |
|---|
| 233 | !!! ---> vaut 0 pour -dztop*g/(r*pt(ig,l)) trop petit |
|---|
| 234 | !!! ---> dans ce cas on utilise le developpement limite ! |
|---|
| 235 | !!! ---> exp(-x) = 1 - x lorsque x --> 0 avec une erreur de x^2 / 2 |
|---|
| 236 | |
|---|
| 237 | IF ( w(ig,l) .eq. 0. ) THEN |
|---|
| 238 | w(ig,l) = ( dztop*g/(r*pt(ig,l)) ) * pplev(ig,l) / g |
|---|
| 239 | ELSE |
|---|
| 240 | w(ig,l) = w(ig,l) * pplev(ig,l) / g |
|---|
| 241 | ENDIF |
|---|
| 242 | ! LK borrowed simple method from Mars model (AS/JF) |
|---|
| 243 | |
|---|
| 244 | !************************************************************** |
|---|
| 245 | cccc Complex method : |
|---|
| 246 | if (dztop.gt.epaisseur(ig,l)) then |
|---|
| 247 | cccc Cas ou on "epuise" la couche l : On calcule le flux |
|---|
| 248 | cccc Venant de dessus en tenant compte de la variation de Vstokes |
|---|
| 249 | |
|---|
| 250 | Ep= epaisseur(ig,l) |
|---|
| 251 | Stra= traversee(ig,l) |
|---|
| 252 | do while(dztop.gt.Ep.and.l+k+1.le.nlay) |
|---|
| 253 | k=k+1 |
|---|
| 254 | dztop= Ep + vstokes(ig,l+k)*(ptimestep -Stra) |
|---|
| 255 | Ep = Ep + epaisseur(ig,l+k) |
|---|
| 256 | Stra = Stra + traversee(ig,l+k) |
|---|
| 257 | enddo |
|---|
| 258 | Ep = Ep - epaisseur(ig,l+k) |
|---|
| 259 | ! ptop=pplev(ig,l+k)*exp(-(dztop-Ep)*g/(r*pt(ig,l+k))) |
|---|
| 260 | ptop=exp(-(dztop-Ep)*g/(r*pt(ig,l+k))) |
|---|
| 261 | IF ( ptop .eq. 1. ) THEN |
|---|
| 262 | !PRINT*, 'newsedim: exposant trop petit ', ig, l |
|---|
| 263 | ptop=pplev(ig,l+k) * ( 1. - (dztop-Ep)*g/(r*pt(ig,l+k))) |
|---|
| 264 | ELSE |
|---|
| 265 | ptop=pplev(ig,l+k) * ptop |
|---|
| 266 | ENDIF |
|---|
| 267 | |
|---|
| 268 | w(ig,l) = (pplev(ig,l) - ptop)/g |
|---|
| 269 | |
|---|
| 270 | endif !!! complex method |
|---|
| 271 | c |
|---|
| 272 | cc write(*,*) 'OK new method l,w =', l, w(ig,l) |
|---|
| 273 | cc write(*,*) 'OK new method dztop =', dztop |
|---|
| 274 | cc if(l.eq.7)write(*,*)'l=7,k,pplev,Ptop',pplev(ig,l),Ptop |
|---|
| 275 | cc if(l.eq.7)write(*,*)'l=7,dztop,Ep',dztop,Ep |
|---|
| 276 | cc if(l.eq.6)write(*,*)'l=6,k, w',k, w(1,l) |
|---|
| 277 | cc if(l.eq.7)write(*,*)'l=7,k, w',k, w(1,l) |
|---|
| 278 | cc if(l.eq.8)write(*,*)'l=8,k, w',k, w(1,l) |
|---|
| 279 | c ************************************************************** |
|---|
| 280 | |
|---|
| 281 | end do |
|---|
| 282 | end do |
|---|
| 283 | |
|---|
| 284 | call vlz_fi(ngrid,nlay,pqi,2.,masse,w,wq) |
|---|
| 285 | c write(*,*) ' newsed: wq(6), wq(7), q(6)', |
|---|
| 286 | c & wq(1,6),wq(1,7),pqi(1,6) |
|---|
| 287 | |
|---|
| 288 | END |
|---|
