Changeset 3187
- Timestamp:
- Jan 26, 2024, 3:58:16 PM (10 months ago)
- File:
-
- 1 edited
-
trunk/LMDZ.PLUTO/libf/phypluto/condense_n2.F90 (modified) (10 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/LMDZ.PLUTO/libf/phypluto/condense_n2.F90
r3184 r3187 124 124 125 125 126 REAL,SAVE :: latcond= 5.9e5126 REAL,SAVE :: latcond=2.5e5 127 127 REAL,SAVE :: ccond 128 128 REAL,SAVE,ALLOCATABLE,DIMENSION(:) :: emisref … … 142 142 IF (firstcall) THEN 143 143 144 ALLOCATE(emisref(ngrid))145 ! Find N2 ice tracer.146 do iq=1,nq147 tracername=noms(iq)148 if (tracername.eq."n2_ice") then149 i_n2ice=iq150 endif151 enddo144 ! ALLOCATE(emisref(ngrid)) 145 ! ! Find N2 ice tracer. 146 ! do iq=1,nq 147 ! tracername=noms(iq) 148 ! if (tracername.eq."n2_ice") then 149 ! i_n2ice=iq 150 ! endif 151 ! enddo 152 152 153 write(*,*) "condense_n2: i_n2ice=",i_n2ice154 155 if((i_n2ice.lt.1))then156 print*,'In condens_cloud but no N2 ice tracer, exiting.'157 print*,'Still need generalisation to arbitrary species!'158 stop159 endif153 ! write(*,*) "condense_n2: i_n2ice=",i_n2ice 154 155 ! if((i_n2ice.lt.1))then 156 ! print*,'In condens_cloud but no N2 ice tracer, exiting.' 157 ! print*,'Still need generalisation to arbitrary species!' 158 ! stop 159 ! endif 160 160 161 161 ccond=cpp/(g*latcond) … … 223 223 224 224 ! Forecast the atmospheric frost temperature 'ztcond' and nucleation temperature 'ztnuc'. 225 DO l=1,nlayer226 DO ig=1,ngrid227 ppn2=gfrac(igas_N2)*pplay(ig,l)228 call get_tcond_n2(ppn2,ztcond(ig,l))229 call get_tnuc_n2(ppn2,ztnuc(ig,l))230 ENDDO231 ENDDO225 ! DO l=1,nlayer 226 ! DO ig=1,ngrid 227 ! ppn2=gfrac(igas_N2)*pplay(ig,l) 228 ! call get_tcond_n2(ppn2,ztcond(ig,l)) 229 ! call get_tnuc_n2(ppn2,ztnuc(ig,l)) 230 ! ENDDO 231 ! ENDDO 232 232 233 233 ! Initialize zq and zt at the beginning of the sub-timestep loop and qsurf. 234 DO ig=1,ngrid235 picen2(ig)=pqsurf(ig,i_n2ice)236 DO l=1,nlayer237 zt(ig,l)=pt(ig,l)238 zq(ig,l,i_n2ice)=pq(ig,l,i_n2ice)239 IF( zq(ig,l,i_n2ice).lt.-1.e-6 ) THEN240 print*,'Uh-oh, zq = ',zq(ig,l,i_n2ice),'at ig,l=',ig,l241 if(l.eq.1)then242 print*,'Perhaps the atmosphere is collapsing on surface...?'243 endif244 END IF245 ENDDO246 ENDDO234 ! DO ig=1,ngrid 235 ! picen2(ig)=pqsurf(ig,i_n2ice) 236 ! DO l=1,nlayer 237 ! zt(ig,l)=pt(ig,l) 238 ! zq(ig,l,i_n2ice)=pq(ig,l,i_n2ice) 239 ! IF( zq(ig,l,i_n2ice).lt.-1.e-6 ) THEN 240 ! print*,'Uh-oh, zq = ',zq(ig,l,i_n2ice),'at ig,l=',ig,l 241 ! if(l.eq.1)then 242 ! print*,'Perhaps the atmosphere is collapsing on surface...?' 243 ! endif 244 ! END IF 245 ! ENDDO 246 ! ENDDO 247 247 248 248 ! Calculate the mass of each atmospheric layer (kg.m-2) 249 do ilay=1,nlayer250 DO ig=1,ngrid251 masse(ig,ilay)=(pplev(ig,ilay) - pplev(ig,ilay+1)) /g252 end do253 end do254 249 ! do ilay=1,nlayer 250 ! DO ig=1,ngrid 251 ! masse(ig,ilay)=(pplev(ig,ilay) - pplev(ig,ilay+1)) /g 252 ! end do 253 ! end do 254 ! 255 255 256 256 !----------------------------------------------------------- … … 259 259 260 260 261 Ntime = 20 ! number of sub-timestep262 subptimestep = ptimestep/float(Ntime)261 ! Ntime = 20 ! number of sub-timestep 262 ! subptimestep = ptimestep/float(Ntime) 263 263 264 264 ! Add the tendencies from other physical processes at each subtimstep. 265 DO it=1,Ntime266 DO l=1,nlayer267 DO ig=1,ngrid268 zt(ig,l) = zt(ig,l) + pdt(ig,l) * subptimestep269 zq(ig,l,i_n2ice) = zq(ig,l,i_n2ice) + pdq(ig,l,i_n2ice) * subptimestep270 END DO271 END DO265 ! DO it=1,Ntime 266 ! DO l=1,nlayer 267 ! DO ig=1,ngrid 268 ! zt(ig,l) = zt(ig,l) + pdt(ig,l) * subptimestep 269 ! zq(ig,l,i_n2ice) = zq(ig,l,i_n2ice) + pdq(ig,l,i_n2ice) * subptimestep 270 ! END DO 271 ! END DO 272 272 273 273 ! Gravitational sedimentation starts. 274 274 275 275 ! Sedimentation computed from radius computed from q in module radii_mod. 276 call n2_reffrad(ngrid,nlayer,nq,zq,reffrad)277 278 DO ilay=1,nlayer279 DO ig=1,ngrid280 281 reff = reffrad(ig,ilay)282 283 284 (pplev(ig,ilay),pt(ig,ilay), &285 reff,vstokes,rho_n2)286 287 !w(ig,ilay,i_n2ice) = 0.0288 w(ig,ilay,i_n2ice) = vstokes * subptimestep * &289 pplev(ig,ilay)/(r*pt(ig,ilay))290 291 END DO292 END DO276 ! call n2_reffrad(ngrid,nlayer,nq,zq,reffrad) 277 ! 278 ! DO ilay=1,nlayer 279 ! DO ig=1,ngrid 280 281 ! reff = reffrad(ig,ilay) 282 283 ! call stokes & 284 ! (pplev(ig,ilay),pt(ig,ilay), & 285 ! reff,vstokes,rho_n2) 286 287 ! !w(ig,ilay,i_n2ice) = 0.0 288 ! w(ig,ilay,i_n2ice) = vstokes * subptimestep * & 289 ! pplev(ig,ilay)/(r*pt(ig,ilay)) 290 291 ! END DO 292 ! END DO 293 293 294 294 ! Computing q after sedimentation 295 call vlz_fi(ngrid,nlayer,zq(1,1,i_n2ice),2.,masse,w(1,1,i_n2ice),wq)295 ! call vlz_fi(ngrid,nlayer,zq(1,1,i_n2ice),2.,masse,w(1,1,i_n2ice),wq) 296 296 297 297 298 298 ! Progressively accumulating the flux to the ground. 299 299 ! Mfallice is the total amount of ice fallen to the ground. 300 DO ig=1,ngrid301 Mfallice(ig) = Mfallice(ig) + wq(ig,i_n2ice)302 END DO300 ! DO ig=1,ngrid 301 ! Mfallice(ig) = Mfallice(ig) + wq(ig,i_n2ice) 302 ! END DO 303 303 304 304 !---------------------------------------------------------- … … 308 308 309 309 310 DO l=nlayer , 1, -1311 DO ig=1,ngrid312 pdtc(ig,l)=0.310 ! DO l=nlayer , 1, -1 311 ! DO ig=1,ngrid 312 ! pdtc(ig,l)=0. 313 313 314 314 ! ztcond-> ztnuc in test beneath to nucleate only when super saturation occurs(JL 2011) 315 IF ((zt(ig,l).LT.ztnuc(ig,l)).or.(zq(ig,l,i_n2ice).gt.1.E-10)) THEN316 pdtc(ig,l) = (ztcond(ig,l) - zt(ig,l))/subptimestep317 pdqc(ig,l,i_n2ice) = pdtc(ig,l)*ccond*g315 ! IF ((zt(ig,l).LT.ztnuc(ig,l)).or.(zq(ig,l,i_n2ice).gt.1.E-10)) THEN 316 ! pdtc(ig,l) = (ztcond(ig,l) - zt(ig,l))/subptimestep 317 ! pdqc(ig,l,i_n2ice) = pdtc(ig,l)*ccond*g 318 318 319 319 ! Case when the ice from above sublimes entirely 320 IF ((zq(ig,l,i_n2ice).lt.-pdqc(ig,l,i_n2ice)*subptimestep) &321 .AND. (zq(ig,l,i_n2ice).gt.0)) THEN322 323 pdqc(ig,l,i_n2ice) = -zq(ig,l,i_n2ice)/subptimestep324 pdtc(ig,l) =-zq(ig,l,i_n2ice)/(ccond*g*subptimestep)325 326 END IF320 ! IF ((zq(ig,l,i_n2ice).lt.-pdqc(ig,l,i_n2ice)*subptimestep) & 321 ! .AND. (zq(ig,l,i_n2ice).gt.0)) THEN 322 323 ! pdqc(ig,l,i_n2ice) = -zq(ig,l,i_n2ice)/subptimestep 324 ! pdtc(ig,l) =-zq(ig,l,i_n2ice)/(ccond*g*subptimestep) 325 326 ! END IF 327 327 328 328 ! Temperature and q after condensation 329 zt(ig,l) = zt(ig,l) + pdtc(ig,l) * subptimestep330 zq(ig,l,i_n2ice) = zq(ig,l,i_n2ice) + pdqc(ig,l,i_n2ice) * subptimestep331 END IF332 333 ENDDO334 ENDDO329 ! zt(ig,l) = zt(ig,l) + pdtc(ig,l) * subptimestep 330 ! zq(ig,l,i_n2ice) = zq(ig,l,i_n2ice) + pdqc(ig,l,i_n2ice) * subptimestep 331 ! END IF 332 333 ! ENDDO 334 ! ENDDO 335 335 336 ENDDO! end of subtimestep loop.336 ! ENDDO! end of subtimestep loop. 337 337 338 338 ! Computing global tendencies after the subtimestep. 339 DO l=1,nlayer340 DO ig=1,ngrid341 pdtc(ig,l) = &342 (zt(ig,l) - (pt(ig,l) + pdt(ig,l)*ptimestep))/ptimestep343 pdqc(ig,l,i_n2ice) = &344 (zq(ig,l,i_n2ice)-(pq(ig,l,i_n2ice)+pdq(ig,l,i_n2ice)*ptimestep))/ptimestep345 END DO346 END DO347 DO ig=1,ngrid348 zfallice(ig) = Mfallice(ig)/ptimestep349 END DO339 ! DO l=1,nlayer 340 ! DO ig=1,ngrid 341 ! pdtc(ig,l) = & 342 ! (zt(ig,l) - (pt(ig,l) + pdt(ig,l)*ptimestep))/ptimestep 343 ! pdqc(ig,l,i_n2ice) = & 344 ! (zq(ig,l,i_n2ice)-(pq(ig,l,i_n2ice)+pdq(ig,l,i_n2ice)*ptimestep))/ptimestep 345 ! END DO 346 ! END DO 347 ! DO ig=1,ngrid 348 ! zfallice(ig) = Mfallice(ig)/ptimestep 349 ! END DO 350 350 351 351 … … 357 357 ! Forecast of ground temperature ztsrf and frost temperature ztcondsol. 358 358 DO ig=1,ngrid 359 picen2(ig)=pqsurf(ig,i_n2ice) 359 360 ppn2=gfrac(igas_N2)*pplay(ig,1) 360 361 call get_tcond_n2(ppn2,ztcondsol(ig)) … … 430 431 DO ig=1,ngrid 431 432 432 IF(latitude(ig).LT.0.) THEN433 icap=2 ! Southern Hemisphere434 ELSE435 icap=1 ! Nortnern hemisphere436 ENDIF433 ! IF(latitude(ig).LT.0.) THEN 434 ! icap=2 ! Southern Hemisphere 435 ! ELSE 436 ! icap=1 ! Nortnern hemisphere 437 ! ENDIF 437 438 438 439 if(.not.picen2(ig).ge.0.) THEN … … 441 442 picen2(ig)=0. 442 443 endif 443 if (picen2(ig) .gt. 1.) then ! N2 Albedo condition changed to ~1 mm coverage. Change by MT2015.444 DO nw=1,L_NSPECTV445 albedo(ig,nw) = albedo_n2_ice_SPECTV(nw)446 ENDDO447 emisref(ig) = emisice(icap)448 else449 DO nw=1,L_NSPECTV450 albedo(ig,nw) = albedo_bareground(ig) ! Note : If you have some water, it will be taken into account in the "hydrol" routine.451 ENDDO452 emisref(ig) = emissiv453 pemisurf(ig) = emissiv454 end if444 ! if (picen2(ig) .gt. 1.) then ! N2 Albedo condition changed to ~1 mm coverage. Change by MT2015. 445 ! DO nw=1,L_NSPECTV 446 ! albedo(ig,nw) = albedo_n2_ice_SPECTV(nw) 447 ! ENDDO 448 ! emisref(ig) = emisice(icap) 449 ! else 450 ! DO nw=1,L_NSPECTV 451 ! albedo(ig,nw) = albedo_bareground(ig) ! Note : If you have some water, it will be taken into account in the "hydrol" routine. 452 ! ENDDO 453 ! emisref(ig) = emissiv 454 ! pemisurf(ig) = emissiv 455 ! end if 455 456 456 457 END DO … … 476 477 implicit none 477 478 478 real p, peff, tcond 479 real, parameter :: ptriple=518000.0 480 481 peff=p 482 483 if(peff.lt.ptriple) then 484 tcond = (-3167.8)/(log(.01*peff)-23.23) ! Fanale's formula. 485 else 486 tcond = 684.2-92.3*log(peff)+4.32*log(peff)**2 ! liquid-vapour transition (based on CRC handbook 2003 data) 487 endif 479 real p, tcond 480 481 IF (p.ge.0.529995) then 482 ! tcond Fray and Schmitt for N2 phase beta (T>35.6 K) FIT TB 483 tcond_n2 = (1.)/(0.01583606505-1.211938776e-3*log(7.963685594e-5*p*vmr)) 484 ELSE 485 ! tcond Fray and Schmitt for N2 phase alpha(T<35.6 K) FIT BT 486 tcond_n2 = (1.)/(1./35.6-1.089633028e-3*log(1.968275338*p*vmr)) 487 ENDIF 488 488 return 489 489 … … 513 513 peff=p/n2supsat 514 514 515 if(peff.lt.ptriple) then 516 tnuc = (-3167.8)/(log(.01*peff)-23.23) ! Fanale's formula 517 else 518 tnuc = 684.2-92.3*log(peff)+4.32*log(peff)**2 519 ! liquid-vapour transition (based on CRC handbook 2003 data) 520 endif 515 !! AF24: this code below is for CO2, needs to be udpated for N2 516 517 ! if(peff.lt.ptriple) then 518 ! tnuc = (-3167.8)/(log(.01*peff)-23.23) ! Fanale's formula 519 ! else 520 ! tnuc = 684.2-92.3*log(peff)+4.32*log(peff)**2 521 ! ! liquid-vapour transition (based on CRC handbook 2003 data) 522 ! endif 521 523 522 524 return
Note: See TracChangeset
for help on using the changeset viewer.
