source: trunk/LMDZ.GENERIC/libf/phystd/initracer.F @ 586

Last change on this file since 586 was 253, checked in by emillour, 13 years ago

Generic GCM

  • Massive update to version 0.7

EM+RW

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