Changeset 3241 for trunk/LMDZ.PLUTO/libf/phypluto/physiq_mod.F90
- Timestamp:
- Feb 27, 2024, 11:40:35 PM (9 months ago)
- File:
-
- 1 edited
-
trunk/LMDZ.PLUTO/libf/phypluto/physiq_mod.F90 (modified) (12 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/LMDZ.PLUTO/libf/phypluto/physiq_mod.F90
r3237 r3241 29 29 use comsaison_h, only: mu0, fract, dist_star, declin, right_ascen 30 30 use comsoil_h, only: nsoilmx, layer, mlayer, inertiedat 31 use geometry_mod, only: latitude, longitude, cell_area 31 use geometry_mod, only: latitude, longitude, cell_area, real_area 32 32 USE comgeomfi_h, only: totarea, totarea_planet 33 33 USE tracer_h, only: noms, mmol, radius, rho_q, qext, & … … 43 43 use mod_phys_lmdz_para, only : is_master 44 44 use planete_mod, only: apoastr, periastr, year_day, peri_day, & 45 obliquit, nres, z0 45 obliquit, nres, z0, adjust, tpal 46 46 use comcstfi_mod, only: pi, g, rcp, r, rad, mugaz, cpp 47 47 use time_phylmdz_mod, only: daysec … … 53 53 lwrite, mass_redistrib, meanOLR, & 54 54 fast,fasthaze,haze,metcloud,monoxcloud,& 55 n2cond,nearn2cond, noseason_day, & 55 n2cond,nearn2cond,noseason_day,conservn2, & 56 kbo,triton,paleo,paleoyears, & 56 57 fast, carbox, methane, UseVdifcPlutold, & 57 aerohaze, haze_proffix, triton, paleo, &58 aerohaze,haze_proffix,source_haze, & 58 59 season, sedimentation,generic_condensation, & 59 60 specOLR, & … … 254 255 !$OMP THREADPRIVATE(day_ini,icount) 255 256 257 !Pluto specific 256 258 REAL,save :: acond,bcond 257 REAL, save:: tcond1p4Pa259 REAL,save :: tcond1p4Pa 258 260 DATA tcond1p4Pa/38/ 261 259 262 260 263 261 264 ! Local variables : 262 265 ! ----------------- 266 ! Tendencies for the paleoclimate mode 267 REAL qsurfyear(ngrid,nq) ! kg.m-2 averaged mass of ice lost/gained in the last Pluto year of the run 268 REAL phisfinew(ngrid) ! geopotential of the bedrock (= phisfi-qsurf/1000*g) 269 REAL qsurfpal(ngrid,nq) ! qsurf after a paleoclimate step : for physdem1 and restartfi 270 REAL phisfipal(ngrid) ! geopotential after a paleoclimate step : for physdem1 and restartfi 271 REAL oblipal ! change of obliquity 272 REAL peri_daypal ! new periday 273 REAL eccpal ! change of eccentricity 274 REAL tpalnew ! change of time 275 REAL adjustnew ! change in N2 ice albedo 276 REAL pdaypal ! new pday = day_ini + step 277 REAL zdt_tot ! time range corresponding to the flux of qsurfyear 278 REAL massacc(nq) ! accumulated mass 279 REAL masslost(nq) ! accumulated mass 280 281 REAL globave ! globalaverage 2D ps 282 REAL globaveice(nq) ! globalaverage 2D ice 283 REAL globavenewice(nq) ! globalaverage 2D ice 284 INTEGER lecttsoil ! lecture of tsoil from proftsoil 285 REAL qsurf1(ngrid,nq) ! saving qsurf to calculate flux over long timescales kg.m-2 286 REAL flusurf(ngrid,nq) ! flux cond/sub kg.m-2.s-1 287 REAL flusurfold(ngrid,nq) ! old flux cond/sub kg.m-2.s-1 288 289 263 290 264 291 ! Aerosol (dust or ice) extinction optical depth at reference wavelength … … 271 298 real omega(ngrid,nlayer) ! omega velocity (Pa/s, >0 when downward) 272 299 273 integer l,ig,ierr,iq,nw,isoil,iesp, igcm_generic_vap, igcm_generic_ice300 integer i,l,ig,ierr,iq,nw,isoil,iesp, igcm_generic_vap, igcm_generic_ice 274 301 logical call_ice_vap_generic ! to call only one time the ice/vap pair of a tracer 275 302 … … 1514 1541 endif! end of if 'tracer' 1515 1542 1543 if (conservn2) then 1544 write(*,*) 'conservn2 tracer' 1545 ! call testconservmass(ngrid,nlayer,pplev(:,1)+ & 1546 ! pdpsrf(:)*ptimestep,qsurf(:,1)) 1547 endif 1548 1549 DO ig=1,ngrid 1550 flusurf(ig,igcm_n2)=(qsurf(ig,igcm_n2)- & 1551 qsurf1(ig,igcm_n2))/ptimestep 1552 flusurfold(ig,igcm_n2)=flusurf(ig,igcm_n2) 1553 if (methane) then 1554 flusurf(ig,igcm_ch4_ice)=(qsurf(ig,igcm_ch4_ice)- & 1555 qsurf1(ig,igcm_ch4_ice))/ptimestep 1556 flusurfold(ig,igcm_ch4_ice)=flusurf(ig,igcm_ch4_ice) 1557 endif 1558 if (carbox) then 1559 flusurf(ig,igcm_co_ice)=(qsurf(ig,igcm_co_ice)- & 1560 qsurf1(ig,igcm_co_ice))/ptimestep 1561 !flusurfold(ig,igcm_co_ice)=flusurf(ig,igcm_co_ice) 1562 endif 1563 ENDDO 1564 1565 !! Special source of haze particle ! 1566 ! todo: should be placed in haze and use tendency of n2 instead of flusurf 1567 IF (source_haze) THEN 1568 ! call hazesource(ngrid,nlayer,nq,ptimestep, & 1569 ! pplev,flusurf,mu0,zdq_source) 1570 1571 DO iq=1, nq 1572 DO l=1,nlayer 1573 DO ig=1,ngrid 1574 pdq(ig,l,iq)=pdq(ig,l,iq)+zdq_source(ig,l,iq) 1575 ENDDO 1576 ENDDO 1577 ENDDO 1578 ENDIF 1579 1516 1580 1517 1581 !------------------------------------------------ 1518 1582 ! VII. Surface and sub-surface soil temperature 1519 1583 !------------------------------------------------ 1584 1585 ! For diagnostic 1586 ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 1587 DO ig=1,ngrid 1588 rho(ig,1) = pplay(ig,1)/(r*pt(ig,1)) 1589 sensiblehf1(ig)=rho(ig,1)*cpp*(0.4/log(zzlay(ig,1)/z0))**2* & 1590 (pu(ig,1)*pu(ig,1)+pv(ig,1)*pv(ig,1))**0.5* & 1591 (tsurf(ig)-pt(ig,1)) 1592 sensiblehf2(ig)=zflubid(ig)-capcal(ig)*zdtsdif(ig) 1593 1594 ENDDO 1520 1595 1521 1596 … … 1702 1777 1703 1778 1779 if (methane) then 1780 IF (fast) then ! zq is the mixing ratio supposingly mixed in all atmosphere 1781 DO ig=1,ngrid 1782 vmr_ch4(ig)=zq(ig,1,igcm_ch4_gas)* & 1783 mmol(igcm_n2)/mmol(igcm_ch4_gas)*100. 1784 ENDDO 1785 ELSE 1786 DO ig=1,ngrid 1787 ! compute vmr methane 1788 vmr_ch4(ig)=qcol(ig,igcm_ch4_gas)* & 1789 g/ps(ig)*mmol(igcm_n2)/mmol(igcm_ch4_gas)*100. 1790 ! compute density methane 1791 DO l=1,nlayer 1792 zrho_ch4(ig,l)=zq(ig,l,igcm_ch4_gas)*rho(ig,l) 1793 ENDDO 1794 ENDDO 1795 ENDIF 1796 endif 1797 1798 if (carbox) then 1799 IF (fast) then 1800 DO ig=1,ngrid 1801 vmr_co(ig)=zq(ig,1,igcm_co_gas)* & 1802 mmol(igcm_n2)/mmol(igcm_co_gas)*100. 1803 ENDDO 1804 ELSE 1805 DO ig=1,ngrid 1806 ! compute vmr CO 1807 vmr_co(ig)=qcol(ig,igcm_co_gas)* & 1808 g/ps(ig)*mmol(igcm_n2)/mmol(igcm_co_gas)*100. 1809 ! compute density CO 1810 DO l=1,nlayer 1811 zrho_co(ig,l)=zq(ig,l,igcm_co_gas)*rho(ig,l) 1812 ENDDO 1813 ENDDO 1814 ENDIF 1815 endif 1816 1817 zrho_haze(:,:)=0. 1818 zdqrho_photprec(:,:)=0. 1819 IF (haze.and.aerohaze) then 1820 DO ig=1,ngrid 1821 DO l=1,nlayer 1822 zrho_haze(ig,l)=zq(ig,l,igcm_haze)*rho(ig,l) 1823 zdqrho_photprec(ig,l)=zdqphot_prec(ig,l)*rho(ig,l) 1824 ENDDO 1825 ENDDO 1826 ENDIF 1827 1828 IF (fasthaze) then 1829 DO ig=1,ngrid 1830 qcol(ig,igcm_haze)=zq(ig,1,igcm_haze)*pplev(ig,1)/g 1831 qcol(ig,igcm_prec_haze)=zq(ig,1,igcm_prec_haze)*pplev(ig,1)/g 1832 ENDDO 1833 ENDIF 1834 1835 ! Info about Ls, declin... 1836 IF (fast) THEN 1837 if (is_master) write(*,*),'Ls=',zls*180./pi,' dec=',declin*180./pi 1838 if (is_master) write(*,*),'zday=',zday,' ps=',globave 1839 IF (lastcall) then 1840 if (is_master) write(*,*),'lastcall' 1841 ENDIF 1842 ELSE 1843 if (is_master) write(*,*),'Ls=',zls*180./pi,'decli=',declin*180./pi,'zday=',zday 1844 ENDIF 1845 1846 lecttsoil=0 ! default value for lecttsoil 1847 call getin_p("lecttsoil",lecttsoil) 1848 IF (lastcall.and.(ngrid.EQ.1).and.(lecttsoil.eq.1)) THEN 1849 ! save tsoil temperature profile for 1D profile 1850 OPEN(13,file='proftsoil.out',form='formatted') 1851 DO i=1,nsoilmx 1852 write(13,*) tsoil(1,i) 1853 ENDDO 1854 CLOSE(13) 1855 ENDIF 1856 1857 1704 1858 if (is_master) print*,'--> Ls =',zls*180./pi 1705 1859 … … 1724 1878 1725 1879 ! endif 1880 if (paleo) then 1881 ! time range for tendencies of ice flux qsurfyear 1882 zdt_tot=year_day ! Last year of simulation 1883 1884 masslost(:)=0. 1885 massacc(:)=0. 1886 1887 DO ig=1,ngrid 1888 ! update new reservoir of ice on the surface 1889 DO iq=1,nq 1890 ! kg/m2 to be sublimed or condensed during paleoyears 1891 qsurfyear(ig,iq)=qsurfyear(ig,iq)* & 1892 paleoyears*365.25/(zdt_tot*daysec/86400.) 1893 1894 ! special case if we sublime the entire reservoir 1895 IF (-qsurfyear(ig,iq).gt.qsurf(ig,iq)) THEN 1896 masslost(iq)=masslost(iq)+real_area(ig)* & 1897 (-qsurfyear(ig,iq)-qsurf(ig,iq)) 1898 qsurfyear(ig,iq)=-qsurf(ig,iq) 1899 ENDIF 1900 1901 IF (qsurfyear(ig,iq).gt.0.) THEN 1902 massacc(iq)=massacc(iq)+real_area(ig)*qsurfyear(ig,iq) 1903 ENDIF 1904 1905 ENDDO 1906 ENDDO 1907 1908 DO ig=1,ngrid 1909 DO iq=1,nq 1910 qsurfpal(ig,iq)=qsurf(ig,iq)+qsurfyear(ig,iq) 1911 IF (qsurfyear(ig,iq).gt.0.) THEN 1912 qsurfpal(ig,iq)=qsurfpal(ig,iq)- & 1913 qsurfyear(ig,iq)*masslost(iq)/massacc(iq) 1914 ENDIF 1915 ENDDO 1916 ENDDO 1917 ! Finally ensure conservation of qsurf 1918 DO iq=1,nq 1919 call globalaverage2d(ngrid,qsurf(:,iq),globaveice(iq)) 1920 call globalaverage2d(ngrid,qsurfpal(:,iq), & 1921 globavenewice(iq)) 1922 IF (globavenewice(iq).gt.0.) THEN 1923 qsurfpal(:,iq)=qsurfpal(:,iq)* & 1924 globaveice(iq)/globavenewice(iq) 1925 ENDIF 1926 ENDDO 1927 1928 ! update new geopotential depending on the ice reservoir 1929 phisfipal(:)=phisfinew(:)+qsurfpal(:,igcm_n2)*g/1000. 1930 !phisfipal(ig)=phisfi(ig) 1931 1932 1933 if (kbo.or.triton) then ! case of Triton : we do not change the orbital parameters 1934 1935 pdaypal=pday ! no increment of pdaypal to keep same evolution of the subsolar point 1936 eccpal=1.-periastr/((periastr+apoastr)/2.) !no change of ecc 1937 peri_daypal=peri_day ! no change 1938 oblipal=obliquit ! no change 1939 tpalnew=tpal 1940 adjustnew=adjust 1941 1942 else ! Pluto 1943 ! update new pday and tpal (Myr) to be set in startfi controle 1944 pdaypal=int(day_ini+paleoyears*365.25/6.3872) 1945 tpalnew=tpal+paleoyears*1.e-6 ! Myrs 1946 1947 ! update new N2 ice adjustment (not tested yet on Pluto) 1948 adjustnew=adjust 1949 1950 ! update milankovitch parameters : obliquity,Lsp,ecc 1951 call calcmilank(tpalnew,oblipal,peri_daypal,eccpal) 1952 !peri_daypal=peri_day 1953 !eccpal=0.009 1954 1955 endif 1956 1957 if (is_master) write(*,*) "Paleo peri=",peri_daypal," tpal=",tpalnew 1958 if (is_master) write(*,*) "Paleo eccpal=",eccpal," tpal=",tpalnew 1959 1960 1726 1961 #ifndef MESOSCALE 1962 ! create restartfi 1963 if (ngrid.ne.1) then 1964 !TODO: import this routine from pluto.old 1965 ! call physdem1pal("restartfi.nc",long,lati,nsoilmx,nq, & 1966 ! ptimestep,pdaypal, & 1967 ! ztime_fin,tsurf,tsoil,emis,q2,qsurfpal, & 1968 ! cell_area,albedodat,tidat,zmea,zstd,zsig, & 1969 ! zgam,zthe,oblipal,eccpal,tpalnew,adjustnew,phisfipal, & 1970 ! peri_daypal) 1971 endif 1972 else ! 'paleo' 1727 1973 1728 1974 if (ngrid.ne.1) then … … 1734 1980 endif 1735 1981 #endif 1982 endif ! end of 'paleo' 1736 1983 ! if(ok_slab_ocean) then 1737 1984 ! call ocean_slab_final!(tslab, seaice) … … 1820 2067 call writediagfi(ngrid,"ps","Surface pressure","Pa",2,ps) 1821 2068 1822 if (.not.fast) then 2069 !! Pluto outputs 2070 call writediagfi(ngrid,"rice_ch4","ch4 ice mass mean radius","m",3,rice_ch4) 2071 call writediagfi(ngrid,"dist_star","dist_star","AU",0,dist_star) 2072 2073 if (fast) then 2074 call writediagfi(ngrid,"globave","surf press","Pa",0,globave) 2075 !AF: TODO which outputs? 2076 else 1823 2077 call writediagfi(ngrid,"temp","temperature","K",3,zt) 1824 2078 call writediagfi(ngrid,"teta","potential temperature","K",3,zh) … … 1915 2169 call writediagfi(ngrid,"zdtdyn","Dyn. heating","T s-1",3,zdtdyn) 1916 2170 2171 !Pluto specific 2172 call writediagfi(ngrid,"zdtc","tendancy T cond N2","K",3,zdtc) 2173 ! call writediagfi(ngrid,"zdtch4cloud","tendancy T ch4cloud","K",3,zdtch4cloud) 2174 ! call writediagfi(ngrid,"zdtcocloud","tendancy T cocloud","K",3,zdtcocloud) 2175 ! call writediagfi(ngrid,"zq1temp_ch4"," "," ",2,zq1temp_ch4) 2176 ! call writediagfi(ngrid,"qsat_ch4"," "," ",2,qsat_ch4) 2177 ! call writediagfi(ngrid,"qsat_ch4_l1"," "," ",2,qsat_ch4_l1) 2178 ! call writediagfi(ngrid,"senshf1","senshf1"," ",2,sensiblehf1) 2179 ! call writediagfi(ngrid,"senshf2","senshf2"," ",2,sensiblehf2) 2180 2181 1917 2182 ! For Debugging. 1918 2183 !call writediagfi(ngrid,'rnat','Terrain type',' ',2,real(rnat)) … … 1954 2219 enddo ! end of 'nq' loop 1955 2220 1956 endif ! end of 'tracer' 2221 !Pluto specific 2222 call writediagfi(ngrid,'n2_iceflux','n2_iceflux',"kg m^-2 s^-1",2,flusurf(1,igcm_n2) ) 2223 ! call writediagfi(ngrid,'haze_reff','haze_reff','m',3,reffrad(1,1,1)) 2224 if (methane) then 2225 call writediagfi(ngrid,'ch4_iceflux','ch4_iceflux',& 2226 "kg m^-2 s^-1",2,flusurf(1,igcm_ch4_ice) ) 2227 call writediagfi(ngrid,"vmr_ch4","vmr_ch4","%",2,vmr_ch4) 2228 if (.not.fast) then 2229 call writediagfi(ngrid,"zrho_ch4","zrho_ch4","kg.m-3",3,zrho_ch4(:,:)) 2230 endif 2231 2232 ! Tendancies 2233 call writediagfi(ngrid,"zdqch4cloud","ch4 cloud","T s-1",& 2234 3,zdqch4cloud(1,1,igcm_ch4_gas)) 2235 call writediagfi(ngrid,"zdqcn2_ch4","zdq condn2 ch4","",& 2236 3,zdqc(:,:,igcm_ch4_gas)) 2237 call writediagfi(ngrid,"zdqdif_ch4","zdqdif ch4","",& 2238 3,zdqdif(:,:,igcm_ch4_gas)) 2239 call writediagfi(ngrid,"zdqsdif_ch4_ice","zdqsdif ch4","",& 2240 2,zdqsdif(:,igcm_ch4_ice)) 2241 call writediagfi(ngrid,"zdqadj_ch4","zdqadj ch4","",& 2242 3,zdqadj(:,:,igcm_ch4_gas)) 2243 call writediagfi(ngrid,"zdqsed_ch4","zdqsed ch4","",& 2244 3,zdqsed(:,:,igcm_ch4_gas)) 2245 call writediagfi(ngrid,"zdqssed_ch4","zdqssed ch4","",& 2246 2,zdqssed(:,igcm_ch4_gas)) 2247 call writediagfi(ngrid,"zdtch4cloud","ch4 cloud","T s-1",& 2248 3,zdtch4cloud) 2249 2250 endif 2251 2252 if (carbox) then 2253 call writediagfi(ngrid,'co_iceflux','co_iceflux',& 2254 "kg m^-2 s^-1",2,flusurf(1,igcm_co_ice) ) 2255 call writediagfi(ngrid,"vmr_co","vmr_co","%",2,vmr_co) 2256 if (.not.fast) THEN 2257 call writediagfi(ngrid,"zrho_co","zrho_co","kg.m-3",3,zrho_co(:,:)) 2258 endif 2259 endif 2260 2261 if (haze) then 2262 ! call writediagfi(ngrid,"zrho_haze","zrho_haze","kg.m-3",3,zrho_haze(:,:)) 2263 call writediagfi(ngrid,"zdqrho_photprec","zdqrho_photprec",& 2264 "kg.m-3.s-1",3,zdqrho_photprec(:,:)) 2265 call writediagfi(ngrid,"zdqphot_prec","zdqphot_prec","",& 2266 3,zdqphot_prec(:,:)) 2267 call writediagfi(ngrid,"zdqhaze_ch4","zdqhaze_ch4","",& 2268 3,zdqhaze(:,:,igcm_ch4_gas)) 2269 call writediagfi(ngrid,"zdqhaze_prec","zdqhaze_prec","",& 2270 3,zdqhaze(:,:,igcm_prec_haze)) 2271 if (igcm_haze.ne.0) then 2272 call writediagfi(ngrid,"zdqhaze_haze","zdqhaze_haze","",& 2273 3,zdqhaze(:,:,igcm_haze)) 2274 call writediagfi(ngrid,"zdqssed_haze","zdqssed haze",& 2275 "kg/m2/s",2,zdqssed(:,igcm_haze)) 2276 endif 2277 call writediagfi(ngrid,"zdqphot_ch4","zdqphot_ch4","",& 2278 3,zdqphot_ch4(:,:)) 2279 call writediagfi(ngrid,"zdqconv_prec","zdqconv_prec","",& 2280 3,zdqconv_prec(:,:)) 2281 ! call writediagfi(ngrid,"zdqhaze_col","zdqhaze col","kg/m2/s", 2282 ! & 2,zdqhaze_col(:)) 2283 endif 2284 2285 if (aerohaze) then 2286 call writediagfi(ngrid,"tau_col",& 2287 "Total aerosol optical depth","opacity",2,tau_col) 2288 endif 2289 2290 endif ! end of 'tracer' 1957 2291 1958 2292
Note: See TracChangeset
for help on using the changeset viewer.
