| 1 | SUBROUTINE initracer(qsurf,co2ice) |
|---|
| 2 | |
|---|
| 3 | IMPLICIT NONE |
|---|
| 4 | c======================================================================= |
|---|
| 5 | c subject: |
|---|
| 6 | c -------- |
|---|
| 7 | c Initialisation related to tracer |
|---|
| 8 | c (transported dust, water, chemical species, ice...) |
|---|
| 9 | c |
|---|
| 10 | c Name of the tracer |
|---|
| 11 | c |
|---|
| 12 | c Test of dimension : |
|---|
| 13 | c Initialise COMMON tracer in tracer.h |
|---|
| 14 | c |
|---|
| 15 | c If water=T : q(iq=nqmx) is the water mass mixing ratio |
|---|
| 16 | c If water=T and iceparty=T : q(iq=nqmx-1) is the ice mass mixing ratio |
|---|
| 17 | |
|---|
| 18 | c If there is transported dust, it uses iq=1 to iq=dustbin |
|---|
| 19 | c If there is no transported dust : dustbin=0 |
|---|
| 20 | c If doubleq=T : q(iq=1) is the dust mass mixing ratio |
|---|
| 21 | c q(iq=2) is the dust number mixing ratio |
|---|
| 22 | |
|---|
| 23 | c If (photochem.or.thermochem) there is "ncomp" chemical species (ncomp |
|---|
| 24 | c is set in aeronomars/chimiedata.h) using the ncomp iq values starting at |
|---|
| 25 | c iq=nqchem_min = dustbin+1 (nqchem_min is defined in inifis.F) |
|---|
| 26 | c |
|---|
| 27 | c |
|---|
| 28 | c author: F.Forget |
|---|
| 29 | c ------ |
|---|
| 30 | c Modifs: Franck Montmessin, Sebastien Lebonnois (june 2003) |
|---|
| 31 | c |
|---|
| 32 | c======================================================================= |
|---|
| 33 | |
|---|
| 34 | |
|---|
| 35 | #include "dimensions.h" |
|---|
| 36 | #include "dimphys.h" |
|---|
| 37 | #include "comcstfi.h" |
|---|
| 38 | #include "callkeys.h" |
|---|
| 39 | #include "tracer.h" |
|---|
| 40 | |
|---|
| 41 | #include "comgeomfi.h" |
|---|
| 42 | #include "watercap.h" |
|---|
| 43 | #include "aerice.h" |
|---|
| 44 | #include "fisice.h" |
|---|
| 45 | #include "chimiedata.h" |
|---|
| 46 | |
|---|
| 47 | |
|---|
| 48 | real qsurf(ngridmx,nqmx) ! tracer on surface (e.g. kg.m-2) |
|---|
| 49 | real co2ice(ngridmx) ! co2 ice mass on surface (e.g. kg.m-2) |
|---|
| 50 | integer iq,ig |
|---|
| 51 | real r0_lift , reff_lift |
|---|
| 52 | |
|---|
| 53 | c----------------------------------------------------------------------- |
|---|
| 54 | c radius(nqmx) ! aerosol particle radius (m) |
|---|
| 55 | c rho_q(nqmx) ! tracer densities (kg.m-3) |
|---|
| 56 | c qext(nqmx) ! Single Scat. Extinction coeff at 0.67 um |
|---|
| 57 | c alpha_lift(nqmx) ! saltation vertical flux/horiz flux ratio (m-1) |
|---|
| 58 | c alpha_devil(nqmx) ! lifting coeeficient by dust devil |
|---|
| 59 | c rho_dust ! Mars dust density |
|---|
| 60 | c rho_ice ! Water ice density |
|---|
| 61 | c doubleq ! if method with mass (iq=1) and number(iq=2) mixing ratio |
|---|
| 62 | c varian ! Characteristic variance of log-normal distribution |
|---|
| 63 | c----------------------------------------------------------------------- |
|---|
| 64 | |
|---|
| 65 | |
|---|
| 66 | c------------------------------------------------------------ |
|---|
| 67 | c Test Dimensions tracers |
|---|
| 68 | c------------------------------------------------------------ |
|---|
| 69 | |
|---|
| 70 | if(photochem.or.thermochem) then |
|---|
| 71 | if (iceparty) then |
|---|
| 72 | if ((nqchem_min+ncomp+1).ne.nqmx) then |
|---|
| 73 | print*,'********* Dimension problem! ********' |
|---|
| 74 | print*,"nqchem_min+ncomp+1).ne.nqmx" |
|---|
| 75 | print*,"ncomp: ",ncomp |
|---|
| 76 | print*,"nqchem_min: ",nqchem_min |
|---|
| 77 | print*,"nqmx: ",nqmx |
|---|
| 78 | print*,'Change ncomp in chimiedata.h' |
|---|
| 79 | endif |
|---|
| 80 | else |
|---|
| 81 | if ((nqchem_min+ncomp).ne.nqmx) then |
|---|
| 82 | print*,'********* Dimension problem! ********' |
|---|
| 83 | print*,"nqchem_min+ncomp).ne.nqmx" |
|---|
| 84 | print*,"ncomp: ",ncomp |
|---|
| 85 | print*,"nqchem_min: ",nqchem_min |
|---|
| 86 | print*,"nqmx: ",nqmx |
|---|
| 87 | print*,'Change ncomp in chimiedata.h' |
|---|
| 88 | STOP |
|---|
| 89 | endif |
|---|
| 90 | endif |
|---|
| 91 | endif |
|---|
| 92 | |
|---|
| 93 | c------------------------------------------------------------ |
|---|
| 94 | c NAME and molar mass of the tracer |
|---|
| 95 | c------------------------------------------------------------ |
|---|
| 96 | |
|---|
| 97 | c noms and mmol vectors: |
|---|
| 98 | if (water) then |
|---|
| 99 | mmol(nqmx) = 18. |
|---|
| 100 | noms(nqmx) = 'h2o' |
|---|
| 101 | end if |
|---|
| 102 | if (iceparty) then |
|---|
| 103 | noms(nqmx-1) = 'ice' |
|---|
| 104 | mmol(nqmx-1) = 18. |
|---|
| 105 | end if |
|---|
| 106 | if(photochem.or.thermochem) then |
|---|
| 107 | do iq=nqchem_min, nqchem_min+ncomp-1 |
|---|
| 108 | noms(iq) = nomchem(iq-nqchem_min+1) |
|---|
| 109 | mmol(iq) = mmolchem(iq-nqchem_min+1) |
|---|
| 110 | enddo |
|---|
| 111 | end if |
|---|
| 112 | if (dustbin.ge.1) then |
|---|
| 113 | do iq=1,dustbin |
|---|
| 114 | noms(iq) = 'dust' |
|---|
| 115 | mmol(iq) = 100. |
|---|
| 116 | enddo |
|---|
| 117 | if (doubleq) then |
|---|
| 118 | noms(1) = 'dust mass mix. ratio' |
|---|
| 119 | noms(dustbin) = 'dust number mix. ratio' |
|---|
| 120 | end if |
|---|
| 121 | end if |
|---|
| 122 | |
|---|
| 123 | c Simulation of CO2 + neutral gaz |
|---|
| 124 | if ((dustbin.eq.0).and.(.not.water)) then |
|---|
| 125 | noms(1) = 'co2' |
|---|
| 126 | mmol(1) = 44 |
|---|
| 127 | if (nqmx.eq.2)then |
|---|
| 128 | noms(nqmx) = 'Ar_N2' |
|---|
| 129 | mmol(nqmx) = 30 |
|---|
| 130 | end if |
|---|
| 131 | end if |
|---|
| 132 | |
|---|
| 133 | |
|---|
| 134 | c------------------------------------------------------------ |
|---|
| 135 | c Initialisation tracers .... |
|---|
| 136 | c------------------------------------------------------------ |
|---|
| 137 | call zerophys(nqmx,rho_q) |
|---|
| 138 | |
|---|
| 139 | rho_dust=2500. ! Mars dust density (kg.m-3) |
|---|
| 140 | rho_ice=920. ! Water ice density (kg.m-3) |
|---|
| 141 | |
|---|
| 142 | if (doubleq) then |
|---|
| 143 | c "doubleq" technique |
|---|
| 144 | c ------------------- |
|---|
| 145 | c (transport of mass and number mixing ratio) |
|---|
| 146 | c iq=1: Q mass mixing ratio, iq=2: N number mixing ratio |
|---|
| 147 | |
|---|
| 148 | if( (nqmx.lt.2).or.(water.and.(nqmx.lt.3)) ) then |
|---|
| 149 | write(*,*) 'nqmx is too low : nqmx=', nqmx |
|---|
| 150 | write(*,*) 'water= ',water,' doubleq= ',doubleq |
|---|
| 151 | end if |
|---|
| 152 | |
|---|
| 153 | varian=0.637 ! Characteristic variance |
|---|
| 154 | qext(1)=3.04 ! reference extinction at 0.67 um for ref dust |
|---|
| 155 | qext(2)=3.04 ! reference extinction at 0.67 um for ref dust |
|---|
| 156 | rho_q(1)=rho_dust |
|---|
| 157 | rho_q(2)=rho_dust |
|---|
| 158 | |
|---|
| 159 | c Intermediate calcul for computing geometric mean radius r0 |
|---|
| 160 | c as a function of mass and number mixing ratio Q and N |
|---|
| 161 | c (r0 = (r3n_q * Q/ N) |
|---|
| 162 | r3n_q = exp(-4.5*varian**2)*(3./4.)/(pi*rho_dust) |
|---|
| 163 | |
|---|
| 164 | c Intermediate calcul for computing effective radius reff |
|---|
| 165 | c from geometric mean radius r0 |
|---|
| 166 | c (reff = ref_r0 * r0) |
|---|
| 167 | ref_r0 = exp(2.5*varian**2) |
|---|
| 168 | |
|---|
| 169 | c lifted dust : |
|---|
| 170 | c ''''''''''' |
|---|
| 171 | reff_lift = 3.e-6 ! Effective radius of lifted dust (m) |
|---|
| 172 | alpha_devil(1)=9.e-9 ! dust devil lift mass coeff |
|---|
| 173 | alpha_lift(1)=3.0e-15 ! Lifted mass coeff |
|---|
| 174 | |
|---|
| 175 | r0_lift = reff_lift/ref_r0 |
|---|
| 176 | alpha_devil(2)= r3n_q * alpha_devil(1)/r0_lift**3 |
|---|
| 177 | alpha_lift(2)= r3n_q * alpha_lift(1)/r0_lift**3 |
|---|
| 178 | |
|---|
| 179 | c Not used: |
|---|
| 180 | radius(1) = 0. |
|---|
| 181 | radius(2) = 0. |
|---|
| 182 | |
|---|
| 183 | else |
|---|
| 184 | |
|---|
| 185 | if (dustbin.gt.1) then |
|---|
| 186 | print*,'ATTENTION:', |
|---|
| 187 | $ ' properties of dust need input in initracer !!!' |
|---|
| 188 | stop |
|---|
| 189 | |
|---|
| 190 | else if (dustbin.eq.1) then |
|---|
| 191 | |
|---|
| 192 | c This will be used for 1 dust particle size: |
|---|
| 193 | c ------------------------------------------ |
|---|
| 194 | radius(1)=3.e-6 |
|---|
| 195 | Qext(1)=3.04 |
|---|
| 196 | alpha_lift(1)=0.0e-6 |
|---|
| 197 | alpha_devil(1)=7.65e-9 |
|---|
| 198 | qextrhor(1)= (3./4.)*Qext(1) / (rho_dust*radius(1)) |
|---|
| 199 | rho_q(1)=rho_dust |
|---|
| 200 | |
|---|
| 201 | endif |
|---|
| 202 | end if ! (doubleq) |
|---|
| 203 | |
|---|
| 204 | c Initialization for photochemistry: |
|---|
| 205 | c --------------------------------- |
|---|
| 206 | if (photochem) then |
|---|
| 207 | do iq=nqchem_min,nqmx |
|---|
| 208 | radius(iq)=0. |
|---|
| 209 | Qext(iq)=0. |
|---|
| 210 | alpha_lift(iq) =0. |
|---|
| 211 | alpha_devil(iq)=0. |
|---|
| 212 | qextrhor(iq)= 0. |
|---|
| 213 | enddo |
|---|
| 214 | endif |
|---|
| 215 | |
|---|
| 216 | c Initialization for water vapor |
|---|
| 217 | c ------------------------------ |
|---|
| 218 | if(water) then |
|---|
| 219 | radius(nqmx)=0. |
|---|
| 220 | Qext(nqmx)=0. |
|---|
| 221 | alpha_lift(nqmx) =0. |
|---|
| 222 | alpha_devil(nqmx)=0. |
|---|
| 223 | qextrhor(nqmx)= 0. |
|---|
| 224 | |
|---|
| 225 | c "Dryness coefficient" controlling the evaporation and |
|---|
| 226 | c sublimation from the ground water ice (close to 1) |
|---|
| 227 | c HERE, the goal is to correct for the fact |
|---|
| 228 | c that the simulated permanent water ice polar caps |
|---|
| 229 | c is larger than the actual cap and the atmospheric |
|---|
| 230 | c opacity not always realistic. |
|---|
| 231 | |
|---|
| 232 | do ig=1,ngridmx |
|---|
| 233 | if (ngridmx.ne.1) watercaptag(ig)=.false. |
|---|
| 234 | dryness(ig) = 1. |
|---|
| 235 | if (activice) pclc(ig)=1. |
|---|
| 236 | enddo |
|---|
| 237 | |
|---|
| 238 | IF (caps) THEN |
|---|
| 239 | c Perennial H20 north cap defined by watercaptag=true (allows surface to be |
|---|
| 240 | c hollowed by sublimation in vdifc). |
|---|
| 241 | c Cloud area fraction (pclc) is defined here. |
|---|
| 242 | do ig=1,ngridmx |
|---|
| 243 | if (lati(ig)*180./pi.gt.84) then |
|---|
| 244 | if (ngridmx.ne.1) watercaptag(ig)=.true. |
|---|
| 245 | dryness(ig) = 1. |
|---|
| 246 | if (activice)then |
|---|
| 247 | pclc(ig)=1. |
|---|
| 248 | print*,'Cloud area ratio : ',pclc(ig),' at lat ' |
|---|
| 249 | $ ,lati(ig)*180./pi |
|---|
| 250 | endif |
|---|
| 251 | c Use the following cap definition for high spatial resolution (latitudinal bin <= 5 deg) |
|---|
| 252 | c if (lati(ig)*180./pi.lt.85.and.long(ig).ge.0) then |
|---|
| 253 | c if (ngridmx.ne.1) watercaptag(ig)=.true. |
|---|
| 254 | c pclc(ig)=.3 |
|---|
| 255 | c dryness(ig) = 1. |
|---|
| 256 | c print*,'Cloud area ratio : ',pclc(ig),' at lat ' |
|---|
| 257 | c endif |
|---|
| 258 | c if (lati(ig)*180./pi.ge.85) then |
|---|
| 259 | c if (ngridmx.ne.1) watercaptag(ig)=.true. |
|---|
| 260 | c dryness(ig) = 1. |
|---|
| 261 | c pclc(ig)=.3 |
|---|
| 262 | c print*,'Cloud area ratio : ',pclc(ig),' at lat ' |
|---|
| 263 | c endif |
|---|
| 264 | endif ! (lati>80 deg) |
|---|
| 265 | end do ! (ngridmx) |
|---|
| 266 | ENDIF ! (caps) |
|---|
| 267 | |
|---|
| 268 | if(iceparty.and.nqmx.ge.2) then |
|---|
| 269 | radius(nqmx-1)=3.e-6 |
|---|
| 270 | rho_q(nqmx-1)=rho_ice |
|---|
| 271 | Qext(nqmx-1)=0. |
|---|
| 272 | alpha_lift(nqmx-1) =0. |
|---|
| 273 | alpha_devil(nqmx-1)=0. |
|---|
| 274 | if (activice) then |
|---|
| 275 | radius(nqmx-1)=rcrystal |
|---|
| 276 | Qext(nqmx-1)=qrefice |
|---|
| 277 | endif |
|---|
| 278 | qextrhor(nqmx-1)= (3./4.)*Qext(nqmx-1) |
|---|
| 279 | $ / (rho_ice*radius(nqmx-1)) |
|---|
| 280 | elseif(iceparty.and.nqmx.lt.2) then |
|---|
| 281 | write(*,*) 'nqmx is too low : nqmx=', nqmx |
|---|
| 282 | write(*,*) 'water= ',water,' iceparty= ',iceparty |
|---|
| 283 | endif |
|---|
| 284 | |
|---|
| 285 | end if ! (water) |
|---|
| 286 | |
|---|
| 287 | c Output for records: |
|---|
| 288 | c ~~~~~~~~~~~~~~~~~~ |
|---|
| 289 | write(*,*) |
|---|
| 290 | Write(*,*) '******** initracer : dust transport parameters :' |
|---|
| 291 | write(*,*) 'alpha_lift = ', alpha_lift |
|---|
| 292 | write(*,*) 'alpha_devil = ', alpha_devil |
|---|
| 293 | write(*,*) 'radius = ', radius |
|---|
| 294 | if(doubleq) then |
|---|
| 295 | write(*,*) 'reff_lift (um) = ', reff_lift |
|---|
| 296 | write(*,*) 'size distribution variance = ', varian |
|---|
| 297 | write(*,*) 'r3n_q , ref_r0 : ', r3n_q , ref_r0 |
|---|
| 298 | end if |
|---|
| 299 | write(*,*) 'Qext = ', qext |
|---|
| 300 | write(*,*) |
|---|
| 301 | |
|---|
| 302 | end |
|---|
