Changeset 1482 for trunk/LMDZ.GENERIC/libf/phystd/physiq.F90
- Timestamp:
- Oct 14, 2015, 5:05:47 PM (9 years ago)
- File:
-
- 1 edited
-
trunk/LMDZ.GENERIC/libf/phystd/physiq.F90 (modified) (18 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/LMDZ.GENERIC/libf/phystd/physiq.F90
r1477 r1482 11 11 use watercommon_h, only : RLVTT, Psat_water,epsi 12 12 use gases_h, only: gnom, gfrac 13 use radcommon_h, only: sigma, glat, grav 13 use radcommon_h, only: sigma, glat, grav, BWNV 14 14 use radii_mod, only: h2o_reffrad, co2_reffrad 15 15 use aerosol_mod, only: iaero_co2, iaero_h2o 16 use surfdat_h, only: phisfi, albedodat,zmea, zstd, zsig, zgam, zthe, &16 use surfdat_h, only: phisfi, zmea, zstd, zsig, zgam, zthe, & 17 17 dryness, watercaptag 18 18 use comdiurn_h, only: coslat, sinlat, coslon, sinlon … … 190 190 ! ---------------------- 191 191 192 integer day_ini ! Initial date of the run (sol since Ls=0). 193 integer icount ! Counter of calls to physiq during the run. 194 real, dimension(:),allocatable,save :: tsurf ! Surface temperature (K). 195 real, dimension(:,:),allocatable,save :: tsoil ! Sub-surface temperatures (K). 196 real, dimension(:),allocatable,save :: albedo ! Surface albedo. 197 198 !$OMP THREADPRIVATE(tsurf,tsoil,albedo) 199 200 real,dimension(:),allocatable,save :: albedo0 ! Initial surface albedo. 201 real,dimension(:),allocatable,save :: rnat ! Defines the type of the grid (ocean,continent,...). By BC. 202 203 !$OMP THREADPRIVATE(albedo0,rnat) 192 integer day_ini ! Initial date of the run (sol since Ls=0). 193 integer icount ! Counter of calls to physiq during the run. 194 real, dimension(:),allocatable,save :: tsurf ! Surface temperature (K). 195 real, dimension(:,:),allocatable,save :: tsoil ! Sub-surface temperatures (K). 196 real, dimension(:,:),allocatable,save :: albedo ! Surface Spectral albedo. By MT2015. 197 real, dimension(:),allocatable,save :: albedo_equivalent ! Spectral Mean albedo. 198 real, dimension(:),allocatable,save :: albedo_snow_SPECTV ! Snow Spectral albedo. 199 real, dimension(:),allocatable,save :: albedo_co2_ice_SPECTV ! CO2 ice Spectral albedo. 200 201 !$OMP THREADPRIVATE(tsurf,tsoil,albedo,albedo_equivalent,albedo_snow_SPECTV,albedo_co2_ice_SPECTV) 202 203 real,dimension(:),allocatable,save :: albedo_bareground ! Bare Ground Albedo. By MT 2015. 204 real,dimension(:),allocatable,save :: rnat ! Defines the type of the grid (ocean,continent,...). By BC. 205 206 !$OMP THREADPRIVATE(albedo_bareground,rnat) 204 207 205 208 real,dimension(:),allocatable,save :: emis ! Thermal IR surface emissivity. … … 228 231 real pw(ngrid,nlayer) ! Vertical velocity (m/s). (NOTE : >0 WHEN DOWNWARDS !!) 229 232 230 integer l,ig,ierr,iq 233 integer l,ig,ierr,iq,nw 231 234 232 235 ! FOR DIAGNOSTIC : 233 236 234 real,dimension(:),allocatable,save :: fluxsurf_lw ! Incident Long Wave (IR) surface flux (W.m-2). 235 real,dimension(:),allocatable,save :: fluxsurf_sw ! incident Short Wave (stellar) surface flux (W.m-2). 236 real,dimension(:),allocatable,save :: fluxtop_lw ! Outgoing LW (IR) flux to space (W.m-2). 237 real,dimension(:),allocatable,save :: fluxabs_sw ! Absorbed SW (stellar) flux (W.m-2). 238 real,dimension(:),allocatable,save :: fluxtop_dn ! Incoming SW (stellar) radiation at the top of the atmosphere (W.m-2). 239 real,dimension(:),allocatable,save :: fluxdyn ! Horizontal heat transport by dynamics (W.m-2). 240 real,dimension(:,:),allocatable,save :: OLR_nu ! Outgoing LW radiation in each band (Normalized to the band width (W/m2/cm-1)). 241 real,dimension(:,:),allocatable,save :: OSR_nu ! Outgoing SW radiation in each band (Normalized to the band width (W/m2/cm-1)). 242 real,dimension(:,:),allocatable,save :: zdtlw ! LW heating tendencies (K/s). 243 real,dimension(:,:),allocatable,save :: zdtsw ! SW heating tendencies (K/s). 244 real,dimension(:),allocatable,save :: sensibFlux ! Turbulent flux given by the atmosphere to the surface (W.m-2). 245 246 !$OMP THREADPRIVATE(fluxsurf_lw,fluxsurf_sw,fluxtop_lw,fluxabs_sw,fluxtop_dn,fluxdyn,OLR_nu,OSR_nu,& 237 real,dimension(:),allocatable,save :: fluxsurf_lw ! Incident Long Wave (IR) surface flux (W.m-2). 238 real,dimension(:),allocatable,save :: fluxsurf_sw ! Incident Short Wave (stellar) surface flux (W.m-2). 239 real,dimension(:),allocatable,save :: fluxsurfabs_sw ! Absorbed Short Wave (stellar) flux by the surface (W.m-2). 240 real,dimension(:),allocatable,save :: fluxtop_lw ! Outgoing LW (IR) flux to space (W.m-2). 241 real,dimension(:),allocatable,save :: fluxabs_sw ! Absorbed SW (stellar) flux (W.m-2). 242 real,dimension(:),allocatable,save :: fluxtop_dn ! Incoming SW (stellar) radiation at the top of the atmosphere (W.m-2). 243 real,dimension(:),allocatable,save :: fluxdyn ! Horizontal heat transport by dynamics (W.m-2). 244 real,dimension(:,:),allocatable,save :: OLR_nu ! Outgoing LW radiation in each band (Normalized to the band width (W/m2/cm-1)). 245 real,dimension(:,:),allocatable,save :: OSR_nu ! Outgoing SW radiation in each band (Normalized to the band width (W/m2/cm-1)). 246 real,dimension(:,:),allocatable,save :: zdtlw ! LW heating tendencies (K/s). 247 real,dimension(:,:),allocatable,save :: zdtsw ! SW heating tendencies (K/s). 248 real,dimension(:),allocatable,save :: sensibFlux ! Turbulent flux given by the atmosphere to the surface (W.m-2). 249 250 !$OMP THREADPRIVATE(fluxsurf_lw,fluxsurf_sw,fluxsurfabs_sw,fluxtop_lw,fluxabs_sw,fluxtop_dn,fluxdyn,OLR_nu,OSR_nu,& 247 251 248 252 !$OMP zdtlw,zdtsw,sensibFlux) … … 376 380 logical clearsky ! For double radiative transfer call. By BC 377 381 378 real fluxsurf_lw1(ngrid), fluxsurf_sw1(ngrid), fluxtop_lw1(ngrid), fluxabs_sw1(ngrid) ! For SW/LW flux diagnostics. 379 real tau_col1(ngrid) ! For aerosol optical depth diagnostic. 380 real OLR_nu1(ngrid,L_NSPECTI), OSR_nu1(ngrid,L_NSPECTV) ! For Outgoing Radiation diagnostics. 382 ! For Clear Sky Case. 383 real fluxsurf_lw1(ngrid), fluxsurf_sw1(ngrid), fluxsurfabs_sw1(ngrid) ! For SW/LW flux. 384 real fluxtop_lw1(ngrid), fluxabs_sw1(ngrid) ! For SW/LW flux. 385 real albedo_equivalent1(ngrid) ! For Equivalent albedo calculation. 386 real tau_col1(ngrid) ! For aerosol optical depth diagnostic. 387 real OLR_nu1(ngrid,L_NSPECTI), OSR_nu1(ngrid,L_NSPECTV) ! For Outgoing Radiation diagnostics. 381 388 real tf, ntf 382 389 … … 435 442 ALLOCATE(tsurf(ngrid)) 436 443 ALLOCATE(tsoil(ngrid,nsoilmx)) 437 ALLOCATE(albedo(ngrid)) 438 ALLOCATE(albedo0(ngrid)) 444 ALLOCATE(albedo(ngrid,L_NSPECTV)) 445 ALLOCATE(albedo_equivalent(ngrid)) 446 ALLOCATE(albedo_snow_SPECTV(L_NSPECTV)) 447 ALLOCATE(albedo_co2_ice_SPECTV(L_NSPECTV)) 448 ALLOCATE(albedo_bareground(ngrid)) 439 449 ALLOCATE(rnat(ngrid)) 440 450 ALLOCATE(emis(ngrid)) … … 457 467 ALLOCATE(fluxsurf_lw(ngrid)) 458 468 ALLOCATE(fluxsurf_sw(ngrid)) 469 ALLOCATE(fluxsurfabs_sw(ngrid)) 459 470 ALLOCATE(fluxtop_lw(ngrid)) 460 471 ALLOCATE(fluxabs_sw(ngrid)) … … 521 532 write (*,*) 'In physiq day_ini =', day_ini 522 533 523 ! Initialize albedo and orbital calculation. 524 ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 525 call surfini(ngrid,nq,qsurf,albedo0) 534 ! Initialize albedo calculation. 535 ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 536 albedo(:,:)=0.0 537 albedo_bareground(:)=0.0 538 albedo_snow_SPECTV(:)=0.0 539 albedo_co2_ice_SPECTV(:)=0.0 540 call surfini(ngrid,nq,qsurf,albedo,albedo_bareground,albedo_snow_SPECTV,albedo_co2_ice_SPECTV) 541 542 543 ! Initialize orbital calculation. 544 ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 526 545 call iniorbit(apoastr,periastr,year_day,peri_day,obliquit) 527 546 528 albedo(:)=albedo0(:)529 547 530 548 if(tlocked)then … … 677 695 if (ngrid.ne.1) then ! Note : no need to create a restart file in 1d. 678 696 call physdem0("restartfi.nc",long,lati,nsoilmx,ngrid,nlayer,nq, & 679 ptimestep,pday+nday,time_phys,area, &680 albedo dat,inertiedat,zmea,zstd,zsig,zgam,zthe)697 ptimestep,pday+nday,time_phys,area, & 698 albedo_bareground,inertiedat,zmea,zstd,zsig,zgam,zthe) 681 699 endif 682 700 701 683 702 endif ! end of 'firstcall' 684 703 … … 848 867 ! standard callcorrk 849 868 clearsky=.false. 850 call callcorrk(ngrid,nlayer,pq,nq,qsurf, & 851 albedo,emis,mu0,pplev,pplay,pt, & 852 tsurf,fract,dist_star,aerosol,muvar, & 853 zdtlw,zdtsw,fluxsurf_lw,fluxsurf_sw,fluxtop_lw, & 854 fluxabs_sw,fluxtop_dn,OLR_nu,OSR_nu, & 855 tau_col,cloudfrac,totcloudfrac, & 856 clearsky,firstcall,lastcall) 857 858 ! Option to call scheme once more for clear regions 869 call callcorrk(ngrid,nlayer,pq,nq,qsurf, & 870 albedo,albedo_equivalent,emis,mu0,pplev,pplay,pt, & 871 tsurf,fract,dist_star,aerosol,muvar, & 872 zdtlw,zdtsw,fluxsurf_lw,fluxsurf_sw, & 873 fluxsurfabs_sw,fluxtop_lw, & 874 fluxabs_sw,fluxtop_dn,OLR_nu,OSR_nu, & 875 tau_col,cloudfrac,totcloudfrac, & 876 clearsky,firstcall,lastcall) 877 878 ! Option to call scheme once more for clear regions 859 879 if(CLFvarying)then 860 880 861 881 ! ---> PROBLEMS WITH ALLOCATED ARRAYS : temporary solution in callcorrk: do not deallocate if CLFvarying ... 862 882 clearsky=.true. 863 call callcorrk(ngrid,nlayer,pq,nq,qsurf, & 864 albedo,emis,mu0,pplev,pplay,pt, & 865 tsurf,fract,dist_star,aerosol,muvar, & 866 zdtlw1,zdtsw1,fluxsurf_lw1,fluxsurf_sw1,fluxtop_lw1, & 867 fluxabs_sw1,fluxtop_dn,OLR_nu1,OSR_nu1, & 868 tau_col1,cloudfrac,totcloudfrac, & 883 call callcorrk(ngrid,nlayer,pq,nq,qsurf, & 884 albedo,albedo_equivalent1,emis,mu0,pplev,pplay,pt, & 885 tsurf,fract,dist_star,aerosol,muvar, & 886 zdtlw1,zdtsw1,fluxsurf_lw1,fluxsurf_sw1, & 887 fluxsurfabs_sw1,fluxtop_lw1, & 888 fluxabs_sw1,fluxtop_dn,OLR_nu1,OSR_nu1, & 889 tau_col1,cloudfrac,totcloudfrac, & 869 890 clearsky,firstcall,lastcall) 870 891 clearsky = .false. ! just in case. 871 872 892 873 893 ! Sum the fluxes and heating rates from cloudy/clear cases 874 894 do ig=1,ngrid 875 895 tf=totcloudfrac(ig) 876 ntf=1.-tf 877 878 fluxsurf_lw(ig) = ntf*fluxsurf_lw1(ig) + tf*fluxsurf_lw(ig) 879 fluxsurf_sw(ig) = ntf*fluxsurf_sw1(ig) + tf*fluxsurf_sw(ig) 880 fluxtop_lw(ig) = ntf*fluxtop_lw1(ig) + tf*fluxtop_lw(ig) 881 fluxabs_sw(ig) = ntf*fluxabs_sw1(ig) + tf*fluxabs_sw(ig) 882 tau_col(ig) = ntf*tau_col1(ig) + tf*tau_col(ig) 896 ntf=1.-tf 897 fluxsurf_lw(ig) = ntf*fluxsurf_lw1(ig) + tf*fluxsurf_lw(ig) 898 fluxsurf_sw(ig) = ntf*fluxsurf_sw1(ig) + tf*fluxsurf_sw(ig) 899 albedo_equivalent(ig) = ntf*albedo_equivalent1(ig) + tf*albedo_equivalent(ig) 900 fluxsurfabs_sw(ig) = ntf*fluxsurfabs_sw1(ig) + tf*fluxsurfabs_sw(ig) 901 fluxtop_lw(ig) = ntf*fluxtop_lw1(ig) + tf*fluxtop_lw(ig) 902 fluxabs_sw(ig) = ntf*fluxabs_sw1(ig) + tf*fluxabs_sw(ig) 903 tau_col(ig) = ntf*tau_col1(ig) + tf*tau_col(ig) 883 904 884 905 zdtlw(ig,1:nlayer) = ntf*zdtlw1(ig,1:nlayer) + tf*zdtlw(ig,1:nlayer) … … 887 908 OSR_nu(ig,1:L_NSPECTV) = ntf*OSR_nu1(ig,1:L_NSPECTV) + tf*OSR_nu(ig,1:L_NSPECTV) 888 909 OLR_nu(ig,1:L_NSPECTI) = ntf*OLR_nu1(ig,1:L_NSPECTI) + tf*OLR_nu(ig,1:L_NSPECTI) 889 enddo910 enddo 890 911 891 912 endif ! end of CLFvarying. … … 896 917 897 918 898 ! Radiative flux from the sky absorbed by the surface (W.m-2) 919 ! Radiative flux from the sky absorbed by the surface (W.m-2). 899 920 GSR=0.0 900 fluxrad_sky(1:ngrid)=emis(1:ngrid)*fluxsurf_lw(1:ngrid)+fluxsurf _sw(1:ngrid)*(1.-albedo(1:ngrid))921 fluxrad_sky(1:ngrid)=emis(1:ngrid)*fluxsurf_lw(1:ngrid)+fluxsurfabs_sw(1:ngrid) 901 922 902 923 !if(noradsurf)then ! no lower surface; SW flux just disappears … … 910 931 911 932 elseif(newtonian)then 912 913 ! II.b Call Newtonian cooling scheme 914 ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 933 934 ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 935 ! II.b Call Newtonian cooling scheme 936 ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 915 937 call newtrelax(ngrid,nlayer,mu0,sinlat,zpopsk,pt,pplay,pplev,dtrad,firstcall) 916 938 … … 928 950 endif 929 951 fluxsurf_sw(1:ngrid) = fluxtop_dn(1:ngrid) 930 fluxrad_sky(1:ngrid) = fluxtop_dn(1:ngrid)*(1.-albedo(1:ngrid)) 952 print*,'------------WARNING---WARNING------------' ! by MT2015. 953 print*,'You are in corrk=false mode, ' 954 print*,'and the surface albedo is taken equal to its value at 0.5 micrometers ...' 955 if ( (10000./BWNV(1) .le. 0.5) .or. (10000./BWNV(L_NSPECTV) .ge. 0.5) ) then 956 print*,'0.5 microns band doesnt match your visible bands ! Abort !' 957 call abort 958 endif 959 DO nw=1, L_NSPECTV-1 960 if ( (10000./BWNV(nw) .ge. 0.5) .and. (10000./BWNV(nw+1) .lt. 0.5) ) then 961 fluxsurfabs_sw(1:ngrid) = fluxtop_dn(1:ngrid)*(1.-albedo(1:ngrid,nw)) 962 fluxrad_sky(1:ngrid) = fluxsurfabs_sw(1:ngrid) 963 endif 964 ENDDO 931 965 fluxtop_lw(1:ngrid) = emis(1:ngrid)*sigma*tsurf(1:ngrid)**4 932 966 … … 950 984 call planetwide_sumval(cpp*massarea(:,:)*zdtsw(:,:)/totarea_planet,dEtotSW) 951 985 call planetwide_sumval(cpp*massarea(:,:)*zdtlw(:,:)/totarea_planet,dEtotLW) 952 call planetwide_sumval(fluxsurf_sw(:)*(1.-albedo(:))*area(:)/totarea_planet,dEtotsSW) !JL13 carefull, albedo can have changed since the last time we called corrk 986 !call planetwide_sumval(fluxsurf_sw(:)*(1.-albedo_equivalent(:))*area(:)/totarea_planet,dEtotsSW) !JL13 carefull, albedo can have changed since the last time we called corrk 987 call planetwide_sumval(fluxsurfabs_sw(:)*area(:)/totarea_planet,dEtotsSW) !JL13 carefull, albedo can have changed since the last time we called corrk 953 988 call planetwide_sumval((fluxsurf_lw(:)*emis(:)-zplanck(:))*area(:)/totarea_planet,dEtotsLW) 954 989 dEzRadsw(:,:)=cpp*mass(:,:)*zdtsw(:,:) … … 1166 1201 pphi,pdt,pdu,pdv,zdtsurf,pu,pv,pq,pdq, & 1167 1202 qsurf(1,igcm_co2_ice),albedo,emis, & 1203 albedo_bareground,albedo_co2_ice_SPECTV, & 1168 1204 zdtc,zdtsurfc,pdpsrf,zduc,zdvc, & 1169 1205 zdqc) … … 1461 1497 if(hydrology)then 1462 1498 1463 call hydrol(ngrid,nq,ptimestep,rnat,tsurf,qsurf,dqsurf,dqs_hyd, & 1464 capcal,albedo0,albedo,mu0,zdtsurf,zdtsurf_hyd,hice,pctsrf_sic, & 1499 call hydrol(ngrid,nq,ptimestep,rnat,tsurf,qsurf,dqsurf,dqs_hyd, & 1500 capcal,albedo,albedo_bareground, & 1501 albedo_snow_SPECTV,albedo_co2_ice_SPECTV, & 1502 mu0,zdtsurf,zdtsurf_hyd,hice,pctsrf_sic, & 1465 1503 sea_ice) 1466 1504 … … 1848 1886 call wstats(ngrid,"ASR","absorbed stellar rad.","W m-2",2,fluxabs_sw) 1849 1887 call wstats(ngrid,"OLR","outgoing longwave rad.","W m-2",2,fluxtop_lw) 1850 call wstats(ngrid,"ALB","Surface albedo"," ",2,albedo) 1888 !call wstats(ngrid,"ALB","Surface albedo"," ",2,albedo_equivalent) 1889 !call wstats(ngrid,"ALB_1st","First Band Surface albedo"," ",2,albedo(:,1)) 1851 1890 call wstats(ngrid,"p","Pressure","Pa",3,pplay) 1852 1891 call wstats(ngrid,"temp","Atmospheric temperature","K",3,zt) … … 1955 1994 if(callrad.and.(.not.newtonian))then 1956 1995 1957 call writediagfi(ngrid,"ALB","Surface albedo"," ",2,albedo) 1996 !call writediagfi(ngrid,"ALB","Surface albedo"," ",2,albedo_equivalent) 1997 !call writediagfi(ngrid,"ALB_1st","First Band Surface albedo"," ",2,albedo(:,1)) 1958 1998 call writediagfi(ngrid,"ISR","incoming stellar rad.","W m-2",2,fluxtop_dn) 1959 1999 call writediagfi(ngrid,"ASR","absorbed stellar rad.","W m-2",2,fluxabs_sw)
Note: See TracChangeset
for help on using the changeset viewer.
