source: trunk/LMDZ.MARS/libf/phymars/co2snow.F @ 2347

Last change on this file since 2347 was 2347, checked in by aslmd, 4 years ago

CO2 cloud microphysics

File size: 5.0 KB
Line 
1c=======================================================================
2c     Program for simulate the impact of the CO2 snow fall on
3c     the surface infrared emission (emissivity)
4c-----------------------------------------------------------------------
5c     Algorithme:
6c     1 - Initialization
7c     2 - Compute the surface emissivity
8c-----------------------------------------------------------------------
9c     Author: F. Forget
10c-----------------------------------------------------------------------
11c     Reference paper:
12c     Francois Forget and James B. Pollack
13c
14c     "Thermal infrared observations of the condensing Martian polar
15c     caps: CO2 ice temperatures and radiative budget"
16c
17c     JOURNAL OF GEOPHYSICAL RESEARCH, VOL. 101, NO. E7, PAGES 16,
18c     865-16,879, JULY 25, 1996
19c=======================================================================
20
21      SUBROUTINE co2snow(ngrid, nlayer, ptimestep, emisref, condsub,
22     &                   pplev, pcondicea, pcondices, pfallice,
23     &                   pemisurf)
24     
25      use surfdat_h, only: iceradius, dtemisice
26      use geometry_mod, only: latitude ! grid point latitudes (rad)
27      use time_phylmdz_mod, only: daysec
28
29      IMPLICIT NONE
30
31#include "callkeys.h"
32
33
34c=======================================================================
35c     Variables
36c=======================================================================
37c     Inputs:
38c     -------
39      integer, intent(in) ::
40     &   ngrid, ! number of atmospheric columns
41     &   nlayer ! number of atmospheric layers
42                       
43      real, intent(in) ::
44     &   ptimestep,               ! timestep of the physics (s)
45     &   emisref(ngrid),          ! grd or ice emissivity without snow
46     &   pplev(ngrid,nlayer+1),   ! interlayer pressure (Pa)
47     &   pcondicea(ngrid,nlayer), ! CO2 condensation rate in layers
48                                  !   (kg/m2/s)
49     &   pcondices(ngrid),        ! CO2 condensation rate on the surface
50                                  !   (kg/m2/s)       
51     &   pfallice(ngrid,nlayer+1) ! falling CO2 ice (kg/m2/s)
52                           
53      logical, intent(in) ::
54     &   condsub(ngrid) ! true if there is CO2 condensation or
55                        ! sublimation in the column
56
57c     Output:
58c     -------
59      real, intent(out) ::
60     &   pemisurf(ngrid) ! surface emissivity
61
62c     local:
63c     ------
64      integer ::
65     &   l,
66     &   ig,
67     &   icap
68
69      real ::
70     &   sumdaer,
71     &   zdemisurf
72
73      real, parameter ::
74     &   alpha = 0.45
75
76c     saved:
77c     ------
78      real, save ::
79     &   Kscat(2) ! Kscat: coefficient for decreasing the surface
80c                 !        emissivity
81c                 !
82c                 ! Kscat = (0.001/3.) * alpha / iceradius
83c                 !          with 0.3 < alpha < 0.6
84c                 !
85c                 ! alpha set to 0.45 (coeff from emis = f (tau)) and
86c                 ! iceradius the mean radius of the scaterring
87c                 ! particles (200.e-6 < iceradius < 10.e-6)
88c                 !
89c                 ! (2) = N and S hemispheres
90
91      logical, save ::
92     &   firstcall = .true.
93c=======================================================================
94c BEGIN
95c=======================================================================
96c 1 - Initialization
97c=======================================================================
98      if (firstcall) then
99        Kscat(1) = (0.001/3.) * alpha / iceradius(1)
100        Kscat(2) = (0.001/3.) * alpha / iceradius(2)
101                   
102        firstcall = .false.
103      end if
104c=======================================================================
105c 2 - Compute the surface emissivity
106c=======================================================================
107      do ig = 1, ngrid
108        if (condsub(ig)) then
109          if (latitude(ig).lt.0.) then
110              icap = 2 ! Southern hemisphere
111          else
112              icap = 1 ! Northern Hemisphere
113          end if
114
115          ! compute zdemisurf using an integrated form for numerical
116          ! safety instead
117          zdemisurf =
118     &     (emisref(ig) - pemisurf(ig)) / (dtemisice(icap) * daysec)
119     &     +
120     &     ( emisref(ig) *
121     &              ( (pemisurf(ig)/emisref(ig))**(-3) +
122     &                3.*Kscat(icap)*pfallice(ig,1)*ptimestep )**(-1/3.)
123     &       - pemisurf(ig)
124     &      ) / ptimestep
125         
126          pemisurf(ig) = pemisurf(ig) + zdemisurf * ptimestep
127          if (pemisurf(ig).lt.0.1) then
128            write(*,*)'ds co2snow: emis < 0.1 !!!'
129            write(*,*)'ig =' , ig
130            write(*,*)'pemisurf(ig)', pemisurf(ig)
131            write(*,*)'zdemisurf*ptimestep', zdemisurf*ptimestep
132          end if
133        else ! if condsub(ig) is false
134           pemisurf(ig) = emisref(ig)
135        end if
136      end do ! ngrid
137c=======================================================================
138c END
139c=======================================================================
140      return
141      end
Note: See TracBrowser for help on using the repository browser.