| 1 | subroutine improvedCO2clouds(ngrid,nlay,ptimestep, |
|---|
| 2 | & pplay,pt,pdt, |
|---|
| 3 | & pq,pdq,pdqcloudco2,pdtcloudco2, |
|---|
| 4 | & nq,tauscaling) |
|---|
| 5 | ! to use 'getin' |
|---|
| 6 | USE comcstfi_h |
|---|
| 7 | USE ioipsl_getincom |
|---|
| 8 | USE updaterad |
|---|
| 9 | use tracer_mod |
|---|
| 10 | !, only: rho_ice_co2, nuiceco2_sed, igcm_co2, |
|---|
| 11 | ! & rho_ice,igcm_h2o_ice, igcm_ccn_number, |
|---|
| 12 | ! & igcm_co2_ice, igcm_dust_mass, |
|---|
| 13 | ! & igcm_dust_number, igcm_ccnco2_mass, |
|---|
| 14 | ! & igcm_ccnco2_number |
|---|
| 15 | use conc_mod, only: mmean |
|---|
| 16 | implicit none |
|---|
| 17 | |
|---|
| 18 | |
|---|
| 19 | c------------------------------------------------------------------ |
|---|
| 20 | c This routine is used to form CO2 clouds when a parcel of the GCM is |
|---|
| 21 | c saturated. It includes the ability to have supersaturation, a |
|---|
| 22 | c computation of the nucleation rates, growthrates and the |
|---|
| 23 | c scavenging of dust particles by clouds. |
|---|
| 24 | c It is worth noting that the amount of dust is computed using the |
|---|
| 25 | c dust optical depth computed in aeropacity.F. That's why |
|---|
| 26 | c the variable called "tauscaling" is used to convert |
|---|
| 27 | c pq(dust_mass) and pq(dust_number), which are relative |
|---|
| 28 | c quantities, to absolute and realistic quantities stored in zq. |
|---|
| 29 | c This has to be done to convert the inputs into absolute |
|---|
| 30 | c values, but also to convert the outputs back into relative |
|---|
| 31 | c values which are then used by the sedimentation and advection |
|---|
| 32 | c schemes. |
|---|
| 33 | |
|---|
| 34 | c Authors of the water ice clouds microphysics |
|---|
| 35 | c J.-B. Madeleine, based on the work by Franck Montmessin |
|---|
| 36 | c (October 2011) |
|---|
| 37 | c T. Navarro, debug,correction, new scheme (October-April 2011) |
|---|
| 38 | c A. Spiga, optimization (February 2012) |
|---|
| 39 | c Adaptation for CO2 clouds by Joachim Audouard (09/16), based on the work |
|---|
| 40 | c of Constantino Listowski |
|---|
| 41 | c------------------------------------------------------------------ |
|---|
| 42 | !#include "dimensions.h" |
|---|
| 43 | !#include "dimphys.h" |
|---|
| 44 | #include "callkeys.h" |
|---|
| 45 | !#include "tracer.h" |
|---|
| 46 | !#include "comgeomfi.h" |
|---|
| 47 | !#include "dimradmars.h" |
|---|
| 48 | #include "microphys.h" |
|---|
| 49 | !#include "microphysCO2.h" |
|---|
| 50 | !#include "conc.h" |
|---|
| 51 | c------------------------------------------------------------------ |
|---|
| 52 | c Inputs: |
|---|
| 53 | |
|---|
| 54 | INTEGER ngrid,nlay |
|---|
| 55 | integer nq ! nombre de traceurs |
|---|
| 56 | REAL ptimestep ! pas de temps physique (s) |
|---|
| 57 | REAL pplay(ngrid,nlay) ! pression au milieu des couches (Pa) |
|---|
| 58 | |
|---|
| 59 | REAL pt(ngrid,nlay) ! temperature at the middle of the |
|---|
| 60 | ! layers (K) |
|---|
| 61 | REAL pdt(ngrid,nlay) ! tendance temperature des autres |
|---|
| 62 | ! param. |
|---|
| 63 | REAL pq(ngrid,nlay,nq) ! traceur (kg/kg) |
|---|
| 64 | REAL pdq(ngrid,nlay,nq) ! tendance avant condensation |
|---|
| 65 | ! (kg/kg.s-1) |
|---|
| 66 | REAL tauscaling(ngrid) ! Convertion factor for qdust and Ndust |
|---|
| 67 | |
|---|
| 68 | real rice(ngrid,nlay) ! Water Ice mass mean radius (m) |
|---|
| 69 | ! used for nucleation of CO2 on ice-coated ccns |
|---|
| 70 | |
|---|
| 71 | c Outputs: |
|---|
| 72 | REAL pdqcloudco2(ngrid,nlay,nq) ! tendance de la condensation |
|---|
| 73 | ! CO2 (kg/kg.s-1) |
|---|
| 74 | ! condensation si igcm_co2_ice |
|---|
| 75 | REAL pdtcloudco2(ngrid,nlay) ! tendance temperature due |
|---|
| 76 | ! a la chaleur latente |
|---|
| 77 | |
|---|
| 78 | c------------------------------------------------------------------ |
|---|
| 79 | c Local variables: |
|---|
| 80 | |
|---|
| 81 | LOGICAL firstcall |
|---|
| 82 | DATA firstcall/.true./ |
|---|
| 83 | SAVE firstcall |
|---|
| 84 | |
|---|
| 85 | REAL*8 derf ! Error function |
|---|
| 86 | !external derf |
|---|
| 87 | |
|---|
| 88 | !REAL*8 massflowrateCO2 |
|---|
| 89 | !external massflowrateCO2 |
|---|
| 90 | |
|---|
| 91 | INTEGER ig,l,i |
|---|
| 92 | |
|---|
| 93 | REAL zq(ngrid,nlay,nq) ! local value of tracers |
|---|
| 94 | REAL zq0(ngrid,nlay,nq) ! local initial value of tracers |
|---|
| 95 | REAL zt(ngrid,nlay) ! local value of temperature |
|---|
| 96 | REAL zqsat(ngrid,nlay) ! saturation vapor pressure for CO2 |
|---|
| 97 | REAL lw !Latent heat of sublimation (J.kg-1) |
|---|
| 98 | REAL l0,l1,l2,l3,l4 |
|---|
| 99 | REAL cste |
|---|
| 100 | REAL dMice ! mass of condensed ice |
|---|
| 101 | DOUBLE PRECISION sumcheck |
|---|
| 102 | DOUBLE PRECISION pco2 ! Co2 vapor partial pressure (Pa) |
|---|
| 103 | DOUBLE PRECISION satu ! Co2 vapor saturation ratio over ice |
|---|
| 104 | DOUBLE PRECISION Mo,No |
|---|
| 105 | DOUBLE PRECISION Rn, Rm, dev2, n_derf, m_derf |
|---|
| 106 | |
|---|
| 107 | ! Radius used by the microphysical scheme (m) |
|---|
| 108 | DOUBLE PRECISION n_aer(nbinco2_cld) ! number concentration volume-1 of particle/each size bin |
|---|
| 109 | DOUBLE PRECISION m_aer(nbinco2_cld) ! mass mixing ratio of particle/each size bin |
|---|
| 110 | |
|---|
| 111 | DOUBLE PRECISION n_aer_h2oice(nbinco2_cld) ! Same - for CO2 nucleation |
|---|
| 112 | DOUBLE PRECISION m_aer_h2oice(nbinco2_cld) ! Same - for CO2 nucleation |
|---|
| 113 | DOUBLE PRECISION rad_h2oice(nbinco2_cld) |
|---|
| 114 | |
|---|
| 115 | c REAL*8 sigco2 ! Co2-ice/air surface tension (N.m) |
|---|
| 116 | c EXTERNAL sigco2 |
|---|
| 117 | |
|---|
| 118 | DOUBLE PRECISION dN,dM, dNh2o, dMh2o, dNN,dMM |
|---|
| 119 | DOUBLE PRECISION rate(nbinco2_cld) ! nucleation rate |
|---|
| 120 | DOUBLE PRECISION rateh2o(nbinco2_cld) ! nucleation rate |
|---|
| 121 | REAL seq |
|---|
| 122 | |
|---|
| 123 | DOUBLE PRECISION riceco2(ngrid,nlay) ! CO2Ice mean radius (m) |
|---|
| 124 | |
|---|
| 125 | REAL rhocloudco2(ngrid,nlay) ! Cloud density (kg.m-3) |
|---|
| 126 | REAL rdust(ngrid,nlay) ! Dust geometric mean radius (m) |
|---|
| 127 | |
|---|
| 128 | c REAL res ! Resistance growth |
|---|
| 129 | DOUBLE PRECISION Ic_rice ! Mass transfer rate CO2 ice crystal |
|---|
| 130 | |
|---|
| 131 | |
|---|
| 132 | c Parameters of the size discretization |
|---|
| 133 | c used by the microphysical scheme |
|---|
| 134 | DOUBLE PRECISION, PARAMETER :: rmin_cld = 1.e-9 ! Minimum radius (m) |
|---|
| 135 | DOUBLE PRECISION, PARAMETER :: rmax_cld = 5.e-5 ! Maximum radius (m) |
|---|
| 136 | DOUBLE PRECISION, PARAMETER :: rbmin_cld =0.099e-9 |
|---|
| 137 | ! Minimum boundary radius (m) |
|---|
| 138 | DOUBLE PRECISION, PARAMETER :: rbmax_cld = 1.e-2 ! Maximum boundary radius (m) |
|---|
| 139 | DOUBLE PRECISION vrat_cld ! Volume ratio |
|---|
| 140 | DOUBLE PRECISION rb_cldco2(nbinco2_cld+1) ! boundary values of each rad_cldco2 bin (m) |
|---|
| 141 | SAVE rb_cldco2 |
|---|
| 142 | |
|---|
| 143 | DOUBLE PRECISION dr_cld(nbinco2_cld) ! width of each rad_cldco2 bin (m) |
|---|
| 144 | DOUBLE PRECISION vol_cld(nbinco2_cld) ! particle volume for each bin (m3) |
|---|
| 145 | |
|---|
| 146 | DOUBLE PRECISION Proba,Masse_atm,drsurdt,reff |
|---|
| 147 | REAL sigma_iceco2 ! Variance of the ice and CCN distributions |
|---|
| 148 | SAVE sigma_iceco2 |
|---|
| 149 | |
|---|
| 150 | DOUBLE PRECISION Niceco2,Qccnco2,Nccnco2 |
|---|
| 151 | |
|---|
| 152 | c---------------------------------- |
|---|
| 153 | c TESTS |
|---|
| 154 | |
|---|
| 155 | INTEGER countcells |
|---|
| 156 | |
|---|
| 157 | LOGICAL test_flag ! flag for test/debuging outputs |
|---|
| 158 | SAVE test_flag |
|---|
| 159 | |
|---|
| 160 | |
|---|
| 161 | REAL satubf(ngrid,nlay),satuaf(ngrid,nlay) |
|---|
| 162 | REAL res_out(ngrid,nlay) |
|---|
| 163 | |
|---|
| 164 | |
|---|
| 165 | c------------------------------------------------------------------ |
|---|
| 166 | |
|---|
| 167 | IF (firstcall) THEN |
|---|
| 168 | !============================================================= |
|---|
| 169 | ! 0. Definition of the size grid |
|---|
| 170 | !============================================================= |
|---|
| 171 | c rad_cldco2 is the primary radius grid used for microphysics computation. |
|---|
| 172 | c The grid spacing is computed assuming a constant volume ratio |
|---|
| 173 | c between two consecutive bins; i.e. vrat_cld. |
|---|
| 174 | c vrat_cld is determined from the boundary values of the size grid: |
|---|
| 175 | c rmin_cld and rmax_cld. |
|---|
| 176 | c The rb_cldco2 array contains the boundary values of each rad_cldco2 bin. |
|---|
| 177 | c dr_cld is the width of each rad_cldco2 bin. |
|---|
| 178 | |
|---|
| 179 | c Volume ratio between two adjacent bins |
|---|
| 180 | ! vrat_cld = log(rmax_cld/rmin_cld) / float(nbinco2_cld-1) *3. |
|---|
| 181 | ! vrat_cld = exp(vrat_cld) |
|---|
| 182 | vrat_cld = log(rmax_cld/rmin_cld) / float(nbinco2_cld-1) *3. |
|---|
| 183 | vrat_cld = exp(vrat_cld) |
|---|
| 184 | c write(*,*) "vrat_cld", vrat_cld |
|---|
| 185 | |
|---|
| 186 | rb_cldco2(1) = rbmin_cld |
|---|
| 187 | rad_cldco2(1) = rmin_cld |
|---|
| 188 | vol_cld(1) = 4./3. * dble(pi) * rmin_cld*rmin_cld*rmin_cld |
|---|
| 189 | ! vol_cld(1) = 4./3. * pi * rmin_cld*rmin_cld*rmin_cld |
|---|
| 190 | |
|---|
| 191 | do i=1,nbinco2_cld-1 |
|---|
| 192 | rad_cldco2(i+1) = rad_cldco2(i) * vrat_cld**(1./3.) |
|---|
| 193 | vol_cld(i+1) = vol_cld(i) * vrat_cld |
|---|
| 194 | enddo |
|---|
| 195 | |
|---|
| 196 | do i=1,nbinco2_cld |
|---|
| 197 | rb_cldco2(i+1)= ( (2.*vrat_cld) / (vrat_cld+1.) )**(1./3.) * |
|---|
| 198 | & rad_cldco2(i) |
|---|
| 199 | dr_cld(i) = rb_cldco2(i+1) - rb_cldco2(i) |
|---|
| 200 | enddo |
|---|
| 201 | rb_cldco2(nbinco2_cld+1) = rbmax_cld |
|---|
| 202 | dr_cld(nbinco2_cld) = rb_cldco2(nbinco2_cld+1) - |
|---|
| 203 | & rb_cldco2(nbinco2_cld) |
|---|
| 204 | |
|---|
| 205 | print*, ' ' |
|---|
| 206 | print*,'Microphysics co2: size bin information:' |
|---|
| 207 | print*,'i,rb_cldco2(i), rad_cldco2(i),dr_cld(i)' |
|---|
| 208 | print*,'-----------------------------------' |
|---|
| 209 | do i=1,nbinco2_cld |
|---|
| 210 | write(*,'(i3,3x,3(e12.6,4x))') i,rb_cldco2(i), rad_cldco2(i), |
|---|
| 211 | & dr_cld(i) |
|---|
| 212 | enddo |
|---|
| 213 | write(*,'(i3,3x,e12.6)') nbinco2_cld+1,rb_cldco2(nbinco2_cld+1) |
|---|
| 214 | print*,'-----------------------------------' |
|---|
| 215 | |
|---|
| 216 | do i=1,nbinco2_cld+1 |
|---|
| 217 | rb_cldco2(i) = log(rb_cldco2(i)) !! we save that so that it is not computed |
|---|
| 218 | !! at each timestep and gridpoint |
|---|
| 219 | enddo |
|---|
| 220 | |
|---|
| 221 | c Contact parameter of co2 ice on dst ( m=cos(theta) ) |
|---|
| 222 | c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
|---|
| 223 | c mteta = 0.952 |
|---|
| 224 | c mtetaco2 = 0.952 |
|---|
| 225 | c write(*,*) 'co2_param contact parameter:', mtetaco2 |
|---|
| 226 | |
|---|
| 227 | c Volume of a co2 molecule (m3) |
|---|
| 228 | vo1 = m0co2 / dble(rho_ice_co2) ! m0co2 et non mco2 |
|---|
| 229 | vo1co2=vo1 ! AJOUT JA |
|---|
| 230 | c Variance of the ice and CCN distributions |
|---|
| 231 | sigma_iceco2 = sqrt(log(1.+nuiceco2_sed)) |
|---|
| 232 | |
|---|
| 233 | |
|---|
| 234 | c write(*,*) 'Variance of ice & CCN distribs :', sigma_iceco2 |
|---|
| 235 | c write(*,*) 'nuice for sedimentation:', nuiceco2_sed |
|---|
| 236 | c write(*,*) 'Volume of a co2 molecule:', vo1co2 |
|---|
| 237 | |
|---|
| 238 | |
|---|
| 239 | test_flag = .false. |
|---|
| 240 | firstcall=.false. |
|---|
| 241 | |
|---|
| 242 | END IF |
|---|
| 243 | |
|---|
| 244 | |
|---|
| 245 | !============================================================= |
|---|
| 246 | ! 1. Initialisation |
|---|
| 247 | !============================================================= |
|---|
| 248 | !cste = 4*pi*rho_ice*ptimestep !not used for co2 |
|---|
| 249 | |
|---|
| 250 | res_out(:,:) = 0 |
|---|
| 251 | rice(:,:) = 1.e-11 |
|---|
| 252 | riceco2(:,:) = 1.e-11 |
|---|
| 253 | |
|---|
| 254 | c Initialize the tendencies |
|---|
| 255 | pdqcloudco2(1:ngrid,1:nlay,1:nq)=0. |
|---|
| 256 | pdtcloudco2(1:ngrid,1:nlay)=0. |
|---|
| 257 | |
|---|
| 258 | c pt temperature layer; pdt dT.s-1 |
|---|
| 259 | c pq traceur kg/kg; pdq tendance idem .s-1 |
|---|
| 260 | zt(1:ngrid,1:nlay) = |
|---|
| 261 | & pt(1:ngrid,1:nlay) + |
|---|
| 262 | & pdt(1:ngrid,1:nlay) * ptimestep |
|---|
| 263 | |
|---|
| 264 | zq(1:ngrid,1:nlay,1:nq) = |
|---|
| 265 | & pq(1:ngrid,1:nlay,1:nq) + |
|---|
| 266 | & pdq(1:ngrid,1:nlay,1:nq) * ptimestep |
|---|
| 267 | |
|---|
| 268 | |
|---|
| 269 | WHERE( zq(1:ngrid,1:nlay,1:nq) < 1.e-30 ) |
|---|
| 270 | & zq(1:ngrid,1:nlay,1:nq) = 1.e-30 |
|---|
| 271 | |
|---|
| 272 | zq0(1:ngrid,1:nlay,1:nq) = zq(1:ngrid,1:nlay,1:nq) |
|---|
| 273 | |
|---|
| 274 | !============================================================= |
|---|
| 275 | ! 2. Compute saturation |
|---|
| 276 | !============================================================= |
|---|
| 277 | |
|---|
| 278 | |
|---|
| 279 | dev2 = 1. / ( sqrt(2.) * sigma_iceco2 ) |
|---|
| 280 | |
|---|
| 281 | call co2sat(ngrid*nlay,zt,pplay,zqsat) !zqsat is psat(co2) |
|---|
| 282 | |
|---|
| 283 | countcells = 0 |
|---|
| 284 | |
|---|
| 285 | c Faire rice co2 update en n-1 puis a chaque microdt, mettre a jour riceco2 |
|---|
| 286 | |
|---|
| 287 | c Main loop over the GCM's grid |
|---|
| 288 | DO l=1,nlay |
|---|
| 289 | DO ig=1,ngrid |
|---|
| 290 | |
|---|
| 291 | |
|---|
| 292 | c Get the partial pressure of co2 vapor and its saturation ratio |
|---|
| 293 | pco2 = zq(ig,l,igcm_co2) * (mmean(ig,l)/44.01) * pplay(ig,l) |
|---|
| 294 | c satu = zq(ig,l,igcm_co2) / zqsat(ig,l) |
|---|
| 295 | satu = pco2 / zqsat(ig,l) |
|---|
| 296 | !============================================================= |
|---|
| 297 | ! 3. Nucleation |
|---|
| 298 | !============================================================= |
|---|
| 299 | |
|---|
| 300 | c call updaterccn(zq(ig,l,igcm_dust_mass), |
|---|
| 301 | c & zq(ig,l,igcm_dust_number),rdust(ig,l),tauscaling(ig)) |
|---|
| 302 | |
|---|
| 303 | IF ( satu .ge. 1d0 ) THEN ! if there is condensation |
|---|
| 304 | write(*,*) |
|---|
| 305 | write(*,*) "l, pco2, satu= ",l,pco2,satu |
|---|
| 306 | c Masse_atm=mmean(ig,l)*1.e-3*pplay(ig,l)/rgp/zt(ig,l) !Kg par couche |
|---|
| 307 | |
|---|
| 308 | |
|---|
| 309 | call updaterccn(zq(ig,l,igcm_dust_mass), |
|---|
| 310 | & zq(ig,l,igcm_dust_number),rdust(ig,l),tauscaling(ig)) |
|---|
| 311 | write(*,*) "Improved, l,Rdust = ",l,rdust(ig,l) |
|---|
| 312 | |
|---|
| 313 | rdust(ig,l)= zq(ig,l,igcm_dust_mass) |
|---|
| 314 | & *0.75/pi/rho_dust |
|---|
| 315 | & / zq(ig,l,igcm_dust_number) |
|---|
| 316 | rdust(ig,l)= rdust(ig,l)**(1./3.) |
|---|
| 317 | write(*,*) "Improved2, l,Rdust = ",l,rdust(ig,l) |
|---|
| 318 | rdust(ig,l)=max(1.e-9,rdust(ig,l)) |
|---|
| 319 | rdust(ig,l)=min(2e-6,rdust(ig,l)) |
|---|
| 320 | write(*,*) "Improved3, l,Rdust = ",l,rdust(ig,l) |
|---|
| 321 | |
|---|
| 322 | c Expand the dust moments into a binned distribution |
|---|
| 323 | Mo = zq(ig,l,igcm_dust_mass)* tauscaling(ig) |
|---|
| 324 | No = zq(ig,l,igcm_dust_number)* tauscaling(ig) |
|---|
| 325 | write(*,*) "dust number, mass = ", |
|---|
| 326 | & zq(ig,l,igcm_dust_number)* tauscaling(ig), |
|---|
| 327 | & zq(ig,l,igcm_dust_mass)* tauscaling(ig) |
|---|
| 328 | c write(*,*) "No, Mo = ",No, Mo |
|---|
| 329 | Rn = rdust(ig,l) |
|---|
| 330 | Rn = -log(Rn) |
|---|
| 331 | Rm = Rn - 3. * sigma_iceco2*sigma_iceco2 |
|---|
| 332 | n_derf = erf( (rb_cldco2(1)+Rn) *dev2) |
|---|
| 333 | m_derf = erf( (rb_cldco2(1)+Rm) *dev2) |
|---|
| 334 | |
|---|
| 335 | do i = 1, nbinco2_cld |
|---|
| 336 | n_aer(i) = -0.5 * No * n_derf !! this ith previously computed |
|---|
| 337 | m_aer(i) = -0.5 * Mo * m_derf !! this ith previously computed |
|---|
| 338 | n_derf = derf((rb_cldco2(i+1)+Rn) *dev2) |
|---|
| 339 | m_derf = derf((rb_cldco2(i+1)+Rm) *dev2) |
|---|
| 340 | n_aer(i) = n_aer(i) + 0.5 * No * n_derf |
|---|
| 341 | m_aer(i) = m_aer(i) + 0.5 * Mo * m_derf |
|---|
| 342 | c write(*,*) "i, rb_cldco2(i) = ",i, rb_cldco2(i),n_aer(i) |
|---|
| 343 | |
|---|
| 344 | enddo |
|---|
| 345 | |
|---|
| 346 | |
|---|
| 347 | sumcheck = 0 |
|---|
| 348 | do i = 1, nbinco2_cld |
|---|
| 349 | sumcheck = sumcheck + n_aer(i) |
|---|
| 350 | enddo |
|---|
| 351 | sumcheck = abs(sumcheck/No - 1) |
|---|
| 352 | if ((sumcheck .gt. 1e-5).and. (1./Rn .gt. rmin_cld)) then |
|---|
| 353 | print*, "WARNING, No sumcheck PROBLEM" |
|---|
| 354 | print*, "sumcheck, No",sumcheck, No |
|---|
| 355 | print*, "min radius, Rn, ig, l", rmin_cld, 1./Rn, ig, l |
|---|
| 356 | print*, "Dust binned distribution", n_aer |
|---|
| 357 | STOP |
|---|
| 358 | endif |
|---|
| 359 | |
|---|
| 360 | sumcheck = 0 |
|---|
| 361 | do i = 1, nbinco2_cld |
|---|
| 362 | sumcheck = sumcheck + m_aer(i) |
|---|
| 363 | enddo |
|---|
| 364 | sumcheck = abs(sumcheck/Mo - 1) |
|---|
| 365 | if ((sumcheck .gt. 1e-5) .and. (1./Rn .gt. rmin_cld)) |
|---|
| 366 | & then |
|---|
| 367 | print*, "WARNING, Mo sumcheck PROBLEM" |
|---|
| 368 | print*, "sumcheck, Mo",sumcheck, Mo |
|---|
| 369 | print*, "min radius, Rm, ig, l", rmin_cld, 1./Rm, ig,l |
|---|
| 370 | print*, "Dust binned distribution", m_aer |
|---|
| 371 | STOP |
|---|
| 372 | endif |
|---|
| 373 | |
|---|
| 374 | c Expand the water ice moments into a binned distribution |
|---|
| 375 | c For now the radius grid's bound are same as for dust |
|---|
| 376 | c min=100 nm and max=10microns |
|---|
| 377 | c might need a change if rice (water) is large (but how large?) - listo |
|---|
| 378 | |
|---|
| 379 | Mo = zq(ig,l,igcm_h2o_ice)* tauscaling(ig) + 1.e-30 |
|---|
| 380 | No = zq(ig,l,igcm_ccn_number)* tauscaling(ig) + 1.e-30 |
|---|
| 381 | Rn = rice(ig,l) |
|---|
| 382 | Rn = -log(Rn) |
|---|
| 383 | Rm = Rn - 3. * sigma_iceco2*sigma_iceco2 |
|---|
| 384 | n_derf = erf( (rb_cldco2(1)+Rn) *dev2) |
|---|
| 385 | m_derf = erf( (rb_cldco2(1)+Rm) *dev2) |
|---|
| 386 | do i = 1, nbinco2_cld |
|---|
| 387 | n_aer_h2oice(i) = -0.5 * No * n_derf !! this ith previously computed |
|---|
| 388 | m_aer_h2oice(i) = -0.5 * Mo * m_derf !! this ith previously computed |
|---|
| 389 | n_derf = derf( (rb_cldco2(i+1)+Rn) *dev2) |
|---|
| 390 | m_derf = derf( (rb_cldco2(i+1)+Rm) *dev2) |
|---|
| 391 | n_aer_h2oice(i) = n_aer(i) + 0.5 * No * n_derf ! vector not really needed - temp var - listo |
|---|
| 392 | m_aer_h2oice(i) = m_aer(i) + 0.5 * Mo * m_derf ! vector not really needed - temp var |
|---|
| 393 | rad_h2oice(i) = ((m_aer_h2oice(i)/rho_ice/ |
|---|
| 394 | & n_aer_h2oice(i) + vol_cld(i))) |
|---|
| 395 | & *0.75/pi**(1./3) |
|---|
| 396 | c write(*,*) "before nuc, i,rad_h2o(i)= ",i,rad_h2oice(i) |
|---|
| 397 | c & ,m_aer(i),n_aer(i) |
|---|
| 398 | enddo |
|---|
| 399 | |
|---|
| 400 | |
|---|
| 401 | |
|---|
| 402 | c Get the rates of nucleation |
|---|
| 403 | call nucleaCO2(dble(pco2),zt(ig,l),dble(satu) |
|---|
| 404 | & ,n_aer,rate,n_aer_h2oice |
|---|
| 405 | & ,rad_h2oice,rateh2o) |
|---|
| 406 | ! regarder rateh20, et mettre = 0 si non nul pour le moment |
|---|
| 407 | dN = 0. |
|---|
| 408 | dM = 0. |
|---|
| 409 | dNh2o = 0. |
|---|
| 410 | dMh2o = 0. |
|---|
| 411 | do i = 1, nbinco2_cld |
|---|
| 412 | !n_aer(i) = n_aer(i)/( 1. + (rate(i)+rateh2o(i))*ptimestep) |
|---|
| 413 | !m_aer(i) = m_aer(i)/( 1. + (rate(i)+rateh2o(i))*ptimestep) |
|---|
| 414 | Proba=1.0-dexp(-1.*ptimestep*rate(i)) |
|---|
| 415 | |
|---|
| 416 | |
|---|
| 417 | c dNh2o = dNh2o + n_aer_h2oice(i) * rateh2o(i) |
|---|
| 418 | c & * ptimestep |
|---|
| 419 | c dMh2o = dMh2o + m_aer_h2oice(i) * rateh2o(i) |
|---|
| 420 | c & *ptimestep |
|---|
| 421 | |
|---|
| 422 | dN = dN + n_aer(i) * Proba |
|---|
| 423 | dM = dM + m_aer(i) * Proba |
|---|
| 424 | c write(*,*) "i, dNi, dN= ",i,n_aer(i)*Proba,dN |
|---|
| 425 | enddo |
|---|
| 426 | |
|---|
| 427 | c do i=1,nbinco2_cld |
|---|
| 428 | c write(*,*) "i n_aer m_aer = ",i,n_aer(i),m_aer(i) |
|---|
| 429 | c enddo |
|---|
| 430 | ! dM masse activée (kg) et dN nb particules par kg d'air |
|---|
| 431 | |
|---|
| 432 | c Update Dust particles |
|---|
| 433 | c For CO2 ice : no subtraction from dust (neither for water ice particles) |
|---|
| 434 | ! zq(ig,l,igcm_dust_mass) = |
|---|
| 435 | ! & zq(ig,l,igcm_dust_mass) - dM/ tauscaling(ig) !max(tauscaling(ig),1.e-10) |
|---|
| 436 | ! zq(ig,l,igcm_dust_number) = |
|---|
| 437 | ! & zq(ig,l,igcm_dust_number) - dN/ tauscaling(ig) !max(tauscaling(ig),1.e-10) |
|---|
| 438 | c write(*,*) " nuclea dM = ",dM/tauscaling(ig), |
|---|
| 439 | c & " nuclea dN = ", dN/tauscaling(ig) |
|---|
| 440 | |
|---|
| 441 | dNN= dN/tauscaling(ig) |
|---|
| 442 | dMM= dM/tauscaling(ig) |
|---|
| 443 | |
|---|
| 444 | dNN=min(dNN,abs(zq0(ig,l,igcm_dust_number))) |
|---|
| 445 | dMM=min(dMM,zq0(ig,l,igcm_dust_mass)) |
|---|
| 446 | |
|---|
| 447 | c Update CCNs for CO2 crystals |
|---|
| 448 | ! WARNING dM dMh2o, interaction nuages eau-co2 -- h20 set to 0 for now |
|---|
| 449 | zq(ig,l,igcm_ccnco2_mass) = |
|---|
| 450 | & zq(ig,l,igcm_ccnco2_mass) + dMM |
|---|
| 451 | |
|---|
| 452 | zq(ig,l,igcm_ccnco2_number) = |
|---|
| 453 | & zq(ig,l,igcm_ccnco2_number) + dNN |
|---|
| 454 | |
|---|
| 455 | zq(ig,l,igcm_dust_mass) = |
|---|
| 456 | & zq(ig,l,igcm_dust_mass) - dMM |
|---|
| 457 | zq(ig,l,igcm_dust_number) = |
|---|
| 458 | & zq(ig,l,igcm_dust_number) - dNN |
|---|
| 459 | |
|---|
| 460 | |
|---|
| 461 | c + enlever les CCN a la distri de dust |
|---|
| 462 | |
|---|
| 463 | write(*,*) "new dust_mass, number =", |
|---|
| 464 | & zq(ig,l,igcm_dust_mass)* tauscaling(ig), |
|---|
| 465 | & zq(ig,l,igcm_dust_number)*tauscaling(ig) |
|---|
| 466 | write(*,*) "new ccn mass, number =", |
|---|
| 467 | & zq(ig,l,igcm_ccnco2_mass)* tauscaling(ig) |
|---|
| 468 | & ,zq(ig,l,igcm_ccnco2_number)*tauscaling(ig) |
|---|
| 469 | |
|---|
| 470 | ENDIF ! of is satu >1 |
|---|
| 471 | !============================================================= |
|---|
| 472 | ! 4. Ice growth: scheme for radius evolution |
|---|
| 473 | !============================================================= |
|---|
| 474 | |
|---|
| 475 | c We trigger crystal growth if and only if there is at least one nuclei (N>1). |
|---|
| 476 | c Indeed, if we are supersaturated and still don't have at least one nuclei, we should better wait |
|---|
| 477 | c to avoid unrealistic value for nuclei radius and so on for cases that remain negligible. |
|---|
| 478 | c IF ( zq(ig,l,igcm_ccnco2_number)*tauscaling(ig).ge. 1.0) THEN |
|---|
| 479 | |
|---|
| 480 | IF (zq(ig,l,igcm_ccnco2_number)*tauscaling(ig) .ge. threshJA) |
|---|
| 481 | & THEN ! we trigger crystal growth |
|---|
| 482 | |
|---|
| 483 | call updaterccn(zq(ig,l,igcm_dust_mass), |
|---|
| 484 | & zq(ig,l,igcm_dust_number),rdust(ig,l),tauscaling(ig)) |
|---|
| 485 | rdust(ig,l)= zq(ig,l,igcm_ccnco2_mass) |
|---|
| 486 | & *0.75/pi/rho_dust |
|---|
| 487 | & / zq(ig,l,igcm_ccnco2_number) |
|---|
| 488 | rdust(ig,l)= rdust(ig,l)**(1./3.) |
|---|
| 489 | |
|---|
| 490 | rdust(ig,l)=max(1.e-10,rdust(ig,l)) |
|---|
| 491 | rdust(ig,l)=min(2e-6,rdust(ig,l)) |
|---|
| 492 | ! rdust(ig,l)=1.e-7 |
|---|
| 493 | |
|---|
| 494 | IF (zq(ig,l,igcm_ccnco2_mass) .lt. 0. .or. |
|---|
| 495 | & zq(ig,l,igcm_ccnco2_number) .lt. 0. .or. |
|---|
| 496 | & zq(ig,l,igcm_dust_mass) .lt. 0. .or. |
|---|
| 497 | & zq(ig,l,igcm_dust_number) .lt. 0. ) THEN |
|---|
| 498 | |
|---|
| 499 | write(*,*) "before growth CCN N,M = " |
|---|
| 500 | & ,zq(ig,l,igcm_ccnco2_number)*tauscaling(ig) |
|---|
| 501 | & ,zq(ig,l,igcm_ccnco2_mass)*tauscaling(ig) |
|---|
| 502 | |
|---|
| 503 | write(*,*) "before growth dust number mass = ", |
|---|
| 504 | & zq(ig,l,igcm_dust_number)*tauscaling(ig), |
|---|
| 505 | & zq(ig,l,igcm_dust_mass)*tauscaling(ig) |
|---|
| 506 | STOP |
|---|
| 507 | END IF |
|---|
| 508 | |
|---|
| 509 | c write(*,*) "reff dN = ",reff,dN |
|---|
| 510 | c reff=reff/dble(dN) |
|---|
| 511 | c if (zq(ig,l,igcm_co2_ice) .le. 1.e-20) then |
|---|
| 512 | c riceco2(ig,l)=reff |
|---|
| 513 | c endif |
|---|
| 514 | |
|---|
| 515 | c write(*,*) "Rdust in improved = ",rdust(ig,l) |
|---|
| 516 | |
|---|
| 517 | riceco2(ig,l)=( zq(ig,l,igcm_co2_ice)*3.0/ |
|---|
| 518 | & (4.0*rho_ice_co2*pi*zq(ig,l,igcm_ccnco2_number) |
|---|
| 519 | & *tauscaling(ig)) +rdust(ig,l)*rdust(ig,l) |
|---|
| 520 | & *rdust(ig,l) )**(1.0/3.0) |
|---|
| 521 | |
|---|
| 522 | ! WATCH OUT: CO2 nuclei is supposed to be dust |
|---|
| 523 | ! only when deriving rhocloud (otherwise would need to keep info on water embedded in co2) - listo |
|---|
| 524 | write(*,*) "Rdust before growth = ",rdust(ig,l) |
|---|
| 525 | write(*,*) "Riceco2 before growth = ",riceco2(ig,l) |
|---|
| 526 | |
|---|
| 527 | !! Niceco2,Qccnco2,Nccnco2 |
|---|
| 528 | Niceco2 = zq(ig,l,igcm_co2_ice) |
|---|
| 529 | Qccnco2 = zq(ig,l,igcm_ccnco2_mass) |
|---|
| 530 | Nccnco2 = zq(ig,l,igcm_ccnco2_number) |
|---|
| 531 | call updaterice_microCO2(Niceco2,Qccnco2,Nccnco2, |
|---|
| 532 | & tauscaling(ig),riceco2(ig,l),rhocloudco2(ig,l)) |
|---|
| 533 | write(*,*) "Riceco2 update before growth = ",riceco2(ig,l) |
|---|
| 534 | |
|---|
| 535 | No = zq(ig,l,igcm_ccnco2_number)* tauscaling(ig) |
|---|
| 536 | & + 1.e-30 |
|---|
| 537 | ! No nb de particules de poussieres mis à l'échelle pour donner une opacité optique |
|---|
| 538 | |
|---|
| 539 | c saturation at equilibrium |
|---|
| 540 | c rice should not be too small, otherwise seq value is not valid |
|---|
| 541 | c seq = exp(2.*sigco2*mco2 / (rho_ice_co2*rgp*zt(ig,l)* |
|---|
| 542 | c & max(riceco2(ig,l),1.e-7))) !Exponant sans unité OK |
|---|
| 543 | |
|---|
| 544 | ccccccc Scheme of microphys. mass growth for CO2 |
|---|
| 545 | |
|---|
| 546 | call massflowrateCO2(pplay(ig,l),zt(ig,l), |
|---|
| 547 | & satu,riceco2(ig,l),mmean(ig,l),Ic_rice) ! Mass transfer rate (kg/s) for a rice particle |
|---|
| 548 | ! Ic_rice mass flux kg.s-1 <0 si croissance ! |
|---|
| 549 | drsurdt=-1.0/(4.0*pi*riceco2(ig,l)* |
|---|
| 550 | & riceco2(ig,l)*rho_ice_co2)*Ic_rice |
|---|
| 551 | dMice = No * Ic_rice * ptimestep ! Kg par kg d'air, <0 si croissance ! |
|---|
| 552 | write(*,*) "dMicev0 in improved = " , dMice |
|---|
| 553 | |
|---|
| 554 | if (dMice .gt. 0) dMice = min(dMice,zq0(ig,l,igcm_co2_ice)) |
|---|
| 555 | if (dMice .lt. 0) dMice = max(dMice,-1.*zq0(ig,l,igcm_co2)) |
|---|
| 556 | |
|---|
| 557 | |
|---|
| 558 | |
|---|
| 559 | |
|---|
| 560 | riceco2(ig,l)=riceco2(ig,l)+drsurdt*ptimestep |
|---|
| 561 | c write(*,*) "riceco2+dr/dt = ", riceco2(ig,l) |
|---|
| 562 | |
|---|
| 563 | write(*,*) "dMice in improved = " , dMice |
|---|
| 564 | |
|---|
| 565 | |
|---|
| 566 | zq(ig,l,igcm_co2_ice) = zq(ig,l,igcm_co2_ice) |
|---|
| 567 | & -dMice |
|---|
| 568 | zq(ig,l,igcm_co2) = zq(ig,l,igcm_co2)+dMice |
|---|
| 569 | c write(*,*) "zq co2 ice = ", zq(ig,l,igcm_co2_ice) |
|---|
| 570 | countcells = countcells + 1 |
|---|
| 571 | |
|---|
| 572 | riceco2(ig,l)=( zq(ig,l,igcm_co2_ice)*3.0/ |
|---|
| 573 | & (4.0*rho_ice_co2*pi*zq(ig,l,igcm_ccnco2_number) |
|---|
| 574 | & *tauscaling(ig)) +rdust(ig,l)*rdust(ig,l) |
|---|
| 575 | & *rdust(ig,l) )**(1.0/3.0) |
|---|
| 576 | write(*,*) "new riceco2 = ",riceco2(ig,l) |
|---|
| 577 | |
|---|
| 578 | !! Niceco2,Qccnco2,Nccnco2 |
|---|
| 579 | Niceco2 = zq(ig,l,igcm_co2_ice) |
|---|
| 580 | Qccnco2 = zq(ig,l,igcm_ccnco2_mass) |
|---|
| 581 | Nccnco2 = zq(ig,l,igcm_ccnco2_number) |
|---|
| 582 | call updaterice_microCO2(Niceco2,Qccnco2,Nccnco2, |
|---|
| 583 | & tauscaling(ig),riceco2(ig,l),rhocloudco2(ig,l)) |
|---|
| 584 | write(*,*) "new riceco2 updaterad= ",riceco2(ig,l) |
|---|
| 585 | |
|---|
| 586 | ! latent heat release |
|---|
| 587 | |
|---|
| 588 | l0=595594. |
|---|
| 589 | l1=903.111 |
|---|
| 590 | l2=-11.5959 |
|---|
| 591 | l3=0.0528288 |
|---|
| 592 | l4=-0.000103183 |
|---|
| 593 | |
|---|
| 594 | lw = l0 + l1 * zt(ig,l) + l2 * zt(ig,l)**2 + |
|---|
| 595 | & l3 * zt(ig,l)**3 + l4 * zt(ig,l)**4 !J.kg-1 |
|---|
| 596 | c write(*,*) "CPP= ",cpp ! = 744.5 |
|---|
| 597 | |
|---|
| 598 | pdtcloudco2(ig,l)= dMice*lw/cpp/ptimestep ! kg par couche * J par kg /J par K / s = K par seconde |
|---|
| 599 | |
|---|
| 600 | c write(*,*) "pdtcloudco2 = ",pdtcloudco2(ig,l) |
|---|
| 601 | |
|---|
| 602 | c Peut etre -1*dMice? |
|---|
| 603 | |
|---|
| 604 | |
|---|
| 605 | !deltaT par condens/subli. qui remplace le dT du CO2 du newcondens pré-constantino |
|---|
| 606 | !PDT should be in K/s. |
|---|
| 607 | !============================================================= |
|---|
| 608 | ! 5. Dust cores released, tendancies, latent heat, etc ... |
|---|
| 609 | !============================================================= |
|---|
| 610 | |
|---|
| 611 | c If all the ice particles sublimate, all the condensation |
|---|
| 612 | c nuclei are released: |
|---|
| 613 | |
|---|
| 614 | c !!! with CO2 ice nuclei: dust/H2O nuclei are not released because |
|---|
| 615 | c they were not subtracted to dust_number |
|---|
| 616 | c Their counter is just set to "0". |
|---|
| 617 | c (see end of section 3.) : On peut les enlever à dust |
|---|
| 618 | |
|---|
| 619 | c interaction ho-co2 ici, dans la mise a jour des traceurs WARNING reflechir |
|---|
| 620 | |
|---|
| 621 | |
|---|
| 622 | ENDIF !of if Nccn>1 |
|---|
| 623 | c if (riceco2(ig,l) .lt. 1.e-9) then |
|---|
| 624 | |
|---|
| 625 | if (zq(ig,l,igcm_co2_ice).le.1.e-20 .or. |
|---|
| 626 | & riceco2(ig,l) .lt. 1.e-9 .or. riceco2(ig,l) |
|---|
| 627 | & .ge. 4.999e-4) then |
|---|
| 628 | ! Reverser les ccn libérés dans les h2o ou dust? |
|---|
| 629 | |
|---|
| 630 | c c ICI: Co2 ice devient vapeur |
|---|
| 631 | zq(ig,l,igcm_co2) = zq(ig,l,igcm_co2) |
|---|
| 632 | & + zq(ig,l,igcm_co2_ice) |
|---|
| 633 | |
|---|
| 634 | zq(ig,l,igcm_dust_mass) = |
|---|
| 635 | & zq(ig,l,igcm_dust_mass) |
|---|
| 636 | & + zq(ig,l,igcm_ccnco2_mass) |
|---|
| 637 | zq(ig,l,igcm_dust_number) |
|---|
| 638 | & = zq(ig,l,igcm_dust_number) |
|---|
| 639 | & + zq(ig,l,igcm_ccnco2_number) |
|---|
| 640 | c c CCNs |
|---|
| 641 | zq(ig,l,igcm_ccnco2_mass) = 0. |
|---|
| 642 | zq(ig,l,igcm_ccnco2_number) =0. |
|---|
| 643 | zq(ig,l,igcm_co2_ice) = 0. |
|---|
| 644 | riceco2(ig,l)=0. |
|---|
| 645 | endif |
|---|
| 646 | |
|---|
| 647 | c write(*,*) "zq co2 end imp= ", zq(i g,l,igcm_co2_ice),satu |
|---|
| 648 | |
|---|
| 649 | |
|---|
| 650 | |
|---|
| 651 | ENDDO ! of ig loop |
|---|
| 652 | ENDDO ! of nlayer loop |
|---|
| 653 | |
|---|
| 654 | |
|---|
| 655 | ! Get cloud tendencies |
|---|
| 656 | pdqcloudco2(1:ngrid,1:nlay,igcm_co2) = |
|---|
| 657 | & (zq(1:ngrid,1:nlay,igcm_co2) - |
|---|
| 658 | & zq0(1:ngrid,1:nlay,igcm_co2))/ptimestep |
|---|
| 659 | |
|---|
| 660 | pdqcloudco2(1:ngrid,1:nlay,igcm_co2_ice) = |
|---|
| 661 | & (zq(1:ngrid,1:nlay,igcm_co2_ice) - |
|---|
| 662 | & zq0(1:ngrid,1:nlay,igcm_co2_ice))/ptimestep |
|---|
| 663 | c do l=1,nlay |
|---|
| 664 | c write(*,*) "end imp",pdqcloudco2(1,l,igcm_co2), |
|---|
| 665 | c & pdqcloudco2(1,l,igcm_co2_ice) |
|---|
| 666 | c enddo |
|---|
| 667 | pdqcloudco2(1:ngrid,1:nlay,igcm_ccnco2_mass) = |
|---|
| 668 | & (zq(1:ngrid,1:nlay,igcm_ccnco2_mass) - |
|---|
| 669 | & zq0(1:ngrid,1:nlay,igcm_ccnco2_mass))/ptimestep |
|---|
| 670 | |
|---|
| 671 | pdqcloudco2(1:ngrid,1:nlay,igcm_ccnco2_number) = |
|---|
| 672 | & (zq(1:ngrid,1:nlay,igcm_ccnco2_number) - |
|---|
| 673 | & zq0(1:ngrid,1:nlay,igcm_ccnco2_number))/ptimestep |
|---|
| 674 | |
|---|
| 675 | |
|---|
| 676 | if (scavenging) then |
|---|
| 677 | |
|---|
| 678 | pdqcloudco2(1:ngrid,1:nlay,igcm_dust_mass) = |
|---|
| 679 | & (zq(1:ngrid,1:nlay,igcm_dust_mass) - |
|---|
| 680 | & zq0(1:ngrid,1:nlay,igcm_dust_mass))/ptimestep |
|---|
| 681 | pdqcloudco2(1:ngrid,1:nlay,igcm_dust_number) = |
|---|
| 682 | & (zq(1:ngrid,1:nlay,igcm_dust_number) - |
|---|
| 683 | & zq0(1:ngrid,1:nlay,igcm_dust_number))/ptimestep |
|---|
| 684 | endif |
|---|
| 685 | |
|---|
| 686 | c call WRITEdiagfi(ngrid,"Improvedb","co2 traceur","kg/kg",1, |
|---|
| 687 | c & zq0(1,:,igcm_co2_ice) ) |
|---|
| 688 | |
|---|
| 689 | c call WRITEdiagfi(ngrid,"Improveda","co2 traceur","kg/kg",1, |
|---|
| 690 | c & zq(1,:,igcm_co2_ice) ) |
|---|
| 691 | |
|---|
| 692 | |
|---|
| 693 | |
|---|
| 694 | |
|---|
| 695 | c call WRITEdiagfi(ngrid,"satuco2","satu co2 in improved","kg/kg",1, |
|---|
| 696 | c & satu) |
|---|
| 697 | |
|---|
| 698 | |
|---|
| 699 | !!!!!!!!!!!!!! TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS |
|---|
| 700 | !!!!!!!!!!!!!! TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS |
|---|
| 701 | !!!!!!!!!!!!!! TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS |
|---|
| 702 | IF (test_flag) then |
|---|
| 703 | |
|---|
| 704 | ! error2d(:) = 0. |
|---|
| 705 | DO l=1,nlay |
|---|
| 706 | DO ig=1,ngrid |
|---|
| 707 | ! error2d(ig) = max(abs(error_out(ig,l)),error2d(ig)) |
|---|
| 708 | satubf(ig,l) = zq0(ig,l,igcm_co2)/zqsat(ig,l) ! att. if zqsat=mmr or psat |
|---|
| 709 | satuaf(ig,l) = zq(ig,l,igcm_co2)/zqsat(ig,l) |
|---|
| 710 | ENDDO |
|---|
| 711 | ENDDO |
|---|
| 712 | |
|---|
| 713 | write(*,*) 'count is ',countcells, ' i.e. ', |
|---|
| 714 | & countcells*100/(nlay*ngrid), '% for microphys computation' |
|---|
| 715 | |
|---|
| 716 | ! IF (ngrid.ne.1) THEN ! 3D |
|---|
| 717 | ! call WRITEDIAGFI(ngrid,"satu","ratio saturation","",3, |
|---|
| 718 | ! & satu_out) |
|---|
| 719 | ! call WRITEDIAGFI(ngrid,"dM","ccn variation","kg/kg",3, |
|---|
| 720 | ! & dM_out) |
|---|
| 721 | ! call WRITEDIAGFI(ngrid,"dN","ccn variation","#",3, |
|---|
| 722 | ! & dN_out) |
|---|
| 723 | ! call WRITEDIAGFI(ngrid,"error","dichotomy max error","%",2, |
|---|
| 724 | ! & error2d) |
|---|
| 725 | ! call WRITEDIAGFI(ngrid,"zqsat","zqsat","kg",3, |
|---|
| 726 | ! & zqsat) |
|---|
| 727 | ! ENDIF |
|---|
| 728 | |
|---|
| 729 | ! IF (ngrid.eq.1) THEN ! 1D |
|---|
| 730 | ! call WRITEDIAGFI(ngrid,"error","incertitude sur glace","%",1, |
|---|
| 731 | ! & error_out) |
|---|
| 732 | ! call WRITEdiagfi(ngrid,"resist","resistance","s/m2",1, |
|---|
| 733 | ! & res_out) |
|---|
| 734 | call WRITEdiagfi(ngrid,"satu_bf","satu before","kg/kg",1, |
|---|
| 735 | & satubf) |
|---|
| 736 | call WRITEdiagfi(ngrid,"satu_af","satu after","kg/kg",1, |
|---|
| 737 | & satuaf) |
|---|
| 738 | call WRITEdiagfi(ngrid,"vapbf","CO2vap before","kg/kg",1, |
|---|
| 739 | & zq0(1,1,igcm_co2)) |
|---|
| 740 | call WRITEdiagfi(ngrid,"vapaf","CO2vap after","kg/kg",1, |
|---|
| 741 | & zq(1,1,igcm_co2)) |
|---|
| 742 | call WRITEdiagfi(ngrid,"icebf","CO2ice before","kg/kg",1, |
|---|
| 743 | & zq0(1,1,igcm_co2_ice)) |
|---|
| 744 | call WRITEdiagfi(ngrid,"iceaf","CO2ice after","kg/kg",1, |
|---|
| 745 | & zq(1,1,igcm_co2_ice)) |
|---|
| 746 | call WRITEdiagfi(ngrid,"ccnbf","ccn before","/kg",1, |
|---|
| 747 | & zq0(1,1,igcm_ccnco2_number)) |
|---|
| 748 | call WRITEdiagfi(ngrid,"ccnaf","ccn after","/kg",1, |
|---|
| 749 | & zq(1,1,igcm_ccnco2_number)) |
|---|
| 750 | c call WRITEDIAGFI(ngrid,"growthrate","growth rate","m^2/s",1, |
|---|
| 751 | c & gr_out) |
|---|
| 752 | c call WRITEDIAGFI(ngrid,"nuclearate","nucleation rate","",1, |
|---|
| 753 | c & rate_out) |
|---|
| 754 | c call WRITEDIAGFI(ngrid,"dM","ccn variation","kg",1, |
|---|
| 755 | c & dM_out) |
|---|
| 756 | c call WRITEDIAGFI(ngrid,"dN","ccn variation","#",1, |
|---|
| 757 | c & dN_out) |
|---|
| 758 | c call WRITEdiagfi(ngrid,"zqsat","p vap sat","kg/kg",1, |
|---|
| 759 | c & zqsat) |
|---|
| 760 | ! call WRITEDIAGFI(ngrid,"satu","ratio saturation","",1, |
|---|
| 761 | ! & satu_out) |
|---|
| 762 | ! call WRITEDIAGFI(ngrid,"rdust_sca","rdust","m",1, |
|---|
| 763 | ! & rdust) |
|---|
| 764 | ! call WRITEDIAGFI(ngrid,"rsedcloud","rsedcloud","m",1, |
|---|
| 765 | ! & rsedcloud) |
|---|
| 766 | ! call WRITEDIAGFI(ngrid,"rhocloud","rhocloud","kg.m-3",1, |
|---|
| 767 | ! & rhocloud) |
|---|
| 768 | ! ENDIF |
|---|
| 769 | |
|---|
| 770 | ENDIF ! endif test_flag |
|---|
| 771 | !!!!!!!!!!!!!! TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS |
|---|
| 772 | !!!!!!!!!!!!!! TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS |
|---|
| 773 | !!!!!!!!!!!!!! TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS |
|---|
| 774 | |
|---|
| 775 | return |
|---|
| 776 | end |
|---|
| 777 | |
|---|
| 778 | |
|---|
| 779 | |
|---|
| 780 | cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc |
|---|
| 781 | cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc |
|---|
| 782 | c The so -called "phi" function is such as phi(r) - phi(r0) = t - t0 |
|---|
| 783 | c It is an analytical solution to the ice radius growth equation, |
|---|
| 784 | c with the approximation of a constant 'reduced' cunningham correction factor |
|---|
| 785 | c (lambda in growthrate.F) taken at radius req instead of rice |
|---|
| 786 | cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc |
|---|
| 787 | cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc |
|---|
| 788 | |
|---|
| 789 | c subroutine phi(rice,req,coeff1,coeff2,time) |
|---|
| 790 | c |
|---|
| 791 | c implicit none |
|---|
| 792 | c |
|---|
| 793 | c ! inputs |
|---|
| 794 | c real rice ! ice radius |
|---|
| 795 | c real req ! ice radius at equilibirum |
|---|
| 796 | c real coeff1 ! coeff for the log |
|---|
| 797 | c real coeff2 ! coeff for the arctan |
|---|
| 798 | c |
|---|
| 799 | c ! output |
|---|
| 800 | c real time |
|---|
| 801 | c |
|---|
| 802 | c !local |
|---|
| 803 | c real var |
|---|
| 804 | c |
|---|
| 805 | c ! 1.73205 is sqrt(3) |
|---|
| 806 | c |
|---|
| 807 | c var = max( |
|---|
| 808 | c & abs(rice-req) / sqrt(rice*rice + rice*req + req*req),1e-30) |
|---|
| 809 | c |
|---|
| 810 | c time = |
|---|
| 811 | c & coeff1 * |
|---|
| 812 | c & log( var ) |
|---|
| 813 | c & + coeff2 * 1.73205 * |
|---|
| 814 | c & atan( (2*rice+req) / (1.73205*req) ) |
|---|
| 815 | c |
|---|
| 816 | c return |
|---|
| 817 | c end |
|---|
| 818 | |
|---|
| 819 | |
|---|
| 820 | |
|---|
