source: trunk/WRF.COMMON/WRFV3/phys/module_gocart_seasalt.F

Last change on this file was 2759, checked in by aslmd, 2 years ago

adding unmodified code from WRFV3.0.1.1, expurged from useless data +1M size
File size: 9.8 KB
Line 
1MODULE GOCART_SEASALT
2
3CONTAINS
4  subroutine gocart_seasalt_driver(ktau,dt,config_flags,julday,alt,t_phy,moist,u_phy,  &
5         v_phy,chem,rho_phy,dz8w,u10,v10,p8w,                  &
6         xland,xlat,xlong,dx,g,emis_seas, &
7         ids,ide, jds,jde, kds,kde,                                        &
8         ims,ime, jms,jme, kms,kme,                                        &
9         its,ite, jts,jte, kts,kte                                         )
10  USE module_configure
11  USE module_state_description
12  USE module_model_constants, ONLY: mwdry
13  IMPLICIT NONE
14   TYPE(grid_config_rec_type),  INTENT(IN   )    :: config_flags
15
16   INTEGER,      INTENT(IN   ) :: julday, ktau,                     &
17                                  ids,ide, jds,jde, kds,kde,               &
18                                  ims,ime, jms,jme, kms,kme,               &
19                                  its,ite, jts,jte, kts,kte
20   REAL, DIMENSION( ims:ime, kms:kme, jms:jme, num_moist ),                &
21         INTENT(IN ) ::                                   moist
22   REAL, DIMENSION( ims:ime, kms:kme, jms:jme, num_chem ),                 &
23         INTENT(INOUT ) ::                                   chem
24   REAL, DIMENSION( ims:ime, 1, jms:jme,num_emis_seas),OPTIONAL,&
25         INTENT(INOUT ) ::                                                 &
26         emis_seas
27   REAL,  DIMENSION( ims:ime , jms:jme )                   ,               &
28          INTENT(IN   ) ::                                                 &
29                                                     u10,                  &
30                                                     v10,                  &
31                                                     xland,                &
32                                                     xlat,                 &
33                                                     xlong
34   REAL,  DIMENSION( ims:ime , kms:kme , jms:jme ),                        &
35          INTENT(IN   ) ::                                                 &
36                                                        alt,               &
37                                                      t_phy,               &
38                                                     dz8w,p8w,             &
39                                              u_phy,v_phy,rho_phy
40
41  REAL, INTENT(IN   ) :: dt,dx,g
42!
43! local variables
44!
45  integer :: ipr,nmx,i,j,k,ndt,imx,jmx,lmx
46  integer,dimension (1,1) :: ilwi
47  real*8, DIMENSION (4) :: tc,bems
48  real*8, dimension (1,1) :: w10m,gwet,airden,airmas
49  real*8, dimension (1) :: dxy
50  real*8 conver,converi
51  conver=1.d-9
52  converi=1.d9
53!
54! number of dust bins
55!
56  imx=1
57  jmx=1
58  lmx=1
59  nmx=4
60  k=kts
61  do j=jts,jte
62  do i=its,ite
63!
64! donṫ do dust over water!!!
65!
66     if(xland(i,j).gt.1.5)then
67     ilwi(1,1)=0
68     tc(1)=chem(i,kts,j,p_seas_1)*conver
69     tc(2)=chem(i,kts,j,p_seas_2)*conver
70     tc(3)=chem(i,kts,j,p_seas_3)*conver
71     tc(4)=chem(i,kts,j,p_seas_4)*conver
72     w10m(1,1)=sqrt(u10(i,j)*u10(i,j)+v10(i,j)*v10(i,j))
73     airmas(1,1)=-(p8w(i,kts+1,j)-p8w(i,kts,j))*dx*dx/g
74!
75! we donṫ trust the u10,v10 values, is model layers are very thin near surface
76!
77     if(dz8w(i,kts,j).lt.12.)w10m=sqrt(u_phy(i,kts,j)*u_phy(i,kts,j)+v_phy(i,kts,j)*v_phy(i,kts,j))
78!
79     dxy(1)=dx*dx
80       ipr=0
81
82    call source_ss( imx,jmx,lmx,nmx, dt, tc,ilwi, dxy, w10m, airmas, bems,ipr)
83     chem(i,kts,j,p_seas_1)=tc(1)*converi
84     chem(i,kts,j,p_seas_2)=tc(2)*converi
85     chem(i,kts,j,p_seas_3)=tc(3)*converi
86     chem(i,kts,j,p_seas_4)=tc(4)*converi
87! for output diagnostics
88     emis_seas(i,1,j,p_edust1)=bems(1)
89     emis_seas(i,1,j,p_edust2)=bems(2)
90     emis_seas(i,1,j,p_edust3)=bems(3)
91     emis_seas(i,1,j,p_edust4)=bems(4)
92     endif
93  enddo
94  enddo
95!
96
97end subroutine gocart_seasalt_driver
98!
99SUBROUTINE source_ss(imx,jmx,lmx,nmx, dt1, tc, &
100                     ilwi, dxy, w10m, airmas, &
101                     bems,ipr)
102
103! ****************************************************************************
104! *  Evaluate the source of each seasalt particles size classes  (kg/m3)
105! *  by soil emission.
106! *  Input:
107! *         SSALTDEN  Sea salt density                               (kg/m3)
108! *         DXY       Surface of each grid cell                     (m2)
109! *         NDT1      Time step                                     (s)
110! *         W10m      Velocity at the anemometer level (10meters)   (m/s)
111! *     
112! *  Output:
113! *         DSRC      Source of each sea salt bins       (kg/timestep/cell)
114! *
115! *
116! * Number flux density: Original formula by Monahan et al. (1986) adapted
117! * by Sunling Gong (JGR 1997 (old) and GBC 2003 (new)).  The new version is
118! * to better represent emission of sub-micron sea salt particles.
119!
120! * dFn/dr = c1*u10**c2/(r**A) * (1+c3*r**c4)*10**(c5*exp(-B**2))
121! * where B = (b1 -log(r))/b2
122! * see c_old, c_new, b_old, b_new below for the constants.
123! * number fluxes are at 80% RH.
124! *
125! * To calculate the flux:
126! * 1) Calculate dFn based on Monahan et al. (1986) and Gong (2003)
127! * 2) Assume that wet radius r at 80% RH = dry radius r_d *frh
128! * 3) Convert particles flux to mass flux :
129! *    dFM/dr_d = 4/3*pi*rho_d*r_d^3 *(dr/dr_d) * dFn/dr
130! *             = 4/3*pi*rho_d*r_d^3 * frh * dFn/dr
131! *               where rho_p is particle density [kg/m3]
132! *    The factor 1.e-18 is to convert in micro-meter r_d^3
133! ****************************************************************************
134 
135
136  USE module_data_gocart_seas
137
138  IMPLICIT NONE
139
140  INTEGER, INTENT(IN)    :: nmx,imx,jmx,lmx,ipr
141  INTEGER, INTENT(IN)    :: ilwi(imx,jmx)
142  REAL*8,    INTENT(IN)    :: dxy(jmx), w10m(imx,jmx)
143  REAL*8,    INTENT(IN)    :: airmas(imx,jmx,lmx)
144  REAL*8,    INTENT(INOUT) :: tc(imx,jmx,lmx,nmx)
145  REAL*8,    INTENT(OUT)   :: bems(imx,jmx,nmx)
146
147  REAL*8 :: c0(5), b0(2)
148!  REAL*8, PARAMETER :: c_old(5)=(/1.373, 3.41, 0.057, 1.05, 1.190/)
149!  REAL*8, PARAMETER :: c_new(5)=(/1.373, 3.41, 0.057, 3.45, 1.607/)
150  ! Change suggested by MC
151  REAL*8, PARAMETER :: c_old(5)=(/1.373, 3.2, 0.057, 1.05, 1.190/)
152  REAL*8, PARAMETER :: c_new(5)=(/1.373, 3.2, 0.057, 3.45, 1.607/)
153  REAL*8, PARAMETER :: b_old(2)=(/0.380, 0.650/)
154  REAL*8, PARAMETER :: b_new(2)=(/0.433, 0.433/)
155  REAL*8, PARAMETER :: dr=5.0D-2 ! um   
156  REAL*8, PARAMETER :: theta=30.0
157  ! Swelling coefficient frh (d rwet / d rd)
158!!!  REAL*8,    PARAMETER :: frh = 1.65
159  REAL*8,    PARAMETER :: frh = 2.d0
160  LOGICAL, PARAMETER :: old=.TRUE., new=.FALSE.
161  REAL*8 :: rho_d, r0, r1, r, r_w, a, b, dfn, r_d, dfm, src
162  INTEGER :: i, j, n, nr, ir
163  REAL :: dt1
164
165
166  REAL*8                  :: tcmw(nmx), ar(nmx), tcvv(nmx)
167  REAL*8                  :: ar_wetdep(nmx), kc(nmx)
168  CHARACTER(LEN=20)     :: tcname(nmx), tcunits(nmx)
169  LOGICAL               :: aerosol(nmx)
170
171
172  REAL*8 :: tc1(imx,jmx,lmx,nmx)
173  REAL*8, TARGET :: tcms(imx,jmx,lmx,nmx) ! tracer mass (kg; kgS for sulfur case)
174  REAL*8, TARGET :: tcgm(imx,jmx,lmx,nmx) ! g/m3
175
176  !----------------------------------------------------------------------- 
177  ! sea salt specific
178  !----------------------------------------------------------------------- 
179! REAL*8, DIMENSION(nmx) :: ra, rb
180! REAL*8 :: ch_ss(nmx,12)
181
182  !----------------------------------------------------------------------- 
183  ! emissions (input)
184  !----------------------------------------------------------------------- 
185  REAL*8 :: e_an(imx,jmx,2,nmx), e_bb(imx,jmx,nmx), &
186          e_ac(imx,jmx,lmx,nmx)
187
188  !----------------------------------------------------------------------- 
189  ! diagnostics (budget)
190  !-----------------------------------------------------------------------
191!  ! tendencies per time step and process
192!  REAL*8, TARGET :: bems(imx,jmx,nmx), bdry(imx,jmx,nmx), bstl(imx,jmx,nmx)
193!  REAL*8, TARGET :: bwet(imx,jmx,nmx), bcnv(imx,jmx,nmx)!
194
195!  ! integrated tendencies per process
196!  REAL*8, TARGET :: tems(imx,jmx,nmx), tstl(imx,jmx,nmx)
197!  REAL*8, TARGET :: tdry(imx,jmx,nmx), twet(imx,jmx,nmx), tcnv(imx,jmx,nmx)
198
199  ! global mass balance per time step
200  REAL*8 :: tmas0(nmx), tmas1(nmx)
201  REAL*8 :: dtems(nmx), dttrp(nmx), dtdif(nmx), dtcnv(nmx)
202  REAL*8 :: dtwet(nmx), dtdry(nmx), dtstl(nmx)
203  REAL*8 :: dtems2(nmx), dttrp2(nmx), dtdif2(nmx), dtcnv2(nmx)
204  REAL*8 :: dtwet2(nmx), dtdry2(nmx), dtstl2(nmx)
205
206  ! detailed integrated budgets for individual emissions
207  REAL*8, TARGET :: ems_an(imx,jmx,nmx),    ems_bb(imx,jmx,nmx), ems_tp(imx,jmx)
208  REAL*8, TARGET :: ems_ac(imx,jmx,lmx,nmx)
209  REAL*8, TARGET :: ems_co(imx,jmx,nmx)
210
211
212  ! executable statements
213
214  DO n = 1,nmx
215!    if(ipr.eq.1)write(0,*)'in seasalt',n,ipr,ilwi
216     bems(:,:,n) = 0.0
217     rho_d = den_seas(n)
218     r0 = ra(n)*frh
219     r1 = rb(n)*frh
220     r = r0
221     nr = INT((r1-r0)/dr+.001)
222!    if(ipr.eq.1.and.n.eq.1)write(0,*)'in seasalt',nr,r1,r0,dr,rho_d
223     DO ir = 1,nr
224        r_w = r + dr*0.5
225        r = r + dr
226        IF (new) THEN
227           a = 4.7*(1.0 + theta*r_w)**(-0.017*r_w**(-1.44))
228           c0 = c_new
229           b0 = b_new
230        ELSE
231           a = 3.0
232           c0 = c_old
233           b0 = b_old
234        END IF
235        !
236        b = (b0(1) - LOG10(r_w))/b0(2)
237        dfn = (c0(1)/r_w**a)*(1.0 + c0(3)*r_w**c0(4))* &
238             10**(c0(5)*EXP(-(b**2)))
239       
240        r_d = r_w/frh*1.0D-6  ! um -> m
241        dfm = 4.0/3.0*pi*r_d**3*rho_d*frh*dfn*dr*dt1
242        DO i = 1,imx
243           DO j = 1,jmx
244!              IF (water(i,j) > 0.0) THEN
245              IF (ilwi(i,j) == 0) THEN
246!                 src = dfm*dxy(j)*water(i,j)*w10m(i,j)**c0(2)
247                 src = dfm*dxy(j)*w10m(i,j)**c0(2)
248!                 src = ch_ss(n,dt(1)%mn)*dfm*dxy(j)*w10m(i,j)**c0(2)
249                 tc(i,j,1,n) = tc(i,j,1,n) + src/airmas(i,j,1)
250!                if(ipr.eq.1)write(0,*)n,dfm,c0(2),dxy(j),w10m(i,j),src,airmas(i,j,1)
251              ELSE
252                 src = 0.0
253              END IF
254              bems(i,j,n) = bems(i,j,n) + src
255           END DO  ! i
256        END DO ! j
257     END DO ! ir
258  END DO ! n
259
260END SUBROUTINE source_ss
261END MODULE GOCART_SEASALT
Note: See TracBrowser for help on using the repository browser.