Changeset 486 for trunk/LMDZ.GENERIC
- Timestamp:
- Dec 21, 2011, 10:09:07 AM (13 years ago)
- Location:
- trunk/LMDZ.GENERIC/libf/phystd
- Files:
-
- 1 deleted
- 11 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/LMDZ.GENERIC/libf/phystd/aeropacity.F90
r305 r486 267 267 do l=1,nlayer 268 268 zp=700./pplay(ig,l) 269 aerosol(ig,l, 1)=(dusttau/700.)*(pplev(ig,l)-pplev(ig,l+1)) &269 aerosol(ig,l,iaer)=(dusttau/700.)*(pplev(ig,l)-pplev(ig,l+1)) & 270 270 *max( exp(.03*(1.-max(zp,1.))) , 1.E-3 ) & 271 271 *QREFvis3d(ig,l,iaer) / QREFvis3d(ig,1,iaer) -
trunk/LMDZ.GENERIC/libf/phystd/callcorrk.F90
r471 r486 211 211 ! Effective radius and variance of the aerosols 212 212 213 ! CO2 ice: 214 DO l=1,nlayer 215 DO ig=1,ngrid 216 reffrad(ig,l,1) = 1.e-4 217 nueffrad(ig,l,1) = 0.1 213 do iaer=1,naerkind 218 214 ! these values will change once the microphysics gets to work 219 215 ! UNLESS tracer=.false., in which case we should be working with 220 216 ! a fixed aerosol layer, and be able to define reffrad in a 221 217 ! .def file. To be improved! 222 ENDDO 223 ENDDO 224 225 ! H2O ice: 226 if(naerkind.eq.2)then 227 DO l=1,nlayer 228 DO ig=1,ngrid 229 reffrad(ig,l,naerkind) = 1.e-5 230 nueffrad(ig,l,naerkind) = 0.1 231 ENDDO 232 ENDDO 233 endif 218 219 if(iaer.eq.1)then ! CO2 ice 220 do l=1,nlayer 221 do ig=1,ngrid 222 reffrad(ig,l,iaer) = 1.e-4 223 nueffrad(ig,l,iaer) = 0.1 224 enddo 225 enddo 226 endif 227 228 if(iaer.eq.2)then ! H2O ice 229 do l=1,nlayer 230 do ig=1,ngrid 231 reffrad(ig,l,iaer) = 1.e-5 232 nueffrad(ig,l,iaer) = 0.1 233 enddo 234 enddo 235 endif 236 237 if(iaer.eq.3)then ! dust 238 do l=1,nlayer 239 do ig=1,ngrid 240 reffrad(ig,l,iaer) = 1.e-5 241 nueffrad(ig,l,iaer) = 0.1 242 enddo 243 enddo 244 endif 245 246 if(iaer.gt.3)then 247 print*,'Error in callcorrk, naerkind is too high.' 248 print*,'The code still needs generalisation to arbitrary' 249 print*,'aerosol kinds and number.' 250 call abort 251 endif 252 253 enddo 234 254 235 255 print*, "callcorrk: Correlated-k data base folder:",trim(datadir) … … 248 268 Cmk= 0.01 * 1.0 / (g * mugaz * 1.672621e-27) ! q_main=1.0 assumed 249 269 250 251 270 if((igcm_h2o_vap.eq.0) .and. varactive)then 252 271 print*,'varactive in callcorrk but no h2o_vap tracer.' … … 268 287 enddo 269 288 enddo 270 271 289 272 290 if(kastprof)then … … 426 444 end do 427 445 428 429 446 ! longwave 430 447 DO nw=1,L_NSPECTI … … 465 482 end do 466 483 end do 467 468 484 469 485 ! test / correct for freaky s. s. albedo values … … 515 531 !tauaero(L_LEVELS+1,iaer) = 0. 516 532 end do 517 !print*,'Note changed tauaero BCs in callcorrk!'518 533 519 534 ! Albedo and emissivity … … 533 548 acosz=mu0(ig) ! cosine of sun incident angle 534 549 endif 535 536 550 537 551 !----------------------------------------------------------------------- … … 591 605 end do 592 606 593 594 595 607 !----------------------------------------------------------------------- 596 608 ! kcm mode only 597 609 if(kastprof)then 598 610 611 ! initial values equivalent to mugaz 599 612 DO l=1,nlayer 600 613 muvarrad(2*l) = mugaz … … 602 615 END DO 603 616 604 do k=1,L_LEVELS605 qvar(k) = 0.0606 end do607 print*,'ASSUMING qH2O=0 EVERYWHERE IN CALLCORRK!'617 !do k=1,L_LEVELS 618 ! qvar(k) = 0.0 619 !end do 620 !print*,'ASSUMING qH2O=0 EVERYWHERE IN CALLCORRK!' 608 621 endif 609 622 … … 633 646 qvar(2*nlayermx+1)=qsurf(ig,i_var)*muvar(ig,1)/mH2O 634 647 635 636 endif 637 648 endif 638 649 639 650 ! Keep values inside limits for which we have radiative transfer coefficients … … 678 689 tmid(L_LEVELS) = tlevrad(L_LEVELS) 679 690 680 681 691 ! test for out-of-bounds pressure 682 692 if(plevrad(3).lt.pgasmin)then … … 695 705 print*,'Minimum temperature is outside the radiative' 696 706 print*,'transfer kmatrix bounds, exiting.' 697 print*,'WARNING, OVERRIDING FOR TEST'698 !call abort707 !print*,'WARNING, OVERRIDING FOR TEST' 708 call abort 699 709 elseif(tlevrad(k).gt.tgasmax)then 700 710 print*,'Maximum temperature is outside the radiative' 701 711 print*,'transfer kmatrix bounds, exiting.' 702 print*,'WARNING, OVERRIDING FOR TEST' 703 !print*, 'T=',pt 704 !call abort 712 !print*,'WARNING, OVERRIDING FOR TEST' 713 call abort 705 714 endif 706 715 enddo … … 735 744 tmid,pmid,taugsurf,qvar,muvarrad) 736 745 737 738 746 call sfluxv(dtauv,tauv,taucumv,albv,dwnv,wbarv,cosbv, & 739 747 acosz,stel_fract,gweight, & 740 748 nfluxtopv,nfluxoutv_nu,nfluxgndv_nu, & 741 !acosz,stel_fract,gweight,nfluxtopv,nfluxgndv_nu, &742 749 fmnetv,fluxupv,fluxdnv,fzerov,taugsurf) 743 750 … … 753 760 end if 754 761 755 756 757 758 762 !----------------------------------------------------------------------- 759 763 ! Longwave … … 846 850 if(specOLR)then 847 851 if(ngrid.ne.1)then 848 !call writediagspecIR(ngrid,"OLR3D","OLR(lon,lat,band)","W m^-2",3,OLR_nu)849 !call writediagspecVI(ngrid,"OSR3D","OSR(lon,lat,band)","W m^-2",3,OSR_nu)852 call writediagspecIR(ngrid,"OLR3D","OLR(lon,lat,band)","W m^-2",3,OLR_nu) 853 call writediagspecVI(ngrid,"OSR3D","OSR(lon,lat,band)","W m^-2",3,OSR_nu) 850 854 endif 851 855 endif … … 892 896 endif 893 897 894 ! !!see physiq.F for explanations about CLFvarying. This is temporary.898 ! see physiq.F for explanations about CLFvarying. This is temporary. 895 899 if (lastcall .and. .not.CLFvarying) then 896 900 IF( ALLOCATED( gasi ) ) DEALLOCATE( gasi ) -
trunk/LMDZ.GENERIC/libf/phystd/callkeys.h
r305 r486 38 38 & , hydrology & 39 39 & , sourceevol & 40 & , icetstep & 40 41 & , albedosnow & 41 42 & , maxicethick & … … 101 102 real tau_relax 102 103 real cloudlvl 104 real icetstep 105 -
trunk/LMDZ.GENERIC/libf/phystd/callsedim.F
r253 r486 53 53 INTEGER l,ig, iq 54 54 55 real zqi(ngridmx,nlayermx ) ! to locally store tracers55 real zqi(ngridmx,nlayermx,nqmx) ! to locally store tracers 56 56 real masse (ngridmx,nlayermx) ! Layer mass (kg.m-2) 57 57 real epaisseur (ngridmx,nlayermx) ! Layer thickness (m) … … 79 79 firstcall=.false. 80 80 ENDIF ! of IF (firstcall) 81 81 82 82 !======================================================================= 83 83 ! Preliminary calculation of the layer characteristics … … 91 91 end do 92 92 93 iq=igcm_h2o_ice 93 !====================================================================== 94 ! Calculate the transport due to sedimentation for each tracer 95 ! [This has been rearranged by L. Kerber to allow the sedimentation 96 ! of general tracers.] 97 98 zqi(1:ngrid,1:nlay,1:nqmx) = 0. 99 do iq=1,nq 100 if((radius(iq).gt.1.e-9).and.(iq.ne.igcm_co2_ice)) then 101 ! (no sedim for gases; co2_ice sedim is done elsewhere) 94 102 95 ! The value of q is updated after the other parameterisations103 ! store locally updated tracers 96 104 97 do l=1,nlay 105 do l=1,nlay 98 106 do ig=1,ngrid 99 ! store locally updated tracers 100 zqi(ig,l)=pq(ig,l,iq)+pdqfi(ig,l,iq)*ptimestep 107 zqi(ig,l,iq)=pq(ig,l,iq)+pdqfi(ig,l,iq)*ptimestep 101 108 enddo 102 109 enddo ! of do l=1,nlay 110 111 !====================================================================== 112 ! Sedimentation 113 !====================================================================== 114 ! Water 115 if (water.and.(iq.eq.igcm_h2o_ice)) then 116 call newsedim(ngrid,nlay,ngrid*nlay,ptimestep, 117 & pplev,masse,epaisseur,pt,rfall,rho_q(iq),zqi,wq) 103 118 104 105 106 !======================================================================= 107 ! Calculate the transport due to sedimentation for each tracer 108 109 call newsedim(ngrid,nlay,ngrid*nlay,ptimestep, 110 & pplev,masse,epaisseur,pt,rfall,rho_q(iq),zqi,wq) 119 ! General Case 120 else 121 call newsedim(ngrid,nlay,1,ptimestep, 122 & pplev,masse,epaisseur,pt,radius(iq),rho_q(iq), 123 & zqi,wq) 124 endif 111 125 112 126 !======================================================================= 113 127 ! Calculate the tendencies 128 !====================================================================== 114 129 115 130 do ig=1,ngrid 116 131 ! Ehouarn: with new way of tracking tracers by name, this is simply 117 pdqs_sed(ig,iq) =wq(ig,1)/ptimestep132 pdqs_sed(ig,iq) = wq(ig,1)/ptimestep 118 133 end do 119 134 120 135 DO l = 1, nlay 121 136 DO ig=1,ngrid 122 pdqsed(ig,l,iq)=(zqi(ig,l )-123 $(pq(ig,l,iq) + pdqfi(ig,l,iq)*ptimestep))/ptimestep137 pdqsed(ig,l,iq)=(zqi(ig,l,iq)- 138 & (pq(ig,l,iq) + pdqfi(ig,l,iq)*ptimestep))/ptimestep 124 139 ENDDO 125 140 ENDDO 126 141 endif ! of no gases no co2_ice 142 enddo ! of do iq=1,nq 127 143 return 128 144 end 129 -
trunk/LMDZ.GENERIC/libf/phystd/condense_cloud.F90
r471 r486 187 187 enddo 188 188 189 write(*,*) "condense_c o2cloud: i_co2ice=",i_co2ice189 write(*,*) "condense_cloud: i_co2ice=",i_co2ice 190 190 191 191 if((i_co2ice.lt.1))then 192 print*,'In condens_co2cloud but no CO2 ice tracer, exiting.' 192 print*,'In condens_cloud but no CO2 ice tracer, exiting.' 193 print*,'Still need generalisation to arbitrary species!' 193 194 stop 194 195 endif 195 196 196 197 ccond=cpp/(g*latcond) 197 print*,'In condens_c o2cloud: ccond=',ccond,' latcond=',latcond198 print*,'In condens_cloud: ccond=',ccond,' latcond=',latcond 198 199 199 200 ! Prepare special treatment if gas is not pure CO2 -
trunk/LMDZ.GENERIC/libf/phystd/inifis.F
r374 r486 472 472 write(*,*) " sourceevol = ",sourceevol 473 473 474 write(*,*) "Ice evolution timestep ?" 475 icetstep=100.0 ! default value 476 call getin("icetstep",icetstep) 477 write(*,*) " icetstep = ",icetstep 478 474 479 write(*,*) "Snow albedo ?" 475 480 albedosnow=0.5 ! default value -
trunk/LMDZ.GENERIC/libf/phystd/newsedim.F
r253 r486 79 79 PRINT*,'ngrid =',ngrid 80 80 PRINT*,'ngridmx =',ngridmx 81 STOP 81 STOP 82 82 ENDIF 83 83 firstcall=.false. … … 111 111 c (stokes law corrected for low pressure by the Cunningham 112 112 c slip-flow correction according to Rossow (Icarus 36, 1-50, 1978) 113 113 114 114 do l=1,nlay 115 115 do ig=1, ngrid … … 121 121 endif 122 122 123 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!124 ! TEMPORARY MODIF BY RDW !!!!125 !rfall=5.e-6126 if(rfall.lt.1.e-7)then127 rfall=1.e-7128 endif129 !if(rfall.gt.5.e-5)then130 ! rfall=5.e-5131 !endif132 133 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!134 135 123 vstokes(ig,l) = b * rho * rfall**2 * 136 124 & (1 + 1.333* ( a*pt(ig,l)/pplev(ig,l) )/rfall) … … 147 135 c va traverser le niveau intercouche l : "dztop" est sa hauteur 148 136 c au dessus de l (m), "ptop" est sa pression (Pa)) 149 150 137 do l=1,nlay 151 138 do ig=1, ngrid … … 154 141 Ep=0 155 142 k=0 156 143 w(ig,l) = 0. !! JF+AS ajout initialisation (LK MARS) 157 144 c ************************************************************** 158 145 c Simple Method 159 w(ig,l) = 160 & (1- exp(-dztop*g/(r*pt(ig,l))))*pplev(ig,l) / g 161 cc write(*,*) 'OK simple method l,w =', l, w(ig,l) 162 cc write(*,*) 'OK simple method dztop =', dztop 163 c ************************************************************** 146 cc w(ig,l) = 147 cc & (1- exp(-dztop*g/(r*pt(ig,l))))*pplev(ig,l) / g 148 cc write(*,*) 'OK simple method l,w =', l, w(ig,l) 149 cc write(*,*) 'OK simple method dztop =', dztop 150 w(ig,l) = 1. - exp(-dztop*g/(r*pt(ig,l))) 151 !!! Diagnostic: JF. Fix: AS. Date: 05/11 152 !!! Probleme arrondi avec la quantite ci-dessus 153 !!! ---> vaut 0 pour -dztop*g/(r*pt(ig,l)) trop petit 154 !!! ---> dans ce cas on utilise le developpement limite ! 155 !!! ---> exp(-x) = 1 - x lorsque x --> 0 avec une erreur de x^2 / 2 156 157 IF ( w(ig,l) .eq. 0. ) THEN 158 w(ig,l) = ( dztop*g/(r*pt(ig,l)) ) * pplev(ig,l) / g 159 ELSE 160 w(ig,l) = w(ig,l) * pplev(ig,l) / g 161 ENDIF 162 ! LK borrowed simple method from Mars model (AS/JF) 163 164 !************************************************************** 164 165 cccc Complex method : 165 166 if (dztop.gt.epaisseur(ig,l)) then 166 167 cccc Cas ou on "epuise" la couche l : On calcule le flux 167 168 cccc Venant de dessus en tenant compte de la variation de Vstokes … … 176 177 enddo 177 178 Ep = Ep - epaisseur(ig,l+k) 178 end if 179 ptop=pplev(ig,l+k)*exp(-(dztop-Ep)*g/(r*pt(ig,l+k))) 180 w(ig,l) = (pplev(ig,l) -Ptop)/g 179 ! ptop=pplev(ig,l+k)*exp(-(dztop-Ep)*g/(r*pt(ig,l+k))) 180 ptop=exp(-(dztop-Ep)*g/(r*pt(ig,l+k))) 181 IF ( ptop .eq. 1. ) THEN 182 !PRINT*, 'newsedim: exposant trop petit ', ig, l 183 ptop=pplev(ig,l+k) * ( 1. - (dztop-Ep)*g/(r*pt(ig,l+k))) 184 ELSE 185 ptop=pplev(ig,l+k) * ptop 186 ENDIF 187 188 w(ig,l) = (pplev(ig,l) - ptop)/g 189 190 endif !!! complex method 181 191 c 182 192 cc write(*,*) 'OK new method l,w =', l, w(ig,l) … … 188 198 cc if(l.eq.8)write(*,*)'l=8,k, w',k, w(1,l) 189 199 c ************************************************************** 200 190 201 end do 191 202 end do … … 195 206 c & wq(1,6),wq(1,7),pqi(1,6) 196 207 197 198 208 RETURN 199 209 END 200 -
trunk/LMDZ.GENERIC/libf/phystd/physiq.F90
r471 r486 1917 1917 qsurf_hist(ig,igcm_h2o_ice) = & 1918 1918 !qsurf_hist(ig,igcm_h2o_ice) + delta_ice*100.0 1919 qsurf_hist(ig,igcm_h2o_ice) + delta_ice* 10.01919 qsurf_hist(ig,igcm_h2o_ice) + delta_ice*icetstep 1920 1920 1921 1921 ! if ice has gone -ve, set to zero -
trunk/LMDZ.GENERIC/libf/phystd/radinc_h.F90
r470 r486 6 6 #include "bands.h" 7 7 8 !====================================================================== C8 !====================================================================== 9 9 ! 10 10 ! RADINC.H RADiation INCludes … … 13 13 ! number of spectral intervals. . . 14 14 ! 15 !====================================================================== C15 !====================================================================== 16 16 17 17 ! RADIATION parameters … … 60 60 integer, parameter :: L_NLEVRAD = llm+1 61 61 62 ! !!! THIS IS SET INsugas_corrk63 ! !!! [use ofallocatable arrays] -- AS 12/201162 ! These are set in sugas_corrk 63 ! [uses allocatable arrays] -- AS 12/2011 64 64 integer :: L_NPREF, L_NTREF, L_REFVAR, L_PINT 65 65 66 !!! THIS ONE DOES NOT CHANGE MUCH... IT CAN BE REGARDED AS A CONSTANT.67 66 integer, parameter :: L_NGAUSS = 17 68 67 … … 72 71 integer, parameter :: NAERKIND = 2 73 72 real, parameter :: L_TAUMAX = 35 74 !integer, parameter :: L_TAUMAX = 3575 !integer, parameter :: L_TAUMAX = 4076 73 77 74 ! For Planck function integration: -
trunk/LMDZ.GENERIC/libf/phystd/rcm1d.F
r305 r486 112 112 REAL cloudfrac(ngridmx,nlayermx) 113 113 REAL hice(ngridmx),totcloudfrac(ngridmx) 114 115 116 c=======================================================================117 114 118 115 c======================================================================= … … 613 610 614 611 DO idt=1,ndt 615 IF (idt.eq.ndt -day_step-1) then !test612 IF (idt.eq.ndt) then !test 616 613 lastcall=.true. 617 614 call stellarlong(day*1.0,zls) -
trunk/LMDZ.GENERIC/libf/phystd/sfluxi.F
r253 r486 31 31 real*8 FZEROI(L_NSPECTI) 32 32 real*8 taugsurf(L_NSPECTI,L_NGAUSS-1), fzero 33 34 ! real*8 BSURFtest ! by RW for test35 33 36 34 real*8 fup_tmp(L_NSPECTI),fdn_tmp(L_NSPECTI) … … 73 71 ! NTS = TSURF*10.0D0-499 74 72 ! NTT = TTOP *10.0D0-499 75 76 ! BSURFtest=0.077 73 78 74 DO 501 NW=1,L_NSPECTI
Note: See TracChangeset
for help on using the changeset viewer.