Changeset 773
- Timestamp:
- Sep 5, 2012, 8:32:32 PM (12 years ago)
- Location:
- trunk/LMDZ.GENERIC
- Files:
-
- 8 edited
-
README (modified) (1 diff)
-
libf/phystd/ave_stelspec.F90 (modified) (4 diffs)
-
libf/phystd/hydrol.F90 (modified) (1 diff)
-
libf/phystd/largescale.F90 (modified) (1 diff)
-
libf/phystd/moistadj.F90 (modified) (1 diff)
-
libf/phystd/physiq.F90 (modified) (4 diffs)
-
libf/phystd/rcm1d.F (modified) (1 diff)
-
libf/phystd/sugas_corrk.F90 (modified) (7 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/LMDZ.GENERIC/README
r751 r773 750 750 - Bug fix in physiq : the size of OLR_nu is L_NSPECTI and not L_NSPECTV 751 751 - A more robust aerosol_mod + iniaerosol : problems with ifort+parallel solved, while still OK with other compilers and seq runs. 752 753 == 05/09/2012 == JL 754 - Correction of the calculation of the solar longitude in tlocked case. 755 -Can now handle any prograde resonance with nres=omega_rot/omega_orb. 756 -Sun now goes westward for the standard 2:1 case, as expected. 757 - In the gray case, the separation between kappa_IR and VI is now set by 758 wave number, independently of the usual IR/VISIBLE calculation separation. 759 i.e. kappa_IR can be used in the calculation of the downward stellar flux 760 if the wavenumber in the band is low enough and vice versa. 761 - In ave_stelspec, stellar flux averaging has been generalized to incorporate 762 very red/blue stellar spectra (great care must however be taken of the band 763 limit used for the corralated k distributions). 764 -Brown dwarf spectra from Allard et al. have been added. 765 -Any Black body temperature can now be used. 766 767 -
trunk/LMDZ.GENERIC/libf/phystd/ave_stelspec.F90
r374 r773 11 11 ! ------- 12 12 ! Robin Wordsworth (2010). 13 ! Generalized to very late spectral types (and Brown dwarfs) Jeremy Leconte (2012) 13 14 ! 14 15 ! Called by … … 34 35 35 36 integer Nfine 36 parameter(Nfine=5000)37 integer,parameter :: Nfineband=200 37 38 integer ifine 38 39 39 real lam(Nfine) 40 real stel_f(Nfine) 41 real band 40 real,allocatable :: lam(:),stel_f(:) 41 real band,lamm,lamp 42 42 real dl 43 43 44 character(len= 50) :: file_id45 character(len= 100) :: file_path44 character(len=100) :: file_id,file_id_lam 45 character(len=200) :: file_path,file_path_lam 46 46 47 47 real lam_temp … … 52 52 STELLAR(:)=0.0 53 53 54 ! load high resolution wavenumber data 55 file_id='/stellar_spectra/lam.txt' 56 file_path=TRIM(datadir)//TRIM(file_id) 57 58 open(110,file=file_path,form='formatted',status='old',iostat=ios) 59 if (ios.ne.0) then ! file not found 60 write(*,*) 'Error from ave_stelspec' 61 write(*,*) 'Data file ',trim(file_id),' not found.' 62 write(*,*)'Check that your path to datagcm:',trim(datadir) 63 write(*,*)' is correct. You can change it in callphys.def with:' 64 write(*,*)' datadir = /absolute/path/to/datagcm' 65 write(*,*)' Also check that there is a ',trim(file_id),' there.' 66 call abort 67 else 68 do ifine=1,Nfine 69 read(110,*) lam(ifine) 70 enddo 71 close(110) 72 endif 73 74 dl=lam(2)-lam(1) 75 76 54 print*,'enter ave_stellspec' 77 55 if(stelbbody)then 78 56 tstellar=stelTbb 79 do ifine=1,Nfine 80 lam_temp=lam(ifine) 81 call blackl(dble(lam_temp*1e-6),dble(tstellar),stel_temp) 82 stel_f(ifine)=stel_temp 83 enddo 57 Nfine=L_NSPECTV*Nfineband 58 do band=1,L_NSPECTV 59 lamm=10000.0/BWNV(band+1) 60 lamp=10000.0/BWNV(band) 61 dl=(lamp-lamm)/(Nfineband) 62 do ifine=1,Nfineband 63 lam_temp=lamm+(lamp-lamm)*(ifine-1.)/(Nfineband) 64 call blackl(dble(lam_temp*1e-6),dble(tstellar),stel_temp) 65 STELLAR(band)=STELLAR(band)+stel_temp*dl 66 enddo 67 end do 68 STELLAR(1:L_NSPECTV)=STELLAR(1:L_NSPECTV)/sum(STELLAR(1:L_NSPECTV)) 84 69 else 85 70 ! load high resolution stellar data 86 if(startype.eq.1)then 71 Select Case(startype) 72 Case(1) 87 73 file_id='/stellar_spectra/sol.txt' 88 74 tstellar=5800. 89 elseif(startype.eq.2)then 75 file_id_lam='/stellar_spectra/lam.txt' 76 Nfine=5000 77 Case(2) 90 78 file_id='/stellar_spectra/gl581.txt' 91 79 tstellar=3200. 92 elseif(startype.eq.3)then 80 file_id_lam='/stellar_spectra/lam.txt' 81 Nfine=5000 82 Case(3) 93 83 file_id='/stellar_spectra/adleo.txt' 94 84 tstellar=3200. 95 elseif(startype.eq.4)then 85 file_id_lam='/stellar_spectra/lam.txt' 86 Nfine=5000 87 Case(4) 96 88 file_id='/stellar_spectra/gj644.txt' 97 89 print*,'Find out tstellar before using this star!' 98 90 call abort 99 elseif(startype.eq.5)then 91 file_id_lam='/stellar_spectra/lam.txt' 92 Nfine=5000 93 Case(5) 100 94 file_id='/stellar_spectra/hd128167.txt' 101 95 tstellar=6700. ! Segura et al. (2003) 102 else 96 file_id_lam='/stellar_spectra/lam.txt' 97 Nfine=5000 98 Case(6) 99 file_id='/stellar_spectra/BD_Teff-1600K.txt' 100 tstellar=1600. ! Segura et al. (2003) 101 file_id_lam='/stellar_spectra/lamBD.txt' 102 Nfine=5000 103 Case(7) 104 file_id='/stellar_spectra/BD_Teff-1000K.txt' 105 tstellar=1000. ! Segura et al. (2003) 106 file_id_lam='/stellar_spectra/lamBD.txt' 107 Nfine=5000 108 Case Default 103 109 print*,'Error: unknown star type chosen' 104 110 call abort 111 End Select 112 113 allocate(lam(Nfine),stel_f(Nfine)) 114 115 file_path_lam=TRIM(datadir)//TRIM(file_id_lam) 116 open(110,file=file_path_lam,form='formatted',status='old',iostat=ios) 117 if (ios.ne.0) then ! file not found 118 write(*,*) 'Error from ave_stelspec' 119 write(*,*) 'Data file ',trim(file_id_lam),' not found.' 120 write(*,*)'Check that your path to datagcm:',trim(datadir) 121 write(*,*)' is correct. You can change it in callphys.def with:' 122 write(*,*)' datadir = /absolute/path/to/datagcm' 123 write(*,*)' Also check that there is a ',trim(file_id_lam),' there.' 124 call abort 125 else 126 do ifine=1,Nfine 127 read(110,*) lam(ifine) 128 enddo 129 close(110) 105 130 endif 106 131 132 133 ! load high resolution wavenumber data 107 134 file_path=TRIM(datadir)//TRIM(file_id) 108 135 open(111,file=file_path,form='formatted',status='old',iostat=ios) … … 121 148 close(111) 122 149 endif 150 151 ! sum data by band 152 band=1 153 Do while(lam(1).lt. real(10000.0/BWNV(band+1))) 154 if (band.gt.L_NSPECTV) exit 155 band=band+1 156 enddo 157 dl=lam(2)-lam(1) 158 STELLAR(band)=STELLAR(band)+stel_f(1)*dl 159 do ifine = 2,Nfine 160 if(lam(ifine) .gt. real(10000.0/BWNV(band)))then 161 band=band-1 162 endif 163 if(band .lt. 1) exit 164 dl=lam(ifine)-lam(ifine-1) 165 STELLAR(band)=STELLAR(band)+stel_f(ifine)*dl 166 end do 167 168 169 STELLAR(1:L_NSPECTV)=STELLAR(1:L_NSPECTV)/sum(STELLAR(1:L_NSPECTV)) 170 deallocate(lam,stel_f) 171 123 172 endif 124 173 125 126 ! sum data by band127 band=1128 do ifine = 1,Nfine129 130 if(lam(Nfine-ifine+1) .lt. real(10000.0/BWNV(band+1)))then131 band=band+1132 endif133 if(band .gt. L_NSPECTV)then134 goto 9999 ! ok, ok, I know they're evil135 endif136 STELLAR(band)=STELLAR(band)+stel_f(Nfine-ifine+1)*dl137 138 end do139 140 9999 continue141 STELLAR=STELLAR/sum(STELLAR)142 143 144 174 end subroutine ave_stelspec -
trunk/LMDZ.GENERIC/libf/phystd/hydrol.F90
r650 r773 194 194 hice(ig) = max(hice(ig),0.0) 195 195 hice(ig) = min(hice(ig),maxicethick) 196 pdtsurf_hyd(ig) = (hice(ig) - hicebis(ig))*RLFTT* & 197 rhowater/pcapcal(ig)/ptimestep 196 pdtsurf_hyd(ig) = (hice(ig) - hicebis(ig))*RLFTT*rhowater/pcapcal(ig)/ptimestep 198 197 albedo(ig) = albedoice 199 198 -
trunk/LMDZ.GENERIC/libf/phystd/largescale.F90
r728 r773 91 91 DO i = 1, ngridmx 92 92 93 if(zt(i).le.15.) zt(i)=15. ! check too low temperatures 93 if(zt(i).le.15.) then 94 print*,'in lsc',i,zt(i) 95 zt(i)=15. ! check too low temperatures 96 endif 94 97 ! call watersat(zt(i),pplay(i,k),zqs(i)) 95 98 call Psat_water(zt(i),pplay(i,k),psat_tmp,zqs(i)) -
trunk/LMDZ.GENERIC/libf/phystd/moistadj.F90
r728 r773 50 50 ! Local variables 51 51 INTEGER i, k, iq, nn 52 INTEGER, PARAMETER :: niter = 652 INTEGER, PARAMETER :: niter = 1 53 53 INTEGER k1, k1p, k2, k2p 54 54 LOGICAL itest(ngridmx) -
trunk/LMDZ.GENERIC/libf/phystd/physiq.F90
r751 r773 642 642 if (tlocked) then 643 643 ztim1=SIN(declin) 644 ztim2=COS(declin)*COS(2.*pi*(zday/year_day) - zls*nres) 645 ztim3=-COS(declin)*SIN(2.*pi*(zday/year_day) - zls*nres) 644 ! ztim2=COS(declin)*COS(2.*pi*(zday/year_day) - zls*nres) 645 ! ztim3=-COS(declin)*SIN(2.*pi*(zday/year_day) - zls*nres) 646 ! JL12 corrects tidally resonant cases. nres=omega_rot/omega_orb 647 ztim2=COS(declin)*COS(2.*pi*(zday/year_day)*(nres-1.)) 648 ztim3=-COS(declin)*SIN(2.*pi*(zday/year_day)*(nres-1.)) 646 649 647 650 call stelang(ngrid,sinlon,coslon,sinlat,coslat, & … … 1017 1020 ! Moist convection 1018 1021 ! ---------------- 1019 ! Re-evaporate cloud water/ice 1020 ! dqevap1(1:ngridmx,1:nlayermx)=0. 1021 ! dtevap1(1:ngridmx,1:nlayermx)=0. 1022 ! call evap(ptimestep,pt,pq,pdq,pdt,dqevap1,dtevap1,qevap,tevap) 1023 1022 1024 1023 dqmoist(1:ngridmx,1:nlayermx,1:nqmx)=0. 1025 1024 dtmoist(1:ngridmx,1:nlayermx)=0. … … 1032 1031 +dqmoist(1:ngridmx,1:nlayermx,igcm_h2o_ice) 1033 1032 pdt(1:ngridmx,1:nlayermx) = pdt(1:ngridmx,1:nlayermx)+dtmoist(1:ngridmx,1:nlayermx) 1034 1035 ! do l=1,nlayermx1036 ! do ig=1,ngridmx1037 ! if(isnan(dqmoist(ig,l,igcm_h2o_ice))) print*,'Nan in dqmoist ice, abort,ig,l',ig,l1038 ! if(isnan(dqmoist(ig,l,igcm_h2o_vap))) print*,'Nan in dqmoist vap, abort,ig,l',ig,l1039 ! if(isnan(dtmoist(ig,l))) print*,'Nan in dtmoist, abort,ig,l',ig,l1040 ! if(isnan(dqmoist(ig,l,igcm_h2o_ice)).or.isnan(dqmoist(ig,l,igcm_h2o_vap)).or.isnan(dtmoist(ig,l))) STOP1041 ! enddo1042 ! enddo1043 ! print*,'after moist',dqmoist(185,1,igcm_h2o_ice),dqmoist(185,1,igcm_h2o_vap),dtmoist(185,1)*86400.1044 1033 1045 1034 !------------------------- … … 1051 1040 enddo 1052 1041 print*,'In moistadj atmospheric energy change =',dEtot,' W m-2' 1053 print*,'In moistadj MAX atmospheric energy change =',MAXVAL(pdt(:,:))*86400.,'K/day' 1042 print*,'In moistadj MAX atmospheric energy change =',MAXVAL(dtmoist(:,:))*ptimestep,'K/step' 1043 print*,'In moistadj MIN atmospheric energy change =',MINVAL(dtmoist(:,:))*ptimestep,'K/step' 1054 1044 1055 1045 ! test energy conservation -
trunk/LMDZ.GENERIC/libf/phystd/rcm1d.F
r751 r773 299 299 !!! - physical constants: nevermind, things are done allright below 300 300 !!! - physical frequency: nevermind, in inifis this is a simple print 301 cpp=1.d-7 !JL because we divide by cpp in inifis, there may be a more elegant solution 301 302 CALL inifis(1,llm,day0,daysec,dtphys, 302 303 . lati,long,area,rad,g,r,cpp) -
trunk/LMDZ.GENERIC/libf/phystd/sugas_corrk.F90
r716 r773 24 24 use radcommon_h, only : tgasref,tgasmin,tgasmax 25 25 use radcommon_h, only : gasv,gasi,FZEROI,FZEROV,gweight 26 use radcommon_h, only : wrefvar 26 use radcommon_h, only : wrefvar,WNOI,WNOV 27 27 use datafile_mod, only: datadir 28 28 … … 49 49 50 50 real*8 x, xi(4), yi(4), ans, p 51 real kappa_IR, kappa_VI 51 ! For gray case (JL12) 52 real kappa_IR, kappa_VI, IR_VI_wnlimit, nVI_limit,nIR_limit 52 53 53 54 integer ngas, igas, jgas … … 274 275 ! Get gaseous k-coefficients and interpolate onto finer pressure grid 275 276 277 if (graybody) then 278 ! constant absorption coefficient in visible 279 write(*,*)"graybody: constant absorption coefficient in visible:" 280 kappa_VI=-100000. 281 call getin("kappa_VI",kappa_VI) 282 write(*,*)" kappa_VI = ",kappa_VI 283 kappa_VI=kappa_VI*1.e4* mugaz * 1.672621e-27 ! conversion from m^2/kg to cm^2/molecule 284 285 ! constant absorption coefficient in IR 286 write(*,*)"graybody: constant absorption coefficient in InfraRed:" 287 kappa_IR=-100000. 288 call getin("kappa_IR",kappa_IR) 289 write(*,*)" kappa_IR = ",kappa_IR 290 kappa_IR=kappa_IR*1.e4* mugaz * 1.672621e-27 ! conversion from m^2/kg to cm^2/molecule 291 292 ! wavelength used to separate IR from VI 293 IR_VI_wnlimit=3000. 294 write(*,*)"graybody: Visible / Infrared separation set at",10000./IR_VI_wnlimit,"um" 295 296 nVI_limit=L_NSPECTV 297 Do nw=1,L_NSPECTV 298 if (WNOV(nw).gt.IR_VI_wnlimit) then 299 nVI_limit=nw-1 300 exit 301 endif 302 End do 303 nIR_limit=L_NSPECTI 304 Do nw=1,L_NSPECTI 305 if (WNOI(nw).gt.IR_VI_wnlimit) then 306 nIR_limit=nw-1 307 exit 308 endif 309 End do 310 write(*,*)"graybody: Visible / Infrared separation set at band ",nIR_limit,nVI_limit 311 312 End if 313 314 276 315 ! VISIBLE 277 if (callgasvis.and.(.not.TRIM(corrkdir).eq.'null').and.(.not. graybody)) then316 if (callgasvis.and.(.not.TRIM(corrkdir).eq.'null').and.(.not.TRIM(corrkdir).eq.'null_LowTeffStar').and.(.not.graybody)) then 278 317 file_id='/corrk_data/'//trim(adjustl(banddir))//'/corrk_gcm_VI.dat' 279 318 file_path=TRIM(datadir)//TRIM(file_id) … … 293 332 294 333 else if (callgasvis.and.graybody) then 295 ! constant absorption coefficient in visible296 write(*,*)"constant absorption coefficient in visible:" 297 kappa_VI = -100000.298 call getin("kappa_VI",kappa_VI) 299 write(*,*)" kappa_VI = ",kappa_VI 300 kappa_VI = kappa_VI * 1.e4 * mugaz * 1.672621e-27 ! conversion from m^2/kg to cm^2/molecule301 gasv8(:,:,:,:,:) =kappa_VI302 334 if(nVI_limit.eq.0) then 335 gasv8(1:L_NTREF,1:L_NPREF,1:L_REFVAR,1:L_NSPECTV,1:L_NGAUSS)=kappa_VI 336 else if (nVI_limit.eq.L_NSPECTV) then 337 gasv8(1:L_NTREF,1:L_NPREF,1:L_REFVAR,1:L_NSPECTV,1:L_NGAUSS)=kappa_IR 338 else 339 gasv8(1:L_NTREF,1:L_NPREF,1:L_REFVAR,1:nVI_limit,1:L_NGAUSS)=kappa_IR 340 gasv8(1:L_NTREF,1:L_NPREF,1:L_REFVAR,nVI_limit+1:L_NSPECTV,1:L_NGAUSS)=kappa_VI 341 end if 303 342 else 304 343 print*,'Visible corrk gaseous absorption is set to zero.' 305 gasv8(:,:,:,:,:) = 0.0 306 344 gasv8(1:L_NTREF,1:L_NPREF,1:L_REFVAR,1:L_NSPECTV,1:L_NGAUSS)=0.0 307 345 endif 308 346 309 347 ! INFRA-RED 310 if ((.not.TRIM(corrkdir).eq.'null').and.(.not. graybody)) then348 if ((.not.TRIM(corrkdir).eq.'null').and.(.not.TRIM(corrkdir).eq.'null_LowTeffStar').and.(.not.graybody)) then 311 349 file_id='/corrk_data/'//trim(adjustl(banddir))//'/corrk_gcm_IR.dat' 312 350 file_path=TRIM(datadir)//TRIM(file_id) … … 364 402 365 403 else if (graybody) then 366 ! constant absorption coefficient in IR367 write(*,*)"constant absorption coefficient in InfraRed:" 368 kappa_IR = -100000.369 call getin("kappa_IR",kappa_IR) 370 write(*,*)" kappa_IR = ",kappa_IR 371 kappa_IR = kappa_IR * 1.e4 * mugaz * 1.672621e-27 ! conversion from m^2/kg to cm^2/molecule 372 gasi8(:,:,:,:,:) = kappa_IR 373 404 if(nIR_limit.eq.0) then 405 gasi8(1:L_NTREF,1:L_NPREF,1:L_REFVAR,1:L_NSPECTI,1:L_NGAUSS)=kappa_VI 406 else if (nIR_limit.eq.L_NSPECTI) then 407 gasi8(1:L_NTREF,1:L_NPREF,1:L_REFVAR,1:L_NSPECTI,1:L_NGAUSS)=kappa_IR 408 else 409 gasi8(1:L_NTREF,1:L_NPREF,1:L_REFVAR,1:nIR_limit,1:L_NGAUSS)=kappa_IR 410 gasi8(1:L_NTREF,1:L_NPREF,1:L_REFVAR,nIR_limit+1:L_NSPECTI,1:L_NGAUSS)=kappa_VI 411 end if 374 412 else 375 413 print*,'Infrared corrk gaseous absorption is set to zero.' 376 gasi8(:,:,:,:,:) = 0.0 377 414 gasi8(1:L_NTREF,1:L_NPREF,1:L_REFVAR,1:L_NSPECTV,1:L_NGAUSS)=0.0 378 415 endif 379 416 … … 389 426 do nw=1,L_NSPECTV 390 427 if(gasv8(nt,np,nh,nw,ng).gt.1.0d-200) then 391 gasv8(nt,np,nh,nw,ng) = & 392 log10(gasv8(nt,np,nh,nw,ng)) 428 gasv8(nt,np,nh,nw,ng) = log10(gasv8(nt,np,nh,nw,ng)) 393 429 else 394 430 gasv8(nt,np,nh,nw,ng) = -200.0 … … 398 434 do nw=1,L_NSPECTI 399 435 if(gasi8(nt,np,nh,nw,ng).gt.1.0d-200) then 400 gasi8(nt,np,nh,nw,ng) = & 401 log10(gasi8(nt,np,nh,nw,ng)) 436 gasi8(nt,np,nh,nw,ng) = log10(gasi8(nt,np,nh,nw,ng)) 402 437 else 403 438 gasi8(nt,np,nh,nw,ng) = -200.0
Note: See TracChangeset
for help on using the changeset viewer.
