Changeset 1016 for trunk/LMDZ.GENERIC
- Timestamp:
- Aug 7, 2013, 4:21:01 PM (11 years ago)
- Location:
- trunk/LMDZ.GENERIC
- Files:
-
- 13 edited
-
README (modified) (2 diffs)
-
libf/phystd/calc_rayleigh.F90 (modified) (11 diffs)
-
libf/phystd/callcorrk.F90 (modified) (6 diffs)
-
libf/phystd/gfluxi.F (modified) (1 diff)
-
libf/phystd/largescale.F90 (modified) (11 diffs)
-
libf/phystd/moistadj.F90 (modified) (6 diffs)
-
libf/phystd/optci.F90 (modified) (1 diff)
-
libf/phystd/optcv.F90 (modified) (1 diff)
-
libf/phystd/physiq.F90 (modified) (11 diffs)
-
libf/phystd/radcommon_h.F90 (modified) (2 diffs)
-
libf/phystd/rain.F90 (modified) (2 diffs)
-
libf/phystd/turbdiff.F90 (modified) (1 diff)
-
libf/phystd/watercommon_h.F90 (modified) (1 diff)
Legend:
- Unmodified
- Added
- Removed
-
trunk/LMDZ.GENERIC/README
r1006 r1016 940 940 == 15/05/2013 == JL 941 941 - correction in radiative scheme to enforce double precision 942 - corrects calculation of ISR 942 943 943 944 == 22/05/2013 == EM … … 991 992 compute the sponge quenching analytically rather than via Forward Euler 992 993 approximation. 994 995 == 07/08/2013 == JL 996 - Water cycle in double precision (largescale+moistadj) 997 - Improved wate rayleigh. 998 - First step for rayleigh with variable species. Now, just need to change optcv. 999 - changed some interpolation indices in callcorrk to limit dependency of OLR on the number of layers -
trunk/LMDZ.GENERIC/libf/phystd/calc_rayleigh.F90
r995 r1016 13 13 ! ------- 14 14 ! Robin Wordsworth (2010) 15 ! Jeremy Leconte (2012): Added option for variable gas. Improved water rayleigh (Bucholtz 1995). 15 16 ! 16 17 ! Called by … … 25 26 26 27 use radinc_h, only: L_NSPECTV 27 use radcommon_h, only: WAVEV, BWNV, DWNV, tstellar, tauray, scalep28 use radcommon_h, only: WAVEV, BWNV, DWNV, tstellar, tauray, taurayvar, scalep 28 29 use gases_h 29 30 … … 39 40 40 41 logical typeknown 41 real*8 tauvar,tau sum,tauwei,bwidth,bstart42 real*8 tauvar,tauvarvar,tausum,tausumvar,tauwei,tauweivar,bwidth,bstart 42 43 double precision df 43 44 … … 54 55 do igas=1,ngasmx 55 56 if(igas.eq.vgas)then 56 print*,' Ignoring ',trim(gnom(igas)),' in Rayleigh scattering '// &57 'as it is variable.' 58 elseif(gfrac(igas).lt.5.e-2)then57 print*,'variable gas is ',trim(gnom(igas)),' in Rayleigh scattering ' 58 endif 59 if((igas/=vgas).and.(gfrac(igas).lt.5.e-2))then 59 60 print*,'Ignoring ',trim(gnom(igas)),' in Rayleigh scattering '// & 60 61 'as its mixing ratio is less than 0.05.' … … 68 69 tauconsti(igas) = (9.81/g)*8.569E-3*scalep/(P0/93.0) 69 70 elseif(igas.eq.igas_H2O)then 70 tauconsti(igas) = (10.0/g)*9.22E-3*scalep/101325.0 71 ! tauconsti(igas) = (10.0/g)*9.22E-3*scalep/101325.0 72 tauconsti(igas) = 1.98017E-10/(g*mugaz*100.) 71 73 elseif(igas.eq.igas_H2)then 72 74 tauconsti(igas) = (10.0/g)*0.0241*scalep/101325.0 … … 83 85 endif 84 86 85 if( gfrac(igas).eq.1.0)then87 if((gfrac(igas).eq.1.0).and.(vgas.eq.0))then 86 88 print*,'Rayleigh scattering is for a pure ',trim(gnom(igas)),' atmosphere.' 87 89 typeknown=.true. … … 101 103 tausum = 0.0 102 104 tauwei = 0.0 105 tausumvar = 0.0 106 tauweivar = 0.0 103 107 bstart = 10000.0/BWNV(N+1) 104 108 bwidth = (10000.0/BWNV(N)) - (10000.0/BWNV(N+1)) … … 107 111 108 112 tauvar=0.0 113 tauvarvar=0.0 109 114 do igas=1,ngasmx 110 if((igas .eq.vgas).or.(gfrac(igas).lt.5.e-2))then115 if((igas/=vgas).and.(gfrac(igas).lt.5.e-2))then 111 116 ! ignore variable gas in Rayleigh calculation 112 117 ! ignore gases of mixing ratio < 0.05 in Rayleigh calculation … … 118 123 tauvari(igas) = (1.0+0.0113/wl**2+0.00013/wl**4)/wl**4 119 124 elseif(igas.eq.igas_H2O)then 120 tauvari(igas) = 1.0/wl**4 ! to be improved... 125 ! tauvari(igas) = 1.0/wl**4 ! to be improved... 126 tauvari(igas) = (5.7918E6/(2.38E2-1/wl**2)+1.679E5/(57.36E0-1/wl**2))**2/wl**4 121 127 elseif(igas.eq.igas_H2)then 122 128 tauvari(igas) = 1.0/wl**4 … … 130 136 endif 131 137 132 tauvar=tauvar+tauconsti(igas)*tauvari(igas)*gfrac(igas) 138 if(igas.eq.vgas) then 139 tauvarvar=tauvarvar+tauconsti(igas)*tauvari(igas) 140 tauvar=tauvar+0.0*0.0*gfrac(igas) 141 else 142 tauvar=tauvar+tauconsti(igas)*tauvari(igas)*gfrac(igas) 143 endif 133 144 134 145 enddo … … 138 149 tauwei=tauwei+df 139 150 tausum=tausum+tauvar*df 151 tauweivar=tauweivar+df 152 tausumvar=tausumvar+tauvarvar*df 140 153 141 154 enddo 142 155 TAURAY(N)=tausum/tauwei 156 TAURAYVAR(N)=tausumvar/tauweivar 143 157 ! we multiply by scalep here (100) because plev, which is used in optcv, 144 158 ! is in units of mBar, so we need to convert. -
trunk/LMDZ.GENERIC/libf/phystd/callcorrk.F90
r962 r1016 460 460 endif 461 461 462 !!! JL13: in the following, I changed some indices in the interpolations so that the model rsults are less dependent on the number of layers 463 !!! the older verions are commented with the commetn !JL13index 464 465 462 466 !----------------------------------------------------------------------- 463 467 ! Water vapour (to be generalised for other gases eventually) … … 468 472 do l=1,nlayer 469 473 qvar(2*l) = pq(ig,nlayer+1-l,i_var) 470 qvar(2*l+1) = (pq(ig,nlayer+1-l,i_var)+pq(ig,max(nlayer-l,1),i_var))/2 471 ! Average approximation as for temperature... 474 qvar(2*l+1) = pq(ig,nlayer+1-l,i_var) 475 !JL13index qvar(2*l+1) = (pq(ig,nlayer+1-l,i_var)+pq(ig,max(nlayer-l,1),i_var))/2 476 !JL13index ! Average approximation as for temperature... 472 477 end do 473 478 qvar(1)=qvar(2) … … 497 502 if(RH.lt.0.0) RH=0.0 498 503 499 ptemp = pplev(ig,1) 500 Ttemp = tsurf(ig) 501 call watersat(Ttemp,ptemp,qsat) 502 503 !qvar(2*nlayermx+1)=qsat ! fully saturated everywhere 504 ! ptemp = pplev(ig,1) 505 ! Ttemp = tsurf(ig) 506 ! call watersat(Ttemp,ptemp,qsat) 507 504 508 qvar(2*nlayermx+1)= RH * qsat ! ~realistic profile (e.g. 80% saturation at ground) 505 !qvar=0.005 ! completely constant profile (JL) 506 509 507 510 else 508 511 do k=1,L_LEVELS … … 528 531 END DO 529 532 530 !do k=1,L_LEVELS 531 ! qvar(k) = 0.0 532 !end do 533 !print*,'ASSUMING qH2O=0 EVERYWHERE IN CALLCORRK!' 534 endif 535 536 537 if(kastprof.and.(ngasmx.gt.1))then 538 533 if(ngasmx.gt.1)then 534 535 DO l=1,nlayer 536 muvarrad(2*l) = muvar(ig,nlayer+2-l) 537 muvarrad(2*l+1) = (muvar(ig,nlayer+2-l) + & 538 muvar(ig,max(nlayer+1-l,1)))/2 539 END DO 540 541 muvarrad(1) = muvarrad(2) 542 muvarrad(2*nlayermx+1)=muvar(ig,1) 543 544 print*,'Recalculating qvar with VARIABLE epsi for kastprof' 545 print*,'Assumes that the variable gas is H2O!!!' 546 print*,'Assumes that there is only one tracer' 547 !i_var=igcm_h2o_vap 548 i_var=1 549 if(nq.gt.1)then 550 print*,'Need 1 tracer only to run kcm1d.e' 551 stop 552 endif 553 do l=1,nlayer 554 vtmp(l)=pq(ig,l,i_var)/(epsi+pq(ig,l,i_var)*(1.-epsi)) 555 !vtmp(l)=pq(ig,l,i_var)*muvar(ig,l+1)/mH2O !JL to be changed 556 end do 557 558 do l=1,nlayer 559 qvar(2*l) = vtmp(nlayer+1-l) 560 qvar(2*l+1) = vtmp(nlayer+1-l) 561 ! qvar(2*l+1) = ( vtmp(nlayer+1-l) + vtmp(max(nlayer-l,1)) )/2 562 end do 563 qvar(1)=qvar(2) 564 565 print*,'Warning: reducing qvar in callcorrk.' 566 print*,'Temperature profile no longer consistent ', & 567 'with saturated H2O. qsat=',satval 568 do k=1,L_LEVELS 569 qvar(k) = qvar(k)*satval 570 end do 571 572 endif 573 else ! if kastprof 539 574 DO l=1,nlayer 540 575 muvarrad(2*l) = muvar(ig,nlayer+2-l) 541 muvarrad(2*l+1) = (muvar(ig,nlayer+2-l) + & 542 muvar(ig,max(nlayer+1-l,1)))/2 576 muvarrad(2*l+1) = (muvar(ig,nlayer+2-l)+muvar(ig,max(nlayer+1-l,1)))/2 543 577 END DO 544 578 545 579 muvarrad(1) = muvarrad(2) 546 muvarrad(2*nlayermx+1)=muvar(ig,1) 547 548 print*,'Recalculating qvar with VARIABLE epsi for kastprof' 549 print*,'Assumes that the variable gas is H2O!!!' 550 print*,'Assumes that there is only one tracer' 551 !i_var=igcm_h2o_vap 552 i_var=1 553 if(nq.gt.1)then 554 print*,'Need 1 tracer only to run kcm1d.e' 555 stop 556 endif 557 do l=1,nlayer 558 vtmp(l)=pq(ig,l,i_var)*muvar(ig,l+1)/mH2O 559 end do 560 561 do l=1,nlayer 562 qvar(2*l) = vtmp(nlayer+1-l) 563 qvar(2*l+1) = ( vtmp(nlayer+1-l) + vtmp(max(nlayer-l,1)) )/2 564 end do 565 qvar(1)=qvar(2) 566 567 print*,'Warning: reducing qvar in callcorrk.' 568 print*,'Temperature profile no longer consistent ', & 569 'with saturated H2O.' 570 do k=1,L_LEVELS 571 qvar(k) = qvar(k)*satval 572 end do 573 580 muvarrad(2*nlayermx+1)=muvar(ig,1) 574 581 endif 575 582 576 583 ! Keep values inside limits for which we have radiative transfer coefficients 577 584 if(L_REFVAR.gt.1)then ! there was a bug here! … … 607 614 608 615 DO l=1,L_NLAYRAD-1 609 tmid(2*l+1) = tlevrad(2*l+1) 610 tmid(2*l+2) = tlevrad(2*l+1) 611 pmid(2*l+1) = plevrad(2*l+1) 612 pmid(2*l+2) = plevrad(2*l+1) 616 tmid(2*l+1) = tlevrad(2*l) 617 tmid(2*l+2) = tlevrad(2*l) 618 pmid(2*l+1) = plevrad(2*l) 619 pmid(2*l+2) = plevrad(2*l) 620 !JL13index tmid(2*l+1) = tlevrad(2*l+1) 621 !JL13index tmid(2*l+2) = tlevrad(2*l+1) 622 !JL13index pmid(2*l+1) = plevrad(2*l+1) 623 !JL13index pmid(2*l+2) = plevrad(2*l+1) 613 624 END DO 614 pmid(L_LEVELS) = plevrad(L_LEVELS) 615 tmid(L_LEVELS) = tlevrad(L_LEVELS) 625 pmid(L_LEVELS) = plevrad(L_LEVELS-1) 626 tmid(L_LEVELS) = tlevrad(L_LEVELS-1) 627 !JL13index pmid(L_LEVELS) = plevrad(L_LEVELS) 628 !JL13index tmid(L_LEVELS) = tlevrad(L_LEVELS) 616 629 617 630 ! test for out-of-bounds pressure … … 635 648 !print*,'WARNING, OVERRIDING FOR TEST' 636 649 call abort 650 !tlevrad(k)=tgasmin 637 651 elseif(tlevrad(k).gt.tgasmax)then 638 652 print*,'Maximum temperature is outside the radiative' 639 653 print*,'transfer kmatrix bounds, exiting.' 640 654 print*,'level,grid,T',k,ig,tlevrad(k) 641 print*,'WARNING, OVERRIDING FOR TEST' 642 !call abort 643 tlevrad(k)=tgasmax 655 !print*,'WARNING, OVERRIDING FOR TEST' 656 call abort 657 !tlevrad(k)=tgasmax 658 endif 659 enddo 660 do k=1,L_NLAYRAD+1 661 if(tmid(k).lt.tgasmin)then 662 print*,'Minimum temperature is outside the radiative' 663 print*,'transfer kmatrix bounds, exiting.' 664 print*,"k=",k," tlevrad(k)=",tlevrad(k) 665 print*,"tgasmin=",tgasmin 666 !print*,'WARNING, OVERRIDING FOR TEST' 667 call abort 668 !tmid(k)=tgasmin 669 elseif(tmid(k).gt.tgasmax)then 670 print*,'Maximum temperature is outside the radiative' 671 print*,'transfer kmatrix bounds, exiting.' 672 print*,'level,grid,T',k,ig,tlevrad(k) 673 !print*,'WARNING, OVERRIDING FOR TEST' 674 call abort 675 !tmid(k)=tgasmax 644 676 endif 645 677 enddo -
trunk/LMDZ.GENERIC/libf/phystd/gfluxi.F
r995 r1016 105 105 ! -- and stabilizes integrations 106 106 NT = int(TLEV(2*L+1)*NTfac) - NTstar+1 107 108 107 !! For deep, opaque, thick first layers (e.g. Saturn) 109 108 !! what is below works much better, not unstable, ... -
trunk/LMDZ.GENERIC/libf/phystd/largescale.F90
r875 r1016 2 2 pdt, pdq, pdtlsc, pdqvaplsc, pdqliqlsc, rneb) 3 3 4 5 ! to use 'getin' 6 use ioipsl_getincom 4 7 use watercommon_h, only : RLVTT, RCPD, RVTMP2, & 5 T_h2O_ice_clouds,T_h2O_ice_liq,Psat_water ,Lcpdqsat_water8 T_h2O_ice_clouds,T_h2O_ice_liq,Psat_waterDP,Lcpdqsat_waterDP 6 9 USE tracer_h 7 10 IMPLICIT none … … 33 36 REAL pplay(ngrid,nlayermx) ! pression au milieu de couche 34 37 REAL pt(ngrid,nlayermx) ! temperature (K) 35 realpq(ngrid,nlayermx,nq) ! tracer mixing ratio (kg/kg)38 REAL pq(ngrid,nlayermx,nq) ! tracer mixing ratio (kg/kg) 36 39 REAL pdt(ngrid,nlayermx) ! physical temperature tenedency (K/s) 37 40 REAL pdq(ngrid,nlayermx,nq)! physical tracer tenedency (K/s) … … 43 46 44 47 ! Options du programme 45 REAL ratqs ! determine largeur de la distribution de vapeur 46 PARAMETER (ratqs=0.2) 48 REAL, SAVE :: ratqs ! determine largeur de la distribution de vapeur 47 49 48 50 ! Variables locales … … 50 52 EXTERNAL CBRT 51 53 INTEGER i, k , nn 52 INTEGER,PARAMETER :: nitermax= 200053 REAL,PARAMETER :: alpha=.5,qthreshold=1.e-654 INTEGER,PARAMETER :: nitermax=5000 55 DOUBLE PRECISION,PARAMETER :: alpha=.1,qthreshold=1.d-8 54 56 ! JL13: if "careful, T<Tmin in psat water" appears often, you may want to stabilise the model by 55 57 ! decreasing alpha and increasing nitermax accordingly 56 REALzt(ngrid), zq(ngrid)57 REALzcond(ngrid),zcond_iter58 REALzdelq(ngrid)59 REALzqs(ngrid), zdqs(ngrid)60 REAL psat_tmp,dlnpsat_tmp58 DOUBLE PRECISION zt(ngrid), zq(ngrid) 59 DOUBLE PRECISION zcond(ngrid),zcond_iter 60 DOUBLE PRECISION zdelq(ngrid) 61 DOUBLE PRECISION zqs(ngrid), zdqs(ngrid) 62 DOUBLE PRECISION local_p,psat_tmp,dlnpsat_tmp,Lcp 61 63 62 64 ! evaporation calculations … … 65 67 REAL tevap(ngrid,nlayermx) 66 68 67 REAL zcor(ngrid), zdelta(ngrid), zcvm5(ngrid) 68 REAL zx_q(ngrid) 69 REAL Nmix_local,zfice 69 DOUBLE PRECISION zx_q(ngrid) 70 LOGICAL,SAVE :: firstcall=.true. 71 72 73 IF (firstcall) THEN 74 75 write(*,*) "value for ratqs? " 76 ratqs=0.2 ! default value 77 call getin("ratqs",ratqs) 78 write(*,*) " ratqs = ",ratqs 79 80 firstcall = .false. 81 ENDIF 70 82 71 83 ! GCM -----> subroutine variables, initialisation of outputs … … 75 87 pdqliqlsc(1:ngrid,1:nlayermx) = 0.0 76 88 rneb(1:ngrid,1:nlayermx) = 0.0 89 Lcp=RLVTT/RCPD 77 90 78 91 … … 93 106 DO i = 1, ngrid 94 107 108 local_p=pplay(i,k) 95 109 if(zt(i).le.15.) then 96 110 print*,'in lsc',i,k,zt(i) 97 111 ! zt(i)=15. ! check too low temperatures 98 112 endif 99 call Psat_water (zt(i),pplay(i,k),psat_tmp,zqs(i))113 call Psat_waterDP(zt(i),local_p,psat_tmp,zqs(i)) 100 114 101 zdelq(i) = MAX(MIN(ratqs * zq(i),1.-zq(i)),1. e-12)115 zdelq(i) = MAX(MIN(ratqs * zq(i),1.-zq(i)),1.d-12) 102 116 rneb(i,k) = (zq(i)+zdelq(i)-zqs(i)) / (2.0*zdelq(i)) 103 ! print*,zq(i),zdelq(i),zqs(i),rneb(i,k)104 117 if (rneb(i,k).lt.0.) then !no clouds 105 118 … … 107 120 zcond(i)=0. 108 121 109 else if (rneb(i,k).gt.1.) then !complete cloud cover, we start without evaporating 110 122 else if ((rneb(i,k).gt.0.99).or.(ratqs.lt.1.e-6)) then !complete cloud cover, we start without evaporating 111 123 rneb(i,k)=1. 112 124 zt(i)=pt(i,k)+pdt(i,k)*ptimestep … … 114 126 dqevap(i,k)=0. 115 127 ! iterative process to stabilize the scheme when large water amounts JL12 116 zcond(i) = 0.0 128 zcond(i) = 0.0d0 117 129 Do nn=1,nitermax 118 call Psat_water (zt(i),pplay(i,k),psat_tmp,zqs(i))119 call Lcpdqsat_water (zt(i),pplay(i,k),psat_tmp,zqs(i),zdqs(i),dlnpsat_tmp)120 zcond_iter = alpha*(zx_q(i)-zqs(i))/(1. +zdqs(i))130 call Psat_waterDP(zt(i),local_p,psat_tmp,zqs(i)) 131 call Lcpdqsat_waterDP(zt(i),local_p,psat_tmp,zqs(i),zdqs(i),dlnpsat_tmp) 132 zcond_iter = alpha*(zx_q(i)-zqs(i))/(1.d0+zdqs(i)) 121 133 !zcond can be negative here 122 134 zx_q(i) = zx_q(i) - zcond_iter 123 135 zcond(i) = zcond(i) + zcond_iter 124 zt(i) = zt(i) + zcond_iter*RLVTT/RCPD 125 if (ABS(zcond_iter/alpha).lt.qthreshold) exit 136 zt(i) = zt(i) + zcond_iter*Lcp 137 if (ABS(zcond_iter/alpha/zqs(i)).lt.qthreshold) exit 138 ! if (ABS(zcond_iter/alpha).lt.qthreshold) exit 139 if (nn.eq.nitermax) print*,'itermax in largescale' 126 140 End do ! niter 127 141 zcond(i)=MAX(zcond(i),-(pq(i,k,igcm_h2o_ice)+pdq(i,k,igcm_h2o_ice)*ptimestep)) … … 129 143 else !standard case 130 144 131 zx_q(i) = (zq(i)+zdelq(i)+zqs(i))/2.0 !water vapor in cloudy sky145 zx_q(i) = (zq(i)+zdelq(i)+zqs(i))/2.0d0 !water vapor in cloudy sky 132 146 ! iterative process to stabilize the scheme when large water amounts JL12 133 zcond(i) = 0.0 147 zcond(i) = 0.0d0 134 148 Do nn=1,nitermax 135 call Lcpdqsat_water (zt(i),pplay(i,k),psat_tmp,zqs(i),zdqs(i),dlnpsat_tmp)136 zcond_iter = MAX(0.0 ,alpha*(zx_q(i)-zqs(i))/(1.+zdqs(i)))149 call Lcpdqsat_waterDP(zt(i),local_p,psat_tmp,zqs(i),zdqs(i),dlnpsat_tmp) 150 zcond_iter = MAX(0.0d0,alpha*(zx_q(i)-zqs(i))/(1.d0+zdqs(i))) 137 151 !zcond always postive! cannot evaporate clouds! 138 152 !this is why we must reevaporate before largescale 139 153 zx_q(i) = zx_q(i) - zcond_iter 140 154 zcond(i) = zcond(i) + zcond_iter 141 if (ABS(zcond_iter/alpha).lt.qthreshold) exit 142 zt(i) = zt(i) + zcond_iter*RLVTT/RCPD*rneb(i,k) 143 call Psat_water(zt(i),pplay(i,k),psat_tmp,zqs(i)) 155 if (ABS(zcond_iter/alpha/zqs(i)).lt.qthreshold) exit 156 ! if (ABS(zcond_iter/alpha).lt.qthreshold) exit 157 zt(i) = zt(i) + zcond_iter*Lcp*rneb(i,k) 158 call Psat_waterDP(zt(i),local_p,psat_tmp,zqs(i)) 159 if (nn.eq.nitermax) print*,'itermax in largescale' 144 160 End do ! niter 145 161 … … 153 169 pdqvaplsc(1:ngrid,k) = dqevap(1:ngrid,k) - zcond(1:ngrid) 154 170 pdqliqlsc(1:ngrid,k) = - pdqvaplsc(1:ngrid,k) 155 pdtlsc(1:ngrid,k) = pdqliqlsc(1:ngrid,k)* RLVTT/RCPD171 pdtlsc(1:ngrid,k) = pdqliqlsc(1:ngrid,k)*real(Lcp) 156 172 157 173 Enddo ! k= nlayermx, 1, -1 158 159 !print*,'qsat=',zqs 160 !print*,'q=',q 161 !print*,'dq=',pdqvaplsc*ptimestep 162 !print*,'dT in LS=',pdtlsc*ptimestep 163 164 !print*,'rice=',rice 165 !print*,'rneb=',rneb 174 166 175 167 176 return -
trunk/LMDZ.GENERIC/libf/phystd/moistadj.F90
r875 r1016 56 56 LOGICAL itest(ngrid) 57 57 REAL delta_q(ngrid, nlayermx) 58 REAL cp_new_t(nlayermx)58 DOUBLE PRECISION :: cp_new_t(nlayermx), v_cptt(ngrid,nlayermx) 59 59 REAL cp_delta_t(nlayermx) 60 REALv_cptj(nlayermx), v_cptjk1, v_ssig61 REAL v_ cptt(ngrid,nlayermx), v_p, v_t, v_zqs,v_cptt2,v_pratio,v_dlnpsat60 DOUBLE PRECISION :: v_cptj(nlayermx), v_cptjk1, v_ssig 61 REAL v_p, v_t, v_zqs,v_cptt2,v_pratio,v_dlnpsat 62 62 REAL zqs(ngrid,nlayermx), zdqs(ngrid,nlayermx),zpsat(ngrid,nlayermx),zdlnpsat(ngrid,nlayermx) 63 63 REAL zq1(ngrid), zq2(ngrid) 64 REALgamcpdz(ngrid,2:nlayermx)65 REALzdp, zdpm64 DOUBLE PRECISION :: gamcpdz(ngrid,2:nlayermx) 65 DOUBLE PRECISION :: zdp, zdpm 66 66 67 67 REAL zsat ! super-saturation 68 68 REAL zflo ! flotabilite 69 69 70 REALlocal_q(ngrid,nlayermx),local_t(ngrid,nlayermx)70 DOUBLE PRECISION :: local_q(ngrid,nlayermx),local_t(ngrid,nlayermx) 71 71 72 72 REAL zdelta, zcor, zcvm5 … … 139 139 v_zqs = (zqs(i,k-1)*zdpm + zqs(i,k)*zdp)/(zdpm+zdp) 140 140 v_cptt2 = (v_cptt(i,k-1)*zdpm + v_cptt(i,k)*zdp)/(zdpm+zdp) 141 v_pratio = ((1. /(1.+zpsat(i,k-1)/pplay(i,k-1)))*zdpm + (1./(1.+zpsat(i,k)/pplay(i,k)))*zdp)/(zdpm+zdp)141 v_pratio = ((1.-zpsat(i,k-1)/pplay(i,k-1))*zdpm + (1.-zpsat(i,k)/pplay(i,k))*zdp)/(zdpm+zdp) 142 142 v_dlnpsat = (zdlnpsat(i,k-1)*zdpm + zdlnpsat(i,k)*zdp)/(zdpm+zdp) 143 gamcpdz(i,k) = v_pratio*( (R/RCPD*v_cptt2*(1.- v_zqs) + RLVTT*v_zqs) * (pplay(i,k-1)-pplay(i,k))/pplev(i,k) ) &144 / (( 1.- v_zqs) + v_zqs * RCPV/RCPD + v_zqs * v_pratio * v_dlnpsat)143 gamcpdz(i,k) = ( (R/RCPD*v_cptt2*(1.- v_zqs) + RLVTT*v_zqs) * (pplay(i,k-1)-pplay(i,k))/pplev(i,k) ) & 144 / (((1.- v_zqs) + v_zqs * RCPV/RCPD)*v_pratio + v_zqs * v_dlnpsat) 145 145 ENDDO 146 146 ENDDO … … 210 210 call Lcpdqsat_water(v_t,v_p,zpsat(i,k),zqs(i,k),zdqs(i,k),zdlnpsat(i,k)) 211 211 212 213 214 ! print*,'i,k,zqs,cpdT=',i,k,zqs(i,k),cp_delta_t(k)215 212 ENDDO 216 213 Enddo ! nn=1,niter … … 240 237 v_zqs = (zqs(i,k-1)*zdpm + zqs(i,k)*zdp)/(zdpm+zdp) 241 238 v_cptt2 = (v_cptt(i,k-1)*zdpm + v_cptt(i,k)*zdp)/(zdpm+zdp) 242 v_pratio = ((1. /(1.+zpsat(i,k-1)/pplay(i,k-1)))*zdpm + (1./(1.+zpsat(i,k)/pplay(i,k)))*zdp)/(zdpm+zdp)239 v_pratio = ((1.-zpsat(i,k-1)/pplay(i,k-1))*zdpm + (1.-zpsat(i,k)/pplay(i,k))*zdp)/(zdpm+zdp) 243 240 v_dlnpsat = (zdlnpsat(i,k-1)*zdpm + zdlnpsat(i,k)*zdp)/(zdpm+zdp) 244 gamcpdz(i,k) = v_pratio*( (R/RCPD*v_cptt2*(1.- v_zqs) + RLVTT*v_zqs) * (pplay(i,k-1)-pplay(i,k))/pplev(i,k) ) &245 / (( 1.- v_zqs) + v_zqs * RCPV/RCPD + v_zqs * v_pratio * v_dlnpsat)241 gamcpdz(i,k) = ( (R/RCPD*v_cptt2*(1.- v_zqs) + RLVTT*v_zqs) * (pplay(i,k-1)-pplay(i,k))/pplev(i,k) ) & 242 / (((1.- v_zqs) + v_zqs * RCPV/RCPD)*v_pratio + v_zqs * v_dlnpsat) 246 243 ENDDO 247 244 … … 271 268 9999 CONTINUE ! loop over all the points 272 269 273 ! print*,'k1=',k1274 ! print*,'k2=',k2275 276 ! print*,'local_t=',local_t277 ! print*,'v_cptt=',v_cptt278 ! print*,'gamcpdz=',gamcpdz279 280 270 !----------------------------------------------------------------------- 281 271 ! Determine the cloud fraction (hypothese: la nebulosite a lieu … … 317 307 ENDDO 318 308 ENDDO 319 320 ! print*,'local_q BEFORE=',local_q321 309 322 310 DO k = 1, nlayermx -
trunk/LMDZ.GENERIC/libf/phystd/optci.F90
r961 r1016 121 121 if(kastprof)then 122 122 dz(k) = dpr(k)*(1000.0d0*8.3145d0/muvar(k))*TMID(K)/(g*PMID(K)) 123 U(k) = (Cmk*mugaz/(muvar(k)))*DPR(k)123 U(k) = Cmk*DPR(k)*mugaz/muvar(k) 124 124 else 125 dz(k) = dpr(k)*R*TMID(K)/(g*PMID(K)) 126 U(k) = Cmk*DPR(k) ! only Cmk line in optci.F 125 dz(k) = dpr(k)*R*TMID(K)/(g*PMID(K))*mugaz/muvar(k) 126 U(k) = Cmk*DPR(k)*mugaz/muvar(k) ! only Cmk line in optci.F 127 !JL13 the mugaz/muvar factor takes into account water meanmolecular weight if water is present 127 128 endif 128 129 -
trunk/LMDZ.GENERIC/libf/phystd/optcv.F90
r919 r1016 108 108 ! if we have continuum opacities, we need dz 109 109 if(kastprof)then 110 dz(k) = dpr(k)*(1000.0 *8.3145/muvar(k))*TMID(K)/(g*PMID(K))111 U(k) = (Cmk*mugaz/(muvar(k)))*DPR(k)110 dz(k) = dpr(k)*(1000.0d0*8.3145d0/muvar(k))*TMID(K)/(g*PMID(K)) 111 U(k) = Cmk*DPR(k)*mugaz/muvar(k) 112 112 else 113 dz(k) = dpr(k)*R*TMID(K)/(g*PMID(K)) 114 U(k) = Cmk*DPR(k) ! only Cmk line in optci.F 113 dz(k) = dpr(k)*R*TMID(K)/(g*PMID(K))*mugaz/muvar(k) 114 U(k) = Cmk*DPR(k)*mugaz/muvar(k) ! only Cmk line in optci.F 115 !JL13 the mugaz/muvar factor takes into account water meanmolecular weight if water is present 115 116 endif 116 117 -
trunk/LMDZ.GENERIC/libf/phystd/physiq.F90
r996 r1016 9 9 10 10 use radinc_h, only : L_NSPECTI,L_NSPECTV,naerkind 11 use watercommon_h, only : RLVTT, Psat_water 11 use watercommon_h, only : RLVTT, Psat_water,epsi 12 12 use gases_h 13 13 use radcommon_h, only: sigma … … 329 329 real dEzRadsw(ngrid,nlayermx),dEzRadlw(ngrid,nlayermx),dEzdiff(ngrid,nlayermx) 330 330 real dEdiffs(ngrid),dEdiff(ngrid) 331 real madjdE(ngrid), lscaledE(ngrid) 331 real madjdE(ngrid), lscaledE(ngrid),madjdEz(ngrid,nlayermx), lscaledEz(ngrid,nlayermx) 332 332 !JL12 conservation test for mean flow kinetic energy has been disabled temporarily 333 333 … … 711 711 call abort 712 712 endif 713 muvar(:,:)=0.0 ! only used for climate evolution studies (kcm1d) for now 713 if(water) then 714 muvar(1:ngrid,1:nlayermx)=mugaz/(1.e0+(1.e0/epsi-1.e0)*pq(1:ngrid,1:nlayermx,igcm_h2o_vap)) 715 muvar(1:ngrid,nlayermx+1)=mugaz/(1.e0+(1.e0/epsi-1.e0)*pq(1:ngrid,nlayermx,igcm_h2o_vap)) 716 ! take into account water effect on mean molecular weight 717 else 718 muvar(1:ngrid,1:nlayermx+1)=mugaz 719 endif 714 720 715 721 ! standard callcorrk … … 1059 1065 dtmoist(1:ngrid,1:nlayermx)=0. 1060 1066 1061 1067 call moistadj(ngrid,nq,pt,pq,pdq,pplev,pplay,dtmoist,dqmoist,ptimestep,rneb_man) 1062 1068 1063 1069 pdq(1:ngrid,1:nlayermx,igcm_h2o_vap) = pdq(1:ngrid,1:nlayermx,igcm_h2o_vap) & … … 1071 1077 if(enertest)then 1072 1078 dEtot=cpp*SUM(massarea(:,:)*dtmoist(:,:))/totarea 1079 madjdEz(:,:)=cpp*mass(:,:)*dtmoist(:,:) 1073 1080 do ig=1,ngrid 1074 1081 madjdE(ig) = cpp*SUM(mass(:,:)*dtmoist(:,:)) … … 1089 1096 ! Large scale condensation/evaporation 1090 1097 ! -------------------------------- 1091 1092 1098 call largescale(ngrid,nq,ptimestep,pplev,pplay,pt,pq,pdt,pdq,dtlscale,dqvaplscale,dqcldlscale,rneb_lsc) 1093 1099 … … 1099 1105 ! test energy conservation 1100 1106 if(enertest)then 1107 lscaledEz(:,:) = cpp*mass(:,:)*dtlscale(:,:) 1101 1108 do ig=1,ngrid 1102 1109 lscaledE(ig) = cpp*SUM(mass(:,:)*dtlscale(:,:)) … … 1534 1541 do l = 1, nlayer 1535 1542 do ig=1,ngrid 1536 ! call watersat(pt(ig,l),pplay(ig,l),qsat(ig,l))1537 1543 call Psat_water(zt(ig,l),pplay(ig,l),psat_tmp,qsat(ig,l)) 1538 1539 1544 RH(ig,l) = zq(ig,l,igcm_h2o_vap) / qsat(ig,l) 1540 1545 enddo … … 1698 1703 call writediagfi(ngrid,"ASR","absorbed stellar rad.","W m-2",2,fluxabs_sw) 1699 1704 call writediagfi(ngrid,"OLR","outgoing longwave rad.","W m-2",2,fluxtop_lw) 1705 ! call writediagfi(ngrid,"ASRcs","absorbed stellar rad (cs).","W m-2",2,fluxabs_sw1) 1706 ! call writediagfi(ngrid,"OLRcs","outgoing longwave rad (cs).","W m-2",2,fluxtop_lw1) 1707 ! call writediagfi(ngrid,"fluxsurfsw","sw surface flux.","W m-2",2,fluxsurf_sw) 1708 ! call writediagfi(ngrid,"fluxsurflw","lw back radiation.","W m-2",2,fluxsurf_lw) 1709 ! call writediagfi(ngrid,"fluxsurfswcs","sw surface flux (cs).","W m-2",2,fluxsurf_sw1) 1710 ! call writediagfi(ngrid,"fluxsurflwcs","lw back radiation (cs).","W m-2",2,fluxsurf_lw1) 1700 1711 call writediagfi(ngrid,"GND","heat flux from ground","W m-2",2,fluxgrd) 1701 1712 call writediagfi(ngrid,"DYN","dynamical heat input","W m-2",2,fluxdyn) … … 1715 1726 endif 1716 1727 if(watercond) then 1728 ! call writediagfi(ngrid,"lscaledEz","heat from largescale","W m-2",3,lscaledEz) 1729 ! call writediagfi(ngrid,"madjdEz","heat from moistadj","W m-2",3,madjdEz) 1717 1730 call writediagfi(ngrid,"lscaledE","heat from largescale","W m-2",2,lscaledE) 1718 1731 call writediagfi(ngrid,"madjdE","heat from moistadj","W m-2",2,madjdE) 1719 call writediagfi(ngrid,"RH","relative humidity"," ",3,RH) 1720 ! call writediagfi(ngrid,"qsatatm","atm qsat"," ",3,qsat) 1732 call writediagfi(ngrid,"qsatatm","atm qsat"," ",3,qsat) 1721 1733 ! call writediagfi(ngrid,"h2o_max_col","maximum H2O column amount","kg.m^-2",2,H2Omaxcol) 1722 1734 endif … … 1759 1771 call writediagfi(ngrid,"CLF","H2O cloud fraction"," ",3,cloudfrac) 1760 1772 call writediagfi(ngrid,"CLFt","H2O column cloud fraction"," ",2,totcloudfrac) 1773 call writediagfi(ngrid,"RH","relative humidity"," ",3,RH) 1761 1774 endif 1762 1775 -
trunk/LMDZ.GENERIC/libf/phystd/radcommon_h.F90
r959 r1016 36 36 ! TAURAY - Array (NSPECTV elements) of the pressure-independent 37 37 ! part of Rayleigh scattering optical depth. 38 ! TAURAYVAR - Array (NSPECTV elements) of the pressure-independent 39 ! part of Rayleigh scattering optical depth for the variable gas. 38 40 ! FZEROI - Fraction of zeros in the IR CO2 k-coefficients, for 39 41 ! each temperature, pressure, and spectral interval … … 61 63 REAL*8 BWNI(L_NSPECTI+1), WNOI(L_NSPECTI), DWNI(L_NSPECTI), WAVEI(L_NSPECTI) 62 64 REAL*8 BWNV(L_NSPECTV+1), WNOV(L_NSPECTV), DWNV(L_NSPECTV), WAVEV(L_NSPECTV) 63 REAL*8 STELLARF(L_NSPECTV), TAURAY(L_NSPECTV) 65 REAL*8 STELLARF(L_NSPECTV), TAURAY(L_NSPECTV), TAURAYVAR(L_NSPECTV) 64 66 65 67 REAL*8 blami(L_NSPECTI+1) -
trunk/LMDZ.GENERIC/libf/phystd/rain.F90
r863 r1016 66 66 INTEGER,PARAMETER :: ninter=5 67 67 68 logical, parameter :: evap_prec=.true.! Does the rain evaporate?68 logical,save :: evap_prec ! Does the rain evaporate? 69 69 70 70 ! for simple scheme … … 142 142 143 143 endif 144 145 write(*,*) "re-evaporate precipitations?" 146 evap_prec=.true. ! default value 147 call getin("evap_prec",evap_prec) 148 write(*,*) " evap_prec = ",evap_prec 144 149 145 150 firstcall = .false. -
trunk/LMDZ.GENERIC/libf/phystd/turbdiff.F90
r952 r1016 543 543 zcq0(ig) = qsat(ig)-dqsat(ig)*ptsrf(ig) 544 544 545 z1(ig) = pcapcal(ig)*ptsrf(ig) +cpp*zflux q(ig,1)*zct(ig,1)*zovExner(ig,1) &545 z1(ig) = pcapcal(ig)*ptsrf(ig) +cpp*zfluxt(ig,1)*zct(ig,1)*zovExner(ig,1) & 546 546 + zdplanck(ig)*ptsrf(ig) + pfluxsrf(ig)*ptimestep & 547 547 + zfluxq(ig,1)*dryness(ig)*RLVTT*((zdq(ig,1)-1.0)*zcq0(ig)+zcq(ig,1)) 548 548 549 z2(ig) = pcapcal(ig) + cpp*zflux q(ig,1)*(1.-zovExner(ig,1)*zdt(ig,1)) &549 z2(ig) = pcapcal(ig) + cpp*zfluxt(ig,1)*(1.-zovExner(ig,1)*zdt(ig,1)) & 550 550 + zdplanck(ig)+zfluxq(ig,1)*dryness(ig)*RLVTT*zdq0(ig)*(1.0-zdq(ig,1)) 551 551 -
trunk/LMDZ.GENERIC/libf/phystd/watercommon_h.F90
r989 r1016 175 175 end subroutine Tsat_water 176 176 177 !================================================================== 178 subroutine Psat_waterDP(T,p,psat,qsat) 179 180 implicit none 181 182 !================================================================== 183 ! Purpose 184 ! ------- 185 ! Compute the saturation vapor pressure and mass mixing ratio at saturation (kg/kg) 186 ! for a given pressure (Pa) and temperature (K) 187 ! DOUBLE PRECISION 188 ! Based on the Tetens formula from L.Li physical parametrization manual 189 ! 190 ! Authors 191 ! ------- 192 ! Jeremy Leconte (2012) 193 ! 194 !================================================================== 195 196 ! input 197 double precision, intent(in) :: T, p 198 199 ! output 200 double precision psat,qsat 201 202 ! JL12 variables for tetens formula 203 double precision,parameter :: Pref_solid_liquid=611.14d0 204 double precision,parameter :: Trefvaporization=35.86D0 205 double precision,parameter :: Trefsublimation=7.66d0 206 double precision,parameter :: Tmin=8.d0 207 double precision,parameter :: r3vaporization=17.269d0 208 double precision,parameter :: r3sublimation=21.875d0 209 210 ! checked vs. old watersat data 14/05/2012 by JL. 211 212 if (T.gt.T_h2O_ice_liq) then 213 psat = Pref_solid_liquid*Exp(r3vaporization*(T-T_h2O_ice_liq)/(T-Trefvaporization)) ! liquid / vapour 214 else if (T.lt.Tmin) then 215 print*, "careful, T<Tmin in psat water" 216 psat = Pref_solid_liquid*Exp(r3sublimation*(Tmin-T_h2O_ice_liq)/(Tmin-Trefsublimation)) ! min psat 217 else 218 psat = Pref_solid_liquid*Exp(r3sublimation*(T-T_h2O_ice_liq)/(T-Trefsublimation)) ! solid / vapour 219 endif 220 if(psat.gt.p) then 221 qsat=1.d0 222 else 223 qsat=epsi*psat/(p-(1.d0-epsi)*psat) 224 endif 225 return 226 end subroutine Psat_waterDP 227 228 229 230 231 !================================================================== 232 subroutine Lcpdqsat_waterDP(T,p,psat,qsat,dqsat,dlnpsat) 233 234 implicit none 235 236 !================================================================== 237 ! Purpose 238 ! ------- 239 ! Compute dqsat=L/cp*d (q_sat)/d T and dlnpsat=L/cp d(ln Psat)/d T 240 ! for a given temperature (K)! 241 ! Based on the Tetens formula from L.Li physical parametrization manual 242 ! 243 ! Authors 244 ! ------- 245 ! Jeremy Leconte (2012) 246 ! 247 !================================================================== 248 249 ! input 250 double precision T, p, psat, qsat 251 252 ! output 253 double precision dqsat,dlnpsat 254 255 ! JL12 variables for tetens formula 256 double precision,parameter :: Pref_solid_liquid=611.14d0 257 double precision,parameter :: Trefvaporization=35.86d0 258 double precision,parameter :: Tmin=8.d0 259 double precision,parameter :: Trefsublimation=7.66d0 260 double precision,parameter :: r3vaporization=17.269d0 261 double precision,parameter :: r3sublimation=21.875d0 262 263 double precision :: dummy 264 265 if (psat.gt.p) then 266 dqsat=0.d0 267 return 268 endif 269 270 if (T.gt.T_h2O_ice_liq) then 271 dummy = r3vaporization*(T_h2O_ice_liq-Trefvaporization)/(T-Trefvaporization)**2 ! liquid / vapour 272 else if (T.lt.Tmin) then 273 print*, "careful, T<Tmin in Lcp psat water" 274 dummy = r3sublimation*(T_h2O_ice_liq-Trefsublimation)/(Tmin-Trefsublimation)**2 ! solid / vapour 275 else 276 dummy = r3sublimation*(T_h2O_ice_liq-Trefsublimation)/(T-Trefsublimation)**2 ! solid / vapour 277 endif 278 279 dqsat=RLVTT/RCPD*qsat*(p/(p-(1.d0-epsi)*psat))*dummy 280 dlnpsat=RLVTT/RCPD*dummy 281 return 282 end subroutine Lcpdqsat_waterDP 283 177 284 178 285 end module watercommon_h
Note: See TracChangeset
for help on using the changeset viewer.
