Changeset 2341 for trunk/LMDZ.MARS/libf/phymars/co2snow.F
- Timestamp:
- Jun 9, 2020, 1:01:31 PM (4 years ago)
- File:
-
- 1 edited
-
trunk/LMDZ.MARS/libf/phymars/co2snow.F (modified) (1 diff)
Legend:
- Unmodified
- Added
- Removed
-
trunk/LMDZ.MARS/libf/phymars/co2snow.F
r1779 r2341 1 SUBROUTINE co2snow (ngrid,nlayer,ptimestep,emisref,condsub 2 & ,pplev,pcondicea,pcondices,pfallice,pemisurf) 1 c======================================================================= 2 c Program for simulate the impact of the CO2 snow fall on 3 c the surface infrared emission (emissivity) 4 c----------------------------------------------------------------------- 5 c Algorithme: 6 c 1 - Initialization 7 c 2 - Compute the surface emissivity 8 c----------------------------------------------------------------------- 9 c Author: F. Forget 10 c----------------------------------------------------------------------- 11 c Reference paper: 12 c Francois Forget and James B. Pollack 13 c 14 c "Thermal infrared observations of the condensing Martian polar 15 c caps: CO2 ice temperatures and radiative budget" 16 c 17 c JOURNAL OF GEOPHYSICAL RESEARCH, VOL. 101, NO. E7, PAGES 16, 18 c 865-16,879, JULY 25, 1996 19 c======================================================================= 3 20 21 SUBROUTINE co2snow(ngrid, nlayer, ptimestep, emisref, condsub, 22 & pplev, pcondicea, pcondices, pfallice, 23 & pemisurf) 24 4 25 use surfdat_h, only: iceradius, dtemisice 5 26 use geometry_mod, only: latitude ! grid point latitudes (rad) 6 27 use time_phylmdz_mod, only: daysec 28 7 29 IMPLICIT NONE 8 9 c=======================================================================10 c Program for simulate the impact of the CO2 snow fall on11 c the surface infrared emission (emissivity) and on12 c the airborne dust13 c F.Forget 199614 c=======================================================================15 16 c Declarations17 c ------------18 30 19 31 #include "callkeys.h" 20 32 21 c Arguments22 c ---------23 33 24 INTEGER,INTENT(IN) :: ngrid ! number of atmospheric columns 25 INTEGER,INTENT(IN) :: nlayer ! number of atmospheric layers 26 REAL,INTENT(IN) :: ptimestep ! timestep of the physics (s) 27 REAL,INTENT(IN) :: emisref(ngrid) ! grd or ice emissivity without snow 28 logical,intent(in) :: condsub(ngrid) ! true if there is CO2 condensation 29 ! or sublimation in the column 30 REAL,INTENT(IN) :: pplev(ngrid,nlayer+1) ! interlayer pressure (Pa) 31 REAL,INTENT(IN) :: pcondicea(ngrid,nlayer) ! CO2 condensation rate (kg/m2/s) 32 REAL,INTENT(IN) :: pcondices(ngrid) ! CO2 condensation rate (kg/m2/s) 33 ! on the surface 34 REAL,INTENT(IN) :: pfallice(ngrid,nlayer+1) ! falling CO2 ice (kg/m2/s) 34 c======================================================================= 35 c Variables 36 c======================================================================= 37 c Inputs: 38 c ------- 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 35 56 36 REAL,INTENT(OUT) :: pemisurf(ngrid) ! surface emissivity 57 c Output: 58 c ------- 59 real, intent(out) :: 60 & pemisurf(ngrid) ! surface emissivity 37 61 38 c local 39 c ----- 40 integer ig , l , icap 62 c local: 63 c ------ 64 integer :: 65 & l, 66 & ig, 67 & icap 41 68 42 REAL zdemisurf ,dtemis 43 REAL sumdaer 69 real :: 70 & sumdaer, 71 & zdemisurf 44 72 45 c saved 46 c ----- 47 REAL,save :: Kscat(2), scaveng 48 LOGICAL,SAVE :: firstcall=.true. 73 real, parameter :: 74 & alpha = 0.45 49 75 50 c -------------- 51 c Initialisation 52 c -------------- 76 c saved: 77 c ------ 78 real, save :: 79 & Kscat(2) ! Kscat: coefficient for decreasing the surface 80 c ! emissivity 81 c ! 82 c ! Kscat = (0.001/3.) * alpha / iceradius 83 c ! with 0.3 < alpha < 0.6 84 c ! 85 c ! alpha set to 0.45 (coeff from emis = f (tau)) and 86 c ! iceradius the mean radius of the scaterring 87 c ! particles (200.e-6 < iceradius < 10.e-6) 88 c ! 89 c ! (2) = N and S hemispheres 53 90 54 ! AS: firstcall OK absolute 91 logical, save :: 92 & firstcall = .true. 93 c======================================================================= 94 c BEGIN 95 c======================================================================= 96 c 1 - Initialization 97 c======================================================================= 55 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 104 c======================================================================= 105 c 2 - Compute the surface emissivity 106 c======================================================================= 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 56 114 57 c Kscat : coefficient for decreasing the surface emissivity 58 c =(0.001/3.)*alpha/iceradius , 59 c with 0.3< alpha < 0.6, set to 0.45 (coeff from emis = f (tau)) 60 c and iceradius the mean radius of the 61 c scaterring particles (200.e-6<iceradius<10.e-6) 62 63 Kscat(1)=(0.001/3.)*0.45/iceradius(1) 64 Kscat(2)=(0.001/3.)*0.45/iceradius(2) 65 66 c Scavenging Ratio (dust concentration in the air / in the snow) 67 scaveng = 100.0 68 69 c Collision Scavenging coefficient (m2.kg-1) 70 c Csca = 2.3 ! not used yet !!!!!!!!!!! 71 firstcall = .false. 72 73 end if 74 75 76 c LOOP on grid points 77 c ------------------- 78 do ig=1,ngrid 79 if (condsub(ig)) then 80 81 c IF (scavenging) then 82 c Airborne Dust 83 c ------------- 84 c sumdaer=0. 85 c do l=nlayer, 1, -1 86 c pdaerosol(ig,l)= -paerosol(ig,l,1)* 87 c & (1-exp(-scaveng*pcondicea(ig,l)*ptimestep*g/ 88 c & (pplev(ig,l)-pplev(ig,l+1))))/ptimestep 89 90 c & - Csca*paerosol(ig,l,1) ! Scavenging by collision 91 c & * 0.5*(pfallice(ig,l)) ! not included 92 93 c test to avoid releasing to much dust when subliming: 94 c if(pdaerosol(ig,l).gt.-sumdaer)pdaerosol(ig,l)=-sumdaer 95 c sumdaer=sumdaer + pdaerosol(ig,l) 96 97 c if (-pdaerosol(ig,l)*ptimestep.gt.paerosol(ig,l,1)) then 98 c write(*,*) 'ds co2snow: aerosol < 0.0 !!!' 99 c write(*,*) 'ig =' , ig 100 c end if 101 c end do 102 c END IF 103 104 c Surface emissivity 105 c ------------------ 106 c dtemis: Time scale for increasing the ice emissivity 107 108 IF(latitude(ig).LT. 0.) THEN 109 icap=2 ! Southern hemisphere 110 ELSE 111 icap=1 ! Northern Hemisphere 112 ENDIF 113 114 zdemisurf = 115 & (emisref(ig)-pemisurf(ig))/(dtemisice(icap)*daysec) 116 c Using directly the diferential equation: 117 c & -Kscat(icap)*emisref(ig)* 118 c & (pemisurf(ig)/emisref(ig))**4 *pfallice(ig,1) 119 c Using an integrated form for numerical safety instead 120 & +(emisref(ig)* ((pemisurf(ig)/emisref(ig))**(-3)+3.*Kscat(icap)* 121 & pfallice(ig,1)*ptimestep)**(-1/3.) -pemisurf(ig))/ptimestep 122 123 124 pemisurf(ig) = pemisurf(ig) + zdemisurf*ptimestep 125 126 if (pemisurf(ig).lt.0.1) then 127 write(*,*) 'ds co2snow: emis < 0.1 !!!' 128 write(*,*) 'ig =' , ig 129 write(*,*)'pemisurf(ig)',pemisurf(ig) 130 write(*,*) 'zdemisurf*ptimestep',zdemisurf*ptimestep 131 end if 132 else 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 133 134 pemisurf(ig) = emisref(ig) 134 end if 135 end do 136 135 end if 136 end do ! ngrid 137 c======================================================================= 138 c END 139 c======================================================================= 137 140 return 138 141 end
Note: See TracChangeset
for help on using the changeset viewer.
