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