Ignore:
Timestamp:
Mar 9, 2023, 9:38:24 AM (22 months ago)
Author:
romain.vande
Message:

Mars PEM:
New Boolean options for following orbital parameters of ob_ex_lsp.asc: var_obl, var_ex, var_lsp.
If using evol_orbit_pem=true, you can specify which parameter to follow.
True by default: Do you want to change the parameter XXX during the PEM run as prescribed in ob_ex_lsp.asc.
If false, it is set to constant (to the value of the tab_cntrl in the start)
RV

Location:
trunk/LMDZ.MARS/libf/dynphy_lonlat/phymars
Files:
1 added
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/LMDZ.MARS/libf/dynphy_lonlat/phymars/newstart.F

    r2832 r2909  
    2929      use surfdat_h, only: phisfi, z0, zmea, zstd, zsig, zgam, zthe,
    3030     &                     albedodat, z0_default, qsurf, tsurf,
    31      &                     emis, hmons, summit, base, watercap
    32       use comsoil_h, only: inertiedat, layer, mlayer, nsoilmx, tsoil
     31     &                     emis, hmons, summit, base, watercap,
     32     &               ini_surfdat_h_slope_var,end_surfdat_h_slope_var
     33      use comsoil_h, only: inertiedat, layer, mlayer, nsoilmx, tsoil,
     34     & ini_comsoil_h_slope_var, end_comsoil_h_slope_var
    3335      use control_mod, only: day_step, iphysiq, anneeref, planet_type
    3436      use geometry_mod, only: longitude,latitude,cell_area
    … …  
    3638      use phyredem, only: physdem0, physdem1
    3739      use iostart, only: open_startphy
    38       use dimradmars_mod, only: albedo
     40      use dimradmars_mod, only: albedo,
     41     & ini_dimradmars_mod_slope_var,end_dimradmars_mod_slope_var
    3942      use dust_param_mod, only: tauscaling
    4043      use turb_mod, only: q2, wstar
    … …  
    5053      USE exner_hyb_m, ONLY: exner_hyb
    5154      USE inichim_newstart_mod, ONLY: inichim_newstart
    52 
     55      use comslope_mod, ONLY: nslope,def_slope,def_slope_mean,
     56     &             subslope_dist,end_comslope_h,ini_comslope_h
    5357      implicit none
    5458
    … …  
    180184      real,allocatable :: coefvmr(:)  ! Correction coefficient when changing composition
    181185      real :: maxq
    182       integer :: iloc(1), iqmax
     186      integer :: iloc(1), iqmax, islope
    183187! sub-grid cloud fraction
    184188      real :: totcloudfrac(ngridmx)
     189!Variable to change the number of subslope
     190      REAL, ALLOCATABLE :: default_def_slope(:)
     191      REAL,ALLOCATABLE :: tsurf_old_slope(:,:)   ! Surface temperature (K)
     192      REAL,ALLOCATABLE :: emis_old_slope(:,:)    ! Thermal IR surface emissivity
     193      REAL,ALLOCATABLE :: qsurf_old_slope(:,:,:) ! tracer on surface (e.g. kg.m-2)
     194      REAL,ALLOCATABLE :: watercap_old_slope(:,:) ! Surface water ice (kg.m-2)
     195      REAL,ALLOCATABLE :: tsoil_old_slope(:,:,:)
     196      REAL,ALLOCATABLE :: albedo_old_slope(:,:,:) ! Surface albedo in each solar band
     197      integer :: iflat
     198      integer :: nslope_old, nslope_new
     199
     200
    185201
    186202c sortie visu pour les champs dynamiques
    … …  
    443459     &        day_ini,time,tsurf,tsoil,albedo,emis,
    444460     &        q2,qsurf,tauscaling,totcloudfrac,
    445      &        wstar,watercap)
     461     &        wstar,watercap,def_slope,def_slope_mean,subslope_dist)
    446462       
    447463        ! copy albedo and soil thermal inertia
    … …  
    536552      write(*,*) 'mons_ice     : put underground ice layer according
    537553     $ to MONS derived data'
     554      write(*,*) 'nslope       : Change the number of subgrid scale
     555     $ slope'
    538556
    539557        write(*,*)
    … …  
    633651          write(*,*) ' new value of the subsurface temperature:',tsud
    634652c         nouvelle temperature sous la calotte permanente
    635           do l=2,nsoilmx
    636                tsoil(ngridmx,l) =  tsud
    637           end do
     653          do islope=1,nslope
     654            do l=2,nsoilmx
     655               tsoil(ngridmx,l,islope) =  tsud
     656            end do
     657          enddo
    638658
    639659
    … …  
    689709                patm = patm + ps(i,j)*aire(i,j)
    690710                airetot= airetot + aire(i,j)
    691                 pcap = pcap + aire(i,j)*qsurf(ig,igcm_co2)*g
     711                DO islope=1,nslope
     712                  pcap = pcap + aire(i,j)*qsurf(ig,igcm_co2,islope)*g
     713     &   *subslope_dist(ig,islope)/cos(def_slope_mean(islope)*pi/180.)
     714                ENDDO
    692715             ENDDO
    693716          ENDDO
    … …  
    776799           DO iq =1, nqtot
    777800             DO ig=1,ngridmx
    778                  qsurf(ig,iq)=0.
     801               DO islope=1,nslope
     802                 qsurf(ig,iq,islope)=0.
     803               ENDDO
    779804             ENDDO
    780805           ENDDO
    … …  
    798823             
    799824             q(1:iip1,1:jjp1,1:llm,iq)=q(1:iip1,1:jjp1,1:llm,iq)*val
    800              qsurf(1:ngridmx,iq)=qsurf(1:ngridmx,iq)*val
     825             qsurf(1:ngridmx,iq,:)=qsurf(1:ngridmx,iq,:)*val
    801826
    802827c       q=x : initialise tracer manually
    … …  
    828853             read(*,*) val
    829854             DO ig=1,ngridmx
    830                  qsurf(ig,iq)=val
     855               DO islope=1,nslope
     856                 qsurf(ig,iq,islope)=val
     857               ENDDO
    831858             ENDDO
    832859
    … …  
    864891               if (ierr.eq.0) then
    865892                 ! initialize tracer values
    866                  qsurf(:,iq)=profile(1)
     893                 do islope=1,nslope
     894                   qsurf(:,iq,islope)=profile(1)
     895                 enddo
    867896                 do l=1,llm
    868897                   q(:,:,l,iq)=profile(l+1)
    … …  
    10131042
    10141043           do ig=1,ngridmx
    1015            qsurf(ig,igcm_h2o_ice)=max(0.,qsurf(ig,igcm_h2o_ice))
     1044             do islope=1,nslope
     1045               qsurf(ig,igcm_h2o_ice,islope)=
     1046     &           max(0.,qsurf(ig,igcm_h2o_ice,islope))
     1047             enddo
    10161048           end do
    10171049
    … …  
    10221054     &      =q(1:iip1,1:jjp1,1:llm,igcm_h2o_ice)* DoverH
    10231055
    1024          qsurf(1:ngridmx,igcm_hdo_ice)
    1025      &      =qsurf(1:ngridmx,igcm_h2o_ice)*DoverH
     1056         qsurf(1:ngridmx,igcm_hdo_ice,:)
     1057     &      =qsurf(1:ngridmx,igcm_h2o_ice,:)*DoverH
    10261058
    10271059
    … …  
    12491281          write(*,*)'also set negative values of surf ice to 0'
    12501282           do ig=1,ngridmx
    1251               qsurf(ig,igcm_h2o_ice)=min(val,qsurf(ig,igcm_h2o_ice))
    1252               qsurf(ig,igcm_h2o_ice)=max(0.,qsurf(ig,igcm_h2o_ice))
     1283             do islope=1,nslope
     1284              qsurf(ig,igcm_h2o_ice,islope)=
     1285     &         min(val,qsurf(ig,igcm_h2o_ice,islope))
     1286              qsurf(ig,igcm_h2o_ice,islope)=
     1287     &         max(0.,qsurf(ig,igcm_h2o_ice,islope))
     1288             enddo
    12531289           end do
    12541290
    … …  
    12611297              write(*,*) 'OK: remove surface ice for |lat|<45'
    12621298              if (abs(rlatu(j)*180./pi).lt.45.) then
    1263                    qsurf(ig,igcm_h2o_ice)=0.
     1299                do islope=1,nslope
     1300                   qsurf(ig,igcm_h2o_ice,islope)=0.
     1301                enddo
    12641302              end if
    12651303           end do
    … …  
    12731311              if(ig.eq.1) j=1
    12741312              if (rlatu(j)*180./pi.gt.80.) then
    1275                    qsurf(ig,igcm_h2o_ice)=1.e5
    1276                    write(*,*) 'ig=',ig,'    H2O ice mass (kg/m2)= ',
    1277      &              qsurf(ig,igcm_h2o_ice)
     1313                do islope=1,nslope
     1314                   qsurf(ig,igcm_h2o_ice,islope)=1.e5
     1315                   write(*,*) 'ig=',ig,', islope=', islope,
     1316     &              '    H2O ice mass (kg/m2)= ',
     1317     &              qsurf(ig,igcm_h2o_ice,islope)
    12781318                   write(*,*)'     ==> Ice mesh South boundary (deg)= ',
    12791319     &              rlatv(j)*180./pi
     1320                enddo
    12801321                end if
    12811322           enddo
    … …  
    12881329               if(ig.eq.1) j=1
    12891330               if (rlatu(j)*180./pi.lt.-80.) then
    1290                    qsurf(ig,igcm_h2o_ice)=1.e5
    1291                    write(*,*) 'ig=',ig,'   H2O ice mass (kg/m2)= ',
    1292      &              qsurf(ig,igcm_h2o_ice)
     1331                do islope=1,nslope
     1332                   qsurf(ig,igcm_h2o_ice,islope)=1.e5
     1333                   write(*,*) 'ig=',ig,', islope=', islope,
     1334     &             '   H2O ice mass (kg/m2)= ',
     1335     &              qsurf(ig,igcm_h2o_ice,islope)
    12931336                   write(*,*)'     ==> Ice mesh North boundary (deg)= ',
    12941337     &              rlatv(j-1)*180./pi
     1338                enddo
    12951339               end if
    12961340           enddo
    … …  
    13071351
    13081352          do ig=1, ngridmx
    1309             tsurf(ig) = Tiso
     1353            do islope=1,nslope
     1354            tsurf(ig,islope) = Tiso
     1355            enddo
    13101356          end do
    13111357          do l=2,nsoilmx
    13121358            do ig=1, ngridmx
    1313               tsoil(ig,l) = Tiso
     1359              do islope=1,nslope
     1360              tsoil(ig,l,islope) = Tiso
     1361              enddo
    13141362            end do
    13151363          end do
    … …  
    13231371        else if (trim(modif) .eq. 'co2ice=0') then
    13241372           do ig=1,ngridmx
    1325               qsurf(ig,igcm_co2)=0
    1326               emis(ig)=emis(ngridmx/2)
     1373              do islope=1,nslope
     1374              qsurf(ig,igcm_co2,islope)=0
     1375              emis(ig,islope)=emis(ngridmx/2,islope)
     1376              enddo
    13271377           end do
    13281378       
    … …  
    16701720        ! recast ith() into ithfi()
    16711721        CALL gr_dyn_fi(nsoilmx,iip1,jjp1,ngridmx,ith,ithfi)
     1722
     1723        else if (trim(modif) .eq. 'nslope') then
     1724          write(*,*) 'set a new number of subgrid scale slope'
     1725          write(*,*) 'Current value=', nslope
     1726          write(*,*) 'Enter value for nslope (ex: 1,5,7)?'
     1727          ierr=1
     1728          do while (ierr.ne.0)
     1729            read(*,*,iostat=ierr) nslope_new
     1730          enddo
     1731
     1732       write(*,*) 'This can take some time...'
     1733       write(*,*) 'You can go grab a coffee and relax a bit'
     1734
     1735      if(nslope.eq.nslope_new) then
     1736        write(*,*) 'The number of subslope you entered is the same'
     1737        write(*,*) 'as the number written in startfi.nc. '
     1738        write(*,*) 'Nothing will be done'
     1739      else
     1740
     1741      nslope_old=nslope
     1742      nslope=nslope_new
     1743
     1744      call end_comslope_h
     1745      call ini_comslope_h(ngridmx,nslope_new)
     1746
     1747      allocate(default_def_slope(nslope_new+1))
     1748      !Sub-grid scale  subslopes
     1749      if (nslope_new.eq.7) then
     1750      default_def_slope(1) = -43.
     1751      default_def_slope(2) = -19.
     1752      default_def_slope(3) = -9.
     1753      default_def_slope(4) = -3.
     1754      default_def_slope(5) = 3.
     1755      default_def_slope(6) = 9.
     1756      default_def_slope(7) = 19.
     1757      default_def_slope(8) = 43.
     1758      elseif (nslope_new.eq.5) then
     1759      default_def_slope(1) = -43.
     1760      default_def_slope(2) = -9.
     1761      default_def_slope(3) = -3.
     1762      default_def_slope(4) = 3.
     1763      default_def_slope(5) = 9.
     1764      default_def_slope(6) = 43.
     1765      elseif (nslope_new.eq.1) then
     1766      default_def_slope(1) = -50.
     1767      default_def_slope(2) = 50.
     1768      endif
     1769
     1770      do islope=1,nslope_new+1
     1771       def_slope(islope) = default_def_slope(islope)
     1772      enddo
     1773
     1774      do islope=1,nslope_new
     1775       def_slope_mean(islope) =(def_slope(islope)+def_slope(islope+1))/2.
     1776      enddo
     1777
     1778       iflat = 1
     1779       DO islope=2,nslope_new
     1780         IF(abs(def_slope_mean(islope)).lt.
     1781     &      abs(def_slope_mean(iflat)))THEN
     1782           iflat = islope
     1783         ENDIF
     1784       ENDDO
     1785
     1786       call subslope_mola(ngridmx,nslope_new,def_slope,subslope_dist)
     1787
     1788! Surfdat related stuff
    16721789       
     1790        allocate(tsurf_old_slope(ngridmx,nslope_old)) 
     1791        allocate(qsurf_old_slope(ngridmx,nqtot,nslope_old)) 
     1792        allocate(emis_old_slope(ngridmx,nslope_old))   
     1793        allocate(watercap_old_slope(ngridmx,nslope_old))
     1794
     1795        tsurf_old_slope(:,:)=tsurf(:,:)
     1796        qsurf_old_slope(:,:,:)=qsurf(:,:,:)
     1797        emis_old_slope(:,:)=emis(:,:)
     1798        watercap_old_slope(:,:)=watercap(:,:)
     1799
     1800        call end_surfdat_h_slope_var
     1801        call ini_surfdat_h_slope_var(ngridmx,nqtot,nslope_new)
     1802
     1803! Comsoil related stuff (tsoil)
     1804
     1805        allocate(tsoil_old_slope(ngridmx,nsoilmx,nslope_old))
     1806
     1807        tsoil_old_slope(:,:,:)=tsoil(:,:,:)
     1808
     1809        call end_comsoil_h_slope_var
     1810        call ini_comsoil_h_slope_var(ngridmx,nslope_new)
     1811
     1812! Dimradmars related stuff (albedo)
     1813
     1814        allocate(albedo_old_slope(ngridmx,2,nslope_old))
     1815        albedo_old_slope(:,:,:)=albedo(:,:,:)
     1816
     1817        call end_dimradmars_mod_slope_var
     1818        call ini_dimradmars_mod_slope_var(ngridmx,nslope_new)
     1819
     1820        if(nslope_old.eq.1 .and. nslope_new.gt.1) then
     1821          do islope=1,nslope_new
     1822             tsurf(:,islope)=tsurf_old_slope(:,1)
     1823             qsurf(:,:,islope)=qsurf_old_slope(:,:,1)
     1824             emis(:,islope)=emis_old_slope(:,1)
     1825             watercap(:,islope)=watercap_old_slope(:,1)
     1826             tsoil(:,:,islope)=tsoil_old_slope(:,:,1)
     1827             albedo(:,:,islope)=albedo_old_slope(:,:,1)
     1828          enddo
     1829        elseif(nslope_new.eq.1) then
     1830             tsurf(:,1)=tsurf_old_slope(:,iflat)
     1831             qsurf(:,:,1)=qsurf_old_slope(:,:,iflat)
     1832             emis(:,1)=emis_old_slope(:,iflat)
     1833             watercap(:,1)=watercap_old_slope(:,iflat)
     1834             tsoil(:,:,1)=tsoil_old_slope(:,:,iflat)
     1835             albedo(:,:,1)=albedo_old_slope(:,:,iflat)
     1836        elseif(nslope_old.eq.5 .and. nslope_new.eq.7) then
     1837          do islope=1,nslope_new
     1838             tsurf(:,islope)=tsurf_old_slope(:,iflat)
     1839             qsurf(:,:,islope)=qsurf_old_slope(:,:,iflat)
     1840             emis(:,islope)=emis_old_slope(:,iflat)
     1841             watercap(:,islope)=watercap_old_slope(:,iflat)
     1842             tsoil(:,:,islope)=tsoil_old_slope(:,:,iflat)
     1843             albedo(:,:,islope)=albedo_old_slope(:,:,iflat)
     1844          enddo
     1845        elseif(nslope_old.eq.7 .and. nslope_new.eq.5) then
     1846          do islope=1,nslope_new
     1847             tsurf(:,islope)=tsurf_old_slope(:,iflat)
     1848             qsurf(:,:,islope)=qsurf_old_slope(:,:,iflat)
     1849             emis(:,islope)=emis_old_slope(:,iflat)
     1850             watercap(:,islope)=watercap_old_slope(:,iflat)
     1851             tsoil(:,:,islope)=tsoil_old_slope(:,:,iflat)
     1852             albedo(:,:,islope)=albedo_old_slope(:,:,iflat)
     1853          enddo
     1854        else
     1855          write(*,*)' Problem choice of nslope'
     1856          write(*,*)' Value not taken into account'
     1857          CALL ABORT
     1858        endif
     1859
     1860        endif !nslope=nslope_new
     1861
    16731862        else
    16741863          write(*,*) '       Unknown (misspelled?) option!!!'
    … …  
    17861975     &              nsoilmx,ngridmx,llm,
    17871976     &              nqtot,dtphys,real(day_ini),0.0,cell_area,
    1788      &              albfi,ithfi)
     1977     &              albfi,ithfi,def_slope,
     1978     &              subslope_dist)
    17891979      call physdem1("restartfi.nc",nsoilmx,ngridmx,llm,nqtot,
    17901980     &              dtphys,hour_ini,
Note: See TracChangeset for help on using the changeset viewer.