Changeset 526 for trunk/LMDZ.GENERIC/libf/phystd/physiq.F90
- Timestamp:
- Feb 13, 2012, 6:30:26 PM (13 years ago)
- File:
-
- 1 edited
-
trunk/LMDZ.GENERIC/libf/phystd/physiq.F90 (modified) (22 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/LMDZ.GENERIC/libf/phystd/physiq.F90
r486 r526 7 7 pdu,pdv,pdt,pdq,pdpsrf,tracerdyn) 8 8 9 use radinc_h, only : naerkind 9 use radinc_h, only : naerkind,L_NSPECTI,L_NSPECTV 10 10 use watercommon_h, only : mx_eau_sol, To, RLVTT, mH2O 11 11 use gases_h … … 116 116 #include "comsaison.h" 117 117 #include "control.h" 118 #include "comg1d.h"119 118 #include "tracer.h" 120 119 #include "watercap.h" … … 190 189 191 190 character*80 fichier 192 integer l,ig,ierr,iq,i, tapphys 191 integer l,ig,ierr,iq,i, tapphys,nw 193 192 194 193 real fluxsurf_lw(ngridmx) ! incident LW (IR) surface flux (W.m-2) … … 199 198 real fluxtop_dn(ngridmx) 200 199 real fluxdyn(ngridmx) ! horizontal heat transport by dynamics 201 save fluxsurf_lw,fluxsurf_sw,fluxtop_lw,fluxabs_sw,fluxtop_dn,fluxdyn 200 real OLR_nu(ngridmx,L_NSPECTI)! Outgoing LW radition in each band (Normalized to the band width (W/m2/cm-1) 201 real OSR_nu(ngridmx,L_NSPECTV)! Outgoing SW radition in each band (Normalized to the band width (W/m2/cm-1) 202 save fluxsurf_lw,fluxsurf_sw,fluxtop_lw,fluxabs_sw,fluxtop_dn,fluxdyn,OLR_nu,OSR_nu 202 203 203 204 … … 326 327 real zpopskNI(ngridmx,nlayermx) 327 328 328 ! included by RW to make 1D saves not every timestep329 integer countG1D330 save countG1D331 329 332 330 ! included by BC for evaporation … … 357 355 ! included by BC for double radiative transfer call 358 356 logical clearsky 359 360 real zdtsw1(ngridmx,nlayermx),zdtsw2(ngridmx,nlayermx) 361 real zdtlw1(ngridmx,nlayermx),zdtlw2(ngridmx,nlayermx) 362 real fluxsurf_lw1(ngridmx), fluxsurf_lw2(ngridmx) 363 real fluxsurf_sw1(ngridmx), fluxsurf_sw2(ngridmx) 364 real fluxtop_lw1(ngridmx), fluxtop_lw2(ngridmx) 365 real fluxabs_sw1(ngridmx), fluxabs_sw2(ngridmx) 366 real tau_col1(ngrid), tau_col2(ngrid) 357 real zdtsw1(ngridmx,nlayermx) 358 real zdtlw1(ngridmx,nlayermx) 359 real fluxsurf_lw1(ngridmx) 360 real fluxsurf_sw1(ngridmx) 361 real fluxtop_lw1(ngridmx) 362 real fluxabs_sw1(ngridmx) 363 real tau_col1(ngrid) 364 real OLR_nu1(ngrid,L_NSPECTI) 365 real OSR_nu1(ngrid,L_NSPECTV) 367 366 real tf, ntf 368 367 … … 423 422 if (firstcall) then 424 423 425 if(ngrid.eq.1)then426 saveG1D=day_step*10427 !saveG1D=1428 !saveG1D=-1 ! to never save429 countG1D=1430 endif431 424 432 425 ! variables set to 0 … … 447 440 ! read startfi (initial parameters) 448 441 ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 449 442 print*,'in physiqu i am in phyetat0' 450 443 call phyetat0("startfi.nc",0,0,nsoilmx,nq, & 451 444 day_ini,time_phys,tsurf,tsoil,emis,q2,qsurf, & … … 681 674 682 675 if (callrad) then 683 if( mod(icount-1,iradia).eq.0 ) then676 if( mod(icount-1,iradia).eq.0.or.lastcall) then 684 677 685 678 ! Compute local stellar zenith angles … … 720 713 endif 721 714 muvar(:,:)=0.0 ! only used for climate evolution for now 722 723 ! Option to call scheme once for clear regions, once for cloudy regions (BC) 715 716 717 ! standard callcorrk 718 clearsky=.false. 719 call callcorrk(ngrid,nlayer,pq,nq,qsurf, & 720 albedo,emis,mu0,pplev,pplay,pt, & 721 tsurf,fract,dist_star,aerosol,cpp3D,muvar, & 722 zdtlw,zdtsw,fluxsurf_lw,fluxsurf_sw,fluxtop_lw, & 723 fluxabs_sw,fluxtop_dn,OLR_nu,OSR_nu, & 724 reffrad,tau_col,cloudfrac,totcloudfrac, & 725 clearsky,firstcall,lastcall) 726 727 728 ! Option to call scheme once more for clear regions 724 729 if(CLFvarying)then 725 730 726 !!! THIS IS NOT A VERY GOOD STRATEGY HERE (everything should be in callcorrk) 727 !!! ---> PROBLEMS WITH ALLOCATED ARRAYS 728 !!! (temporary solution in callcorrk: do not deallocate if CLFvarying ...) 729 731 !!! ---> PROBLEMS WITH ALLOCATED ARRAYS 732 !!! (temporary solution in callcorrk: do not deallocate if CLFvarying ...) 730 733 clearsky=.true. 731 734 call callcorrk(ngrid,nlayer,pq,nq,qsurf, & … … 733 736 tsurf,fract,dist_star,aerosol,cpp3D,muvar, & 734 737 zdtlw1,zdtsw1,fluxsurf_lw1,fluxsurf_sw1,fluxtop_lw1, & 735 fluxabs_sw1,fluxtop_dn,reffrad,tau_col1, & 736 cloudfrac,totcloudfrac,clearsky,firstcall,lastcall) 737 738 clearsky=.false. 739 call callcorrk(ngrid,nlayer,pq,nq,qsurf, & 740 albedo,emis,mu0,pplev,pplay,pt, & 741 tsurf,fract,dist_star,aerosol,cpp3D,muvar, & 742 zdtlw2,zdtsw2,fluxsurf_lw2,fluxsurf_sw2,fluxtop_lw2, & 743 fluxabs_sw2,fluxtop_dn,reffrad,tau_col2, & 744 cloudfrac,totcloudfrac,clearsky,firstcall,lastcall) 738 fluxabs_sw1,fluxtop_dn,OLR_nu1,OSR_nu1, & 739 reffrad,tau_col1,cloudfrac,totcloudfrac, & 740 clearsky,firstcall,lastcall) 741 clearsky = .false. !! just in case. 745 742 746 743 ! Sum the fluxes and heating rates from cloudy/clear cases … … 749 746 ntf=1.-tf 750 747 751 fluxsurf_lw(ig) = ntf*fluxsurf_lw1(ig) + tf*fluxsurf_lw 2(ig)752 fluxsurf_sw(ig) = ntf*fluxsurf_sw1(ig) + tf*fluxsurf_sw 2(ig)753 fluxtop_lw(ig) = ntf*fluxtop_lw1(ig) + tf*fluxtop_lw 2(ig)754 fluxabs_sw(ig) = ntf*fluxabs_sw1(ig) + tf*fluxabs_sw 2(ig)755 tau_col(ig) = ntf*tau_col1(ig) + tf*tau_col 2(ig)748 fluxsurf_lw(ig) = ntf*fluxsurf_lw1(ig) + tf*fluxsurf_lw(ig) 749 fluxsurf_sw(ig) = ntf*fluxsurf_sw1(ig) + tf*fluxsurf_sw(ig) 750 fluxtop_lw(ig) = ntf*fluxtop_lw1(ig) + tf*fluxtop_lw(ig) 751 fluxabs_sw(ig) = ntf*fluxabs_sw1(ig) + tf*fluxabs_sw(ig) 752 tau_col(ig) = ntf*tau_col1(ig) + tf*tau_col(ig) 756 753 757 754 do l=1,nlayer 758 zdtlw(ig,l) = ntf*zdtlw1(ig,l) + tf*zdtlw 2(ig,l)759 zdtsw(ig,l) = ntf*zdtsw1(ig,l) + tf*zdtsw 2(ig,l)755 zdtlw(ig,l) = ntf*zdtlw1(ig,l) + tf*zdtlw(ig,l) 756 zdtsw(ig,l) = ntf*zdtsw1(ig,l) + tf*zdtsw(ig,l) 760 757 enddo 758 759 do nw=1,L_NSPECTV 760 OSR_nu(ig,nw) = ntf*OSR_nu1(ig,nw) + tf*OSR_nu(ig,nw) 761 enddo 762 do nw=1,L_NSPECTI 763 OLR_nu(ig,nw) = ntf*OLR_nu1(ig,nw) + tf*OLR_nu(ig,nw) 764 enddo 765 761 766 enddo 762 767 763 ! Otherwise standard callcorrk 764 else 765 766 clearsky=.false. 767 call callcorrk(ngrid,nlayer,pq,nq,qsurf, & 768 albedo,emis,mu0,pplev,pplay,pt, & 769 tsurf,fract,dist_star,aerosol,cpp3D,muvar, & 770 zdtlw,zdtsw,fluxsurf_lw,fluxsurf_sw,fluxtop_lw, & 771 fluxabs_sw,fluxtop_dn,reffrad,tau_col, & 772 cloudfrac,totcloudfrac,clearsky,firstcall,lastcall) 773 774 endif 768 endif !CLFvarying 775 769 776 770 ! Radiative flux from the sky absorbed by the surface (W.m-2) … … 872 866 print*,'In corrk SW atmospheric energy change =',dEtotSW,' W m-2' 873 867 print*,'In corrk LW atmospheric energy change =',dEtotLW,' W m-2' 868 print*,'atmospheric energy change (SW+LW) =',dEtotLW+dEtotSW,' W m-2' 874 869 print*,'In corrk SW surface energy change =',dEtotsSW,' W m-2' 875 870 print*,'In corrk LW surface energy change =',dEtotsLW,' W m-2' 871 print*,'surface energy change (SW+LW) =',dEtotsLW+dEtotsSW,' W m-2' 876 872 endif 877 873 !------------------------- … … 884 880 885 881 if (calldifv) then 886 882 887 883 do ig=1,ngrid 888 884 zflubid(ig)=fluxrad(ig)+fluxgrd(ig) … … 1743 1739 OLR = OLR + area(ig)*fluxtop_lw(ig) 1744 1740 GND = GND + area(ig)*fluxgrd(ig) 1741 if(.not.callsoil) GND=GND+ area(ig)*fluxrad(ig) 1745 1742 DYN = DYN + area(ig)*fluxdyn(ig) 1746 1743 … … 1765 1762 if(enertest)then 1766 1763 print*,'SW energy balance SW++ - ASR = ',dEtotSW+dEtotsSW-ASR/totarea,' W m-2' 1764 print*,'LW energy balance LW++ + ***ASR*** = ',dEtotLW+dEtotsLW+ASR/totarea,' W m-2' 1767 1765 endif 1768 1766 1769 1767 if(meanOLR)then 1770 if((ngridmx.gt.1) .or. ( countG1D.eq.saveG1D))then1768 if((ngridmx.gt.1) .or. (mod(icount-1,nint(ecritphy)).eq.0))then 1771 1769 ! to record global radiative balance 1772 1770 open(92,file="rad_bal.out",form='formatted',access='append') … … 1857 1855 1858 1856 if(meanOLR)then 1859 if((ngridmx.gt.1) .or. ( countG1D.eq.saveG1D))then1857 if((ngridmx.gt.1) .or. (mod(icount-1,nint(ecritphy)).eq.0))then 1860 1858 ! to record global water balance 1861 1859 open(98,file="h2o_bal.out",form='formatted',access='append') … … 1891 1889 1892 1890 1893 if (ngrid.ne.1) then1894 1891 print*,'' 1895 1892 print*,'--> Ls =',zls*180./pi … … 1977 1974 ! "Solar radiative flux to space","W.m-2",2, & 1978 1975 ! fluxtop_sw_tot) 1976 1977 call wstats(ngrid,"ISR","incoming stellar rad.","W m-2",2,fluxtop_dn) 1978 call wstats(ngrid,"ASR","absorbed stellar rad.","W m-2",2,fluxabs_sw) 1979 call wstats(ngrid,"OLR","outgoing longwave rad.","W m-2",2,fluxtop_lw) 1980 1979 1981 call wstats(ngrid,"temp","Atmospheric temperature","K",3,zt) 1980 1982 call wstats(ngrid,"u","Zonal (East-West) wind","m.s-1",3,zu) … … 1984 1986 1985 1987 if (tracer) then 1988 do iq=1,nq 1989 call wstats(ngrid,noms(iq),noms(iq),'kg/kg',3,zq(1,1,iq)) 1990 call wstats(ngridmx,trim(noms(iq))//'_surf',trim(noms(iq))//'_surf', & 1991 'kg m^-2',2,qsurf(1,iq) ) 1992 1993 call wstats(ngridmx,trim(noms(iq))//'_col',trim(noms(iq))//'_col', & 1994 'kg m^-2',2,qcol(1,iq) ) 1995 call wstats(ngridmx,trim(noms(iq))//'_reff', & 1996 trim(noms(iq))//'_reff', & 1997 'm',3,reffrad(1,1,iq)) 1998 end do 1986 1999 if (water) then 1987 2000 vmr=zq(1:ngridmx,1:nlayermx,igcm_h2o_vap)*mugaz/mmol(igcm_h2o_vap) … … 2041 2054 2042 2055 ! Output aerosols 2043 call writediagfi(ngridmx,'CO2ice_reff','CO2ice_reff','m',3,reffrad(1,1,1))2044 call writediagfi(ngridmx,'H2Oice_reff','H2Oice_reff','m',3,reffrad(1,1,2))2045 call writediagfi(ngridmx,'CO2ice_reffcol','CO2ice_reffcol','um kg m^-2',2,reffcol(1,1))2046 call writediagfi(ngridmx,'H2Oice_reffcol','H2Oice_reffcol','um kg m^-2',2,reffcol(1,2))2056 if (igcm_co2_ice.ne.0) call writediagfi(ngridmx,'CO2ice_reff','CO2ice_reff','m',3,reffrad(1,1,1)) 2057 if (igcm_h2o_ice.ne.0) call writediagfi(ngridmx,'H2Oice_reff','H2Oice_reff','m',3,reffrad(1,1,2)) 2058 if (igcm_co2_ice.ne.0) call writediagfi(ngridmx,'CO2ice_reffcol','CO2ice_reffcol','um kg m^-2',2,reffcol(1,1)) 2059 if (igcm_h2o_ice.ne.0) call writediagfi(ngridmx,'H2Oice_reffcol','H2Oice_reffcol','um kg m^-2',2,reffcol(1,2)) 2047 2060 2048 2061 ! Output tracers … … 2093 2106 endif 2094 2107 2095 ! ---------------------------------------------------------- 2096 ! Output in netcdf file "diagsoil.nc" for subterranean 2097 ! variables (output every "ecritphy", as for writediagfi) 2098 ! ---------------------------------------------------------- 2099 2100 ! Write soil temperature 2101 !call writediagsoil(ngrid,"soiltemp","Soil temperature","K",3,tsoil) 2102 2103 else ! if(ngrid.eq.1) 2104 2105 ! ---------------------------------------------------------------------- 2106 ! Output in grads file "g1d" (ONLY when using testphys1d) 2107 ! ---------------------------------------------------------------------- 2108 2109 if(countG1D.eq.saveG1D)then 2110 2111 call writeg1d(ngrid,1,tsurf,'tsurf','K') 2112 call writeg1d(ngrid,1,ps,'ps','Pa') 2113 call writeg1d(ngrid,nlayer,zt,'T','K') 2114 call writeg1d(ngrid,nlayer,pq(1,1,1),'q','kg/kg') 2115 call writeg1d(ngrid,nlayer,aerosol,'aerosol','SI') 2116 call writeg1d(ngrid,nlayer,reffrad,'reffrad','SI') 2117 call writeg1d(ngrid,nlayer,zdtlw,'dtlw','SI') 2118 call writeg1d(ngrid,nlayer,zdtsw,'dtsw','SI') 2119 call writeg1d(ngrid,nlayer,zdtdyn,'dtdyn','SI') 2120 2121 ! radiation balance 2122 call writeg1d(ngrid,1,ISR/totarea,'ISR','W m-2') 2123 call writeg1d(ngrid,1,ASR/totarea,'ASR','W m-2') 2124 call writeg1d(ngrid,1,OLR/totarea,'OLR','W m-2') 2125 2126 if(tracer) then 2127 do iq=1,nq 2128 call writeg1d(ngrid,1,qsurf(1,iq),trim(noms(iq))//'_s','kg m^-2') 2129 call writeg1d(ngrid,1,qcol(1,iq),trim(noms(iq))//'_c','kg m^-2') 2130 call writeg1d(ngrid,nlayer,zq(1,1,iq),noms(iq),'kg/kg') 2131 end do 2132 2133 if(waterrain)then 2134 CALL writeg1d(ngrid,1,zdqsrain,'rainfall','kg m-2 s-1') 2135 CALL writeg1d(ngrid,1,zdqssnow,'snowfall','kg m-2 s-1') 2136 endif 2137 2138 end if 2139 2140 countG1D=1 2141 else 2142 countG1D=countG1D+1 2143 endif ! if time to save 2144 2145 endif ! if(ngrid.ne.1) 2108 ! output spectrum 2109 if(specOLR.and.corrk)then 2110 call writediagspecIR(ngrid,"OLR3D","OLR(lon,lat,band)","W/m^2/cm^-1",3,OLR_nu) 2111 call writediagspecVI(ngrid,"OSR3D","OSR(lon,lat,band)","W/m^2/cm^-1",3,OSR_nu) 2112 endif 2113 2146 2114 2147 2115 icount=icount+1
Note: See TracChangeset
for help on using the changeset viewer.
