| 1 | SUBROUTINE initracer() |
|---|
| 2 | |
|---|
| 3 | |
|---|
| 4 | IMPLICIT NONE |
|---|
| 5 | c======================================================================= |
|---|
| 6 | c subject: |
|---|
| 7 | c -------- |
|---|
| 8 | c Initialization related to tracer |
|---|
| 9 | c (transported dust, water, chemical species, ice...) |
|---|
| 10 | c |
|---|
| 11 | c Name of the tracer |
|---|
| 12 | c |
|---|
| 13 | c Test of dimension : |
|---|
| 14 | c Initialize COMMON tracer in tracer.h, using tracer names provided |
|---|
| 15 | c by the dynamics in "advtrac.h" |
|---|
| 16 | c |
|---|
| 17 | c author: F.Forget |
|---|
| 18 | c ------ |
|---|
| 19 | c Ehouarn Millour (oct. 2008) identify tracers by their names |
|---|
| 20 | c======================================================================= |
|---|
| 21 | |
|---|
| 22 | |
|---|
| 23 | #include "dimensions.h" |
|---|
| 24 | #include "dimphys.h" |
|---|
| 25 | #include "comcstfi.h" |
|---|
| 26 | #include "callkeys.h" |
|---|
| 27 | #include "tracer.h" |
|---|
| 28 | #include "advtrac.h" |
|---|
| 29 | #include "comgeomfi.h" |
|---|
| 30 | #include "watercap.h" |
|---|
| 31 | #include "fisice.h" |
|---|
| 32 | |
|---|
| 33 | ! real qsurf(ngridmx,nqmx) ! tracer on surface (e.g. kg.m-2) |
|---|
| 34 | ! real co2ice(ngridmx) ! co2 ice mass on surface (e.g. kg.m-2) |
|---|
| 35 | character(len=20) :: txt ! to store some text |
|---|
| 36 | integer iq,ig,count |
|---|
| 37 | real r0_lift , reff_lift |
|---|
| 38 | ! logical :: oldnames ! =.true. if old tracer naming convention (q01,...) |
|---|
| 39 | |
|---|
| 40 | |
|---|
| 41 | c----------------------------------------------------------------------- |
|---|
| 42 | c radius(nqmx) ! aerosol particle radius (m) |
|---|
| 43 | c rho_q(nqmx) ! tracer densities (kg.m-3) |
|---|
| 44 | c qext(nqmx) ! Single Scat. Extinction coeff at 0.67 um |
|---|
| 45 | c alpha_lift(nqmx) ! saltation vertical flux/horiz flux ratio (m-1) |
|---|
| 46 | c alpha_devil(nqmx) ! lifting coeeficient by dust devil |
|---|
| 47 | c rho_dust ! Mars dust density |
|---|
| 48 | c rho_ice ! Water ice density |
|---|
| 49 | c doubleq ! if method with mass (iq=1) and number(iq=2) mixing ratio |
|---|
| 50 | c varian ! Characteristic variance of log-normal distribution |
|---|
| 51 | c----------------------------------------------------------------------- |
|---|
| 52 | |
|---|
| 53 | ! Initialization: get tracer names from the dynamics and check if we are |
|---|
| 54 | ! using 'old' tracer convention ('q01',q02',...) |
|---|
| 55 | ! or new convention (full tracer names) |
|---|
| 56 | ! check if tracers have 'old' names |
|---|
| 57 | |
|---|
| 58 | ! copy tracer names from dynamics |
|---|
| 59 | do iq=1,nqmx |
|---|
| 60 | noms(iq)=tnom(iq) |
|---|
| 61 | enddo |
|---|
| 62 | |
|---|
| 63 | |
|---|
| 64 | print*,'nqmx=',nqmx |
|---|
| 65 | ! Identify tracers by their names: (and set corresponding values of mmol) |
|---|
| 66 | ! 0. initialize tracer indexes to zero: |
|---|
| 67 | ! NB: igcm_* indexes are commons in 'tracer.h' |
|---|
| 68 | do iq=1,nqmx |
|---|
| 69 | igcm_dustbin(iq)=0 |
|---|
| 70 | enddo |
|---|
| 71 | igcm_dust_mass=0 |
|---|
| 72 | igcm_dust_number=0 |
|---|
| 73 | igcm_h2o_vap=0 |
|---|
| 74 | igcm_h2o_ice=0 |
|---|
| 75 | igcm_co2=0 |
|---|
| 76 | igcm_co=0 |
|---|
| 77 | igcm_o=0 |
|---|
| 78 | igcm_o1d=0 |
|---|
| 79 | igcm_o2=0 |
|---|
| 80 | igcm_o3=0 |
|---|
| 81 | igcm_h=0 |
|---|
| 82 | igcm_h2=0 |
|---|
| 83 | igcm_oh=0 |
|---|
| 84 | igcm_ho2=0 |
|---|
| 85 | igcm_h2o2=0 |
|---|
| 86 | igcm_n2=0 |
|---|
| 87 | igcm_ar=0 |
|---|
| 88 | igcm_ar_n2=0 |
|---|
| 89 | igcm_co2_ice=0 |
|---|
| 90 | |
|---|
| 91 | write(*,*) 'initracer: noms() ', noms |
|---|
| 92 | |
|---|
| 93 | |
|---|
| 94 | print*,'Setting dustbin = 0 in initracer.F' |
|---|
| 95 | dustbin=0 |
|---|
| 96 | |
|---|
| 97 | ! 1. find dust tracers |
|---|
| 98 | count=0 |
|---|
| 99 | ! if (dustbin.gt.0) then |
|---|
| 100 | ! do iq=1,nqmx |
|---|
| 101 | ! txt=" " |
|---|
| 102 | ! write(txt,'(a4,i2.2)')'dust',count+1 |
|---|
| 103 | ! if (noms(iq).eq.txt) then |
|---|
| 104 | ! count=count+1 |
|---|
| 105 | ! igcm_dustbin(count)=iq |
|---|
| 106 | ! mmol(iq)=100. |
|---|
| 107 | ! endif |
|---|
| 108 | ! enddo !do iq=1,nqmx |
|---|
| 109 | ! endif ! of if (dustbin.gt.0) |
|---|
| 110 | |
|---|
| 111 | |
|---|
| 112 | ! if (doubleq) then |
|---|
| 113 | ! do iq=1,nqmx |
|---|
| 114 | ! if (noms(iq).eq."dust_mass") then |
|---|
| 115 | ! igcm_dust_mass=iq |
|---|
| 116 | ! count=count+1 |
|---|
| 117 | ! endif |
|---|
| 118 | ! if (noms(iq).eq."dust_number") then |
|---|
| 119 | ! igcm_dust_number=iq |
|---|
| 120 | ! count=count+1 |
|---|
| 121 | ! endif |
|---|
| 122 | ! enddo |
|---|
| 123 | ! endif ! of if (doubleq) |
|---|
| 124 | ! 2. find chemistry and water tracers |
|---|
| 125 | do iq=1,nqmx |
|---|
| 126 | if (noms(iq).eq."co2") then |
|---|
| 127 | igcm_co2=iq |
|---|
| 128 | mmol(igcm_co2)=44. |
|---|
| 129 | count=count+1 |
|---|
| 130 | ! write(*,*) 'co2: count=',count |
|---|
| 131 | endif |
|---|
| 132 | if (noms(iq).eq."co2_ice") then |
|---|
| 133 | igcm_co2_ice=iq |
|---|
| 134 | mmol(igcm_co2_ice)=44. |
|---|
| 135 | count=count+1 |
|---|
| 136 | ! write(*,*) 'co2_ice: count=',count |
|---|
| 137 | endif |
|---|
| 138 | if (noms(iq).eq."h2o_vap") then |
|---|
| 139 | igcm_h2o_vap=iq |
|---|
| 140 | mmol(igcm_h2o_vap)=18. |
|---|
| 141 | count=count+1 |
|---|
| 142 | ! write(*,*) 'h2o_vap: count=',count |
|---|
| 143 | endif |
|---|
| 144 | if (noms(iq).eq."h2o_ice") then |
|---|
| 145 | igcm_h2o_ice=iq |
|---|
| 146 | mmol(igcm_h2o_ice)=18. |
|---|
| 147 | count=count+1 |
|---|
| 148 | ! write(*,*) 'h2o_ice: count=',count |
|---|
| 149 | endif |
|---|
| 150 | enddo ! of do iq=1,nqmx |
|---|
| 151 | |
|---|
| 152 | ! check that we identified all tracers: |
|---|
| 153 | if (count.ne.nqmx) then |
|---|
| 154 | write(*,*) "initracer: found only ",count," tracers" |
|---|
| 155 | write(*,*) " expected ",nqmx |
|---|
| 156 | do iq=1,count |
|---|
| 157 | write(*,*)' ',iq,' ',trim(noms(iq)) |
|---|
| 158 | enddo |
|---|
| 159 | stop |
|---|
| 160 | else |
|---|
| 161 | write(*,*) "initracer: found all expected tracers, namely:" |
|---|
| 162 | do iq=1,nqmx |
|---|
| 163 | write(*,*)' ',iq,' ',trim(noms(iq)) |
|---|
| 164 | enddo |
|---|
| 165 | endif |
|---|
| 166 | |
|---|
| 167 | |
|---|
| 168 | c------------------------------------------------------------ |
|---|
| 169 | c Initialisation tracers .... |
|---|
| 170 | c------------------------------------------------------------ |
|---|
| 171 | call zerophys(nqmx,rho_q) |
|---|
| 172 | |
|---|
| 173 | rho_dust=2500. ! Mars dust density (kg.m-3) |
|---|
| 174 | rho_ice=920. ! Water ice density (kg.m-3) |
|---|
| 175 | rho_co2=1620. ! CO2 ice density (kg.m-3) |
|---|
| 176 | |
|---|
| 177 | |
|---|
| 178 | |
|---|
| 179 | c$$$ if (doubleq) then |
|---|
| 180 | c$$$c "doubleq" technique |
|---|
| 181 | c$$$c ------------------- |
|---|
| 182 | c$$$c (transport of mass and number mixing ratio) |
|---|
| 183 | c$$$c iq=1: Q mass mixing ratio, iq=2: N number mixing ratio |
|---|
| 184 | c$$$ |
|---|
| 185 | c$$$ if( (nqmx.lt.2).or.(water.and.(nqmx.lt.3)) ) then |
|---|
| 186 | c$$$ write(*,*)'initracer: nqmx is too low : nqmx=', nqmx |
|---|
| 187 | c$$$ write(*,*)'water= ',water,' doubleq= ',doubleq |
|---|
| 188 | c$$$ end if |
|---|
| 189 | c$$$ |
|---|
| 190 | c$$$ varian=0.637 ! Characteristic variance |
|---|
| 191 | c$$$ qext(igcm_dust_mass)=3.04 ! reference extinction at 0.67 um for ref dust |
|---|
| 192 | c$$$ qext(igcm_dust_number)=3.04 ! reference extinction at 0.67 um for ref dust |
|---|
| 193 | c$$$ rho_q(igcm_dust_mass)=rho_dust |
|---|
| 194 | c$$$ rho_q(igcm_dust_number)=rho_dust |
|---|
| 195 | c$$$ |
|---|
| 196 | c$$$c Intermediate calcul for computing geometric mean radius r0 |
|---|
| 197 | c$$$c as a function of mass and number mixing ratio Q and N |
|---|
| 198 | c$$$c (r0 = (r3n_q * Q/ N) |
|---|
| 199 | c$$$ r3n_q = exp(-4.5*varian**2)*(3./4.)/(pi*rho_dust) |
|---|
| 200 | c$$$ |
|---|
| 201 | c$$$c Intermediate calcul for computing effective radius reff |
|---|
| 202 | c$$$c from geometric mean radius r0 |
|---|
| 203 | c$$$c (reff = ref_r0 * r0) |
|---|
| 204 | c$$$ ref_r0 = exp(2.5*varian**2) |
|---|
| 205 | c$$$ |
|---|
| 206 | c$$$c lifted dust : |
|---|
| 207 | c$$$c ''''''''''' |
|---|
| 208 | c$$$ reff_lift = 3.e-6 ! Effective radius of lifted dust (m) |
|---|
| 209 | c$$$ alpha_devil(igcm_dust_mass)=9.e-9 ! dust devil lift mass coeff |
|---|
| 210 | c$$$ alpha_lift(igcm_dust_mass)=3.0e-15 ! Lifted mass coeff |
|---|
| 211 | c$$$ |
|---|
| 212 | c$$$ r0_lift = reff_lift/ref_r0 |
|---|
| 213 | c$$$ alpha_devil(igcm_dust_number)=r3n_q* |
|---|
| 214 | c$$$ & alpha_devil(igcm_dust_mass)/r0_lift**3 |
|---|
| 215 | c$$$ alpha_lift(igcm_dust_number)=r3n_q* |
|---|
| 216 | c$$$ & alpha_lift(igcm_dust_mass)/r0_lift**3 |
|---|
| 217 | c$$$ |
|---|
| 218 | c$$$c Not used: |
|---|
| 219 | c$$$ radius(igcm_dust_mass) = 0. |
|---|
| 220 | c$$$ radius(igcm_dust_number) = 0. |
|---|
| 221 | c$$$ |
|---|
| 222 | c$$$ else |
|---|
| 223 | |
|---|
| 224 | print*,dustbin |
|---|
| 225 | |
|---|
| 226 | if (dustbin.gt.1) then |
|---|
| 227 | print*,'ATTENTION:', |
|---|
| 228 | $ ' properties of dust need input in initracer !!!' |
|---|
| 229 | stop |
|---|
| 230 | |
|---|
| 231 | else if (dustbin.eq.1) then |
|---|
| 232 | |
|---|
| 233 | c This will be used for 1 dust particle size: |
|---|
| 234 | c ------------------------------------------ |
|---|
| 235 | radius(igcm_dustbin(1))=3.e-6 |
|---|
| 236 | Qext(igcm_dustbin(1))=3.04 |
|---|
| 237 | alpha_lift(igcm_dustbin(1))=0.0e-6 |
|---|
| 238 | alpha_devil(igcm_dustbin(1))=7.65e-9 |
|---|
| 239 | qextrhor(igcm_dustbin(1))=(3./4.)*Qext(igcm_dustbin(1)) |
|---|
| 240 | & /(rho_dust*radius(igcm_dustbin(1))) |
|---|
| 241 | rho_q(igcm_dustbin(1))=rho_dust |
|---|
| 242 | |
|---|
| 243 | endif |
|---|
| 244 | ! end if ! (doubleq) |
|---|
| 245 | |
|---|
| 246 | c Initialization for water vapor |
|---|
| 247 | c ------------------------------ |
|---|
| 248 | if(water) then |
|---|
| 249 | radius(igcm_h2o_vap)=0. |
|---|
| 250 | Qext(igcm_h2o_vap)=0. |
|---|
| 251 | alpha_lift(igcm_h2o_vap) =0. |
|---|
| 252 | alpha_devil(igcm_h2o_vap)=0. |
|---|
| 253 | qextrhor(igcm_h2o_vap)= 0. |
|---|
| 254 | |
|---|
| 255 | c "Dryness coefficient" controlling the evaporation and |
|---|
| 256 | c sublimation from the ground water ice (close to 1) |
|---|
| 257 | c HERE, the goal is to correct for the fact |
|---|
| 258 | c that the simulated permanent water ice polar caps |
|---|
| 259 | c is larger than the actual cap and the atmospheric |
|---|
| 260 | c opacity not always realistic. |
|---|
| 261 | |
|---|
| 262 | do ig=1,ngridmx |
|---|
| 263 | if (ngridmx.ne.1) watercaptag(ig)=.false. |
|---|
| 264 | dryness(ig) = 1. |
|---|
| 265 | enddo |
|---|
| 266 | |
|---|
| 267 | ! IF (caps) THEN |
|---|
| 268 | c Perennial H20 north cap defined by watercaptag=true (allows surface to be |
|---|
| 269 | c hollowed by sublimation in vdifc). |
|---|
| 270 | ! do ig=1,ngridmx |
|---|
| 271 | ! if (lati(ig)*180./pi.gt.84) then |
|---|
| 272 | ! if (ngridmx.ne.1) watercaptag(ig)=.true. |
|---|
| 273 | ! dryness(ig) = 1. |
|---|
| 274 | c Use the following cap definition for high spatial resolution (latitudinal bin <= 5 deg) |
|---|
| 275 | c if (lati(ig)*180./pi.lt.85.and.long(ig).ge.0) then |
|---|
| 276 | c if (ngridmx.ne.1) watercaptag(ig)=.true. |
|---|
| 277 | c dryness(ig) = 1. |
|---|
| 278 | c endif |
|---|
| 279 | c if (lati(ig)*180./pi.ge.85) then |
|---|
| 280 | c if (ngridmx.ne.1) watercaptag(ig)=.true. |
|---|
| 281 | c dryness(ig) = 1. |
|---|
| 282 | c endif |
|---|
| 283 | ! endif ! (lati>80 deg) |
|---|
| 284 | ! end do ! (ngridmx) |
|---|
| 285 | ! ENDIF ! (caps) |
|---|
| 286 | |
|---|
| 287 | ! if(iceparty.and.(nqmx.ge.2)) then |
|---|
| 288 | |
|---|
| 289 | radius(igcm_h2o_ice)=3.e-6 |
|---|
| 290 | rho_q(igcm_h2o_ice)=rho_ice |
|---|
| 291 | Qext(igcm_h2o_ice)=0. |
|---|
| 292 | ! alpha_lift(igcm_h2o_ice) =0. |
|---|
| 293 | ! alpha_devil(igcm_h2o_ice)=0. |
|---|
| 294 | qextrhor(igcm_h2o_ice)= (3./4.)*Qext(igcm_h2o_ice) |
|---|
| 295 | $ / (rho_ice*radius(igcm_h2o_ice)) |
|---|
| 296 | |
|---|
| 297 | |
|---|
| 298 | |
|---|
| 299 | ! elseif(iceparty.and.(nqmx.lt.2)) then |
|---|
| 300 | ! write(*,*) 'nqmx is too low : nqmx=', nqmx |
|---|
| 301 | ! write(*,*) 'water= ',water,' iceparty= ',iceparty |
|---|
| 302 | ! endif |
|---|
| 303 | |
|---|
| 304 | end if ! (water) |
|---|
| 305 | |
|---|
| 306 | c Output for records: |
|---|
| 307 | c ~~~~~~~~~~~~~~~~~~ |
|---|
| 308 | write(*,*) |
|---|
| 309 | Write(*,*) '******** initracer : dust transport parameters :' |
|---|
| 310 | write(*,*) 'alpha_lift = ', alpha_lift |
|---|
| 311 | write(*,*) 'alpha_devil = ', alpha_devil |
|---|
| 312 | write(*,*) 'radius = ', radius |
|---|
| 313 | ! if(doubleq) then |
|---|
| 314 | ! write(*,*) 'reff_lift (um) = ', reff_lift |
|---|
| 315 | ! write(*,*) 'size distribution variance = ', varian |
|---|
| 316 | ! write(*,*) 'r3n_q , ref_r0 : ', r3n_q , ref_r0 |
|---|
| 317 | ! end if |
|---|
| 318 | write(*,*) 'Qext = ', qext |
|---|
| 319 | write(*,*) |
|---|
| 320 | |
|---|
| 321 | end |
|---|
