Changeset 4040 for trunk/LMDZ.PLUTO/libf


Ignore:
Timestamp:
Feb 4, 2026, 10:38:48 PM (12 hours ago)
Author:
tbertrand
Message:

PLUTO PCM:
Implementing paleo mode (to be continued)
TB

Location:
trunk/LMDZ.PLUTO/libf/phypluto
Files:
4 edited

Legend:

Unmodified
Added
Removed
  • trunk/LMDZ.PLUTO/libf/phypluto/initracer.F90

    r3951 r4040  
    187187      igcm_eddy5e8=0
    188188      write(*,*) 'initracer: noms() ', noms
    189       rho_n2=1000        ! n2 ice
     189
     190      ! Densities :
     191      rho_n2=1000.           ! n2 ice
    190192      rho_ch4_ice=520.       ! rho ch4 ice
    191       rho_co_ice=520.       ! rho ch4 ice
    192       rho_haze=800.     ! haze
     193      rho_co_ice=520.        ! rho CO ice
     194      rho_haze=800.          ! haze
    193195
    194196      ! 1. find dust tracers
     
    200202          igcm_n2=iq
    201203          mmol(igcm_n2)=28.
     204          rho_q(iq)=rho_n2
    202205          count=count+1
    203206          write(*,*) 'Tracer ',count,' = n2'
  • trunk/LMDZ.PLUTO/libf/phypluto/physiq_mod.F90

    r4026 r4040  
    2323      use aerosol_mod, only: i_haze, haze_prof
    2424      use surfdat_h, only: phisfi, zmea, zstd, zsig, zgam, zthe, &
    25                            dryness
     25                           dryness, qsurfyear, phisfibed
    2626      use comdiurn_h, only: coslat, sinlat, coslon, sinlon
    2727      use comsaison_h, only: mu0, fract, dist_star, declin, right_ascen
     
    2929      use geometry_mod, only: latitude, longitude, cell_area, &
    3030                          cell_area_for_lonlat_outputs
     31      use control_mod, only: iphysiq
    3132      USE comgeomfi_h, only: totarea, totarea_planet
    3233      USE tracer_h, only: noms, mmol, radius, rho_q, qext, &
     
    252253! ----------------------
    253254      integer,save :: day_ini                                      ! Initial date of the run (sol since Ls=0).
     255      real,save    :: time_phys                                    ! Initial time of day of the run
    254256      integer,save :: icount                                       ! Counter of calls to physiq during the run.
    255257!$OMP THREADPRIVATE(day_ini,icount)
    256258
    257       !Pluto specific
     259      ! Pluto specific
    258260      REAL,save  :: acond,bcond
    259261      REAL,save  :: tcond1p4Pa
     
    262264! Local variables :
    263265! -----------------
    264       !     Tendencies for the paleoclimate mode
    265       REAL qsurfyear(ngrid,nq)  ! kg.m-2 averaged mass of ice lost/gained in the last Pluto year of the run
    266       REAL phisfinew(ngrid)       ! geopotential of the bedrock (= phisfi-qsurf/1000*g)
    267       REAL qsurfpal(ngrid,nq)           ! qsurf after a paleoclimate step : for physdem1 and restartfi
    268       REAL phisfipal(ngrid)               ! geopotential after a paleoclimate step : for physdem1 and restartfi
     266      ! Variables for the paleoclimate mode
     267      REAL qsurfpal(ngrid,nq)        ! qsurf after a paleoclimate step : for physdem1 and restartfi
     268      REAL phisfipal(ngrid)          ! geopotential after a paleoclimate step : for physdem1 and restartfi
    269269      REAL oblipal                   ! change of obliquity
    270270      REAL peri_daypal               ! new periday
     
    274274      REAL pdaypal                   ! new pday = day_ini + step
    275275      REAL zdt_tot                   ! time range corresponding to the flux of qsurfyear
    276       REAL massacc(nq)             ! accumulated mass
    277       REAL masslost(nq)            ! accumulated mass
    278 
    279       REAL globave                   ! globalaverage 2D ps
    280       REAL globaveice(nq)          ! globalaverage 2D ice
    281       REAL globavenewice(nq)          ! globalaverage 2D ice
    282       INTEGER lecttsoil     ! lecture of tsoil from proftsoil
    283       REAL qsurf1(ngrid,nq)      ! saving qsurf to calculate flux over long timescales kg.m-2
    284       REAL flusurf(ngrid,nq)     ! flux cond/sub kg.m-2.s-1
    285       REAL flusurfold(ngrid,nq)  ! old flux cond/sub kg.m-2.s-1
    286       REAL zplev(ngrid,nlayer+1),zplay(ngrid,nlayer)
    287       REAL vecnull(ngrid,nq)      ! null vector used to check conservation of tracer
     276      REAL massacc(nq)               ! accumulated mass
     277      REAL masslost(nq)              ! accumulated mass
     278      REAL qsurf1(ngrid,nq)          ! saving qsurf to calculate flux over long timescales kg.m-2
     279      REAL flusurf(ngrid,nq)         ! flux cond/sub kg.m-2.s-1
     280      REAL flusurfold(ngrid,nq)      ! old flux cond/sub kg.m-2.s-1
     281      REAL globaveice(nq)            ! globalaverage 2D ice
     282      REAL globavenewice(nq)         ! globalaverage 2D ice
    288283
    289284      REAL,SAVE :: ptime0    ! store the first time
    290285      REAL dstep
    291286      REAL,SAVE :: glastep=20   ! step in pluto day to spread glacier
     287
     288      ! Other variables
     289      REAL globave                   ! globalaverage 2D ps
     290      INTEGER lecttsoil              ! lecture of tsoil from proftsoil
     291      REAL zplev(ngrid,nlayer+1),zplay(ngrid,nlayer) ! pressure density levels
     292      REAL vecnull(ngrid,nq)         ! null vector used to check conservation of tracer
    292293
    293294!     Aerosol (dust or ice) extinction optical depth  at reference wavelength
     
    341342      real zdqsed(ngrid,nlayer,nq)    ! Callsedim routine.
    342343      real zdqmr(ngrid,nlayer,nq)     ! Mass_redistribution routine.
    343       REAL,allocatable,save :: zdqchim(:,:,:) ! Calchim_asis routine
    344       REAL,allocatable,save :: zdqschim(:,:)  ! Calchim_asis routine
    345 !$OMP THREADPRIVATE(zdqchim,zdqschim)
    346344
    347345      !! PLUTO variables
     
    438436      real reff(ngrid,nlayer)                       ! Effective dust radius (used if doubleq=T).
    439437      real vmr(ngrid,nlayer)                        ! volume mixing ratio
    440       real time_phys
    441438
    442439      real ISR,ASR,OLR,GND,DYN,GSR,Ts1,Ts2,Ts3,TsS ! for Diagnostic.
     
    654651         endif
    655652
     653!        Initialize paleo variables
     654!        ~~~~~~~~~~~~~~~~~~~~~~~~~~
     655         if (paleo) then
     656           IF (.not. ALLOCATED(qsurfyear)) ALLOCATE(qsurfyear(ngrid,nq))
     657           IF (.not. ALLOCATED(phisfibed)) ALLOCATE(phisfibed(ngrid))
     658           write(*,*) 'Paleo time tpal = ',tpal
     659           qsurfyear(:,:)=0.
     660           DO ig=1,ngrid
     661             phisfibed(ig)=phisfi(ig)-qsurf(ig,igcm_n2)*g/rho_q(igcm_n2) ! topo of bedrock below the ice
     662           ENDDO
     663         endif
     664
    656665!        Initialize correlated-k.
    657666!        ~~~~~~~~~~~~~~~~~~~~~~~~
     
    785794      enddo
    786795
    787       !Altitude of top interface (nlayer+1), using the thicknesss of the level below the top one. LT22
     796      !Altitude of top interface (nlayer+1), using the thickness of the level below the top one. LT22
    788797
    789798      zzlev(1:ngrid,nlayer+1) = 2*zzlev(1:ngrid,nlayer)-zzlev(1:ngrid,nlayer-1)
     
    19221931           IF (paleo) then
    19231932             call spreadglacier_paleo(ngrid,nq,qsurf, &
    1924                                     phisfinew,dstep,tsurf)
     1933                                    phisfibed,dstep,tsurf)
    19251934           else
    19261935             call spreadglacier_simple(ngrid,nq,qsurf,dstep)
     
    22242233
    22252234            ! update new geopotential depending on the ice reservoir
    2226             phisfipal(:)=phisfinew(:)+qsurfpal(:,igcm_n2)*g/1000.
     2235            phisfipal(:)=phisfibed(:)+qsurfpal(:,igcm_n2)*g/rho_q(igcm_n2)
    22272236            !phisfipal(ig)=phisfi(ig)
    22282237
     
    22542263            ! create restartfi
    22552264            if (ngrid.ne.1) then
     2265               ztime_restart = ptime + ptimestep/(iphysiq*daysec)
    22562266               call physdem0pal("restartfi.nc",longitude,latitude,nsoilmx,ngrid,nlayer,nq, &
    2257                          ptimestep,pdaypal,time_phys,cell_area,          &
     2267                         ptimestep,pdaypal,ztime_restart,cell_area,          &
    22582268                         albedo_bareground,zmea,zstd,zsig,zgam,zthe,     &
    22592269                         oblipal,eccpal,tpalnew,adjustnew,phisfipal,peri_daypal) 
    2260  
    2261                  !call physdem1pal("restartfi.nc",long,lati,nsoilmx,nq, &
    2262                !      ptimestep,pdaypal, &
    2263                !      ztime_restart,tsurf,tsoil,emis,q2,qsurfpal, &
    2264                !      cell_area,albedodat,therm_inertia,zmea,zstd,zsig, &
    2265                !      zgam,zthe,oblipal,eccpal,tpalnew,adjustnew,phisfipal,  &
    2266                !      peri_daypal)
    22672270            endif
    2268          else ! 'paleo'
    2269 
    22702271
    22712272         endif ! end of 'paleo'
     
    22862287!              in 'restart'. Between now and the writing of 'restart',
    22872288!              there will have been the itau=itau+1 instruction and
    2288 !              a reset of 'time' (lastacll = .true. when itau+1= itaufin)
     2289!              a reset of 'time' (lastcall = .true. when itau+1= itaufin)
    22892290!              thus we store for time=time+dtvr
    22902291
     
    23052306         endif ! ngrid
    23062307      endif ! is_omp_master
    2307 
     2308      if (lastcall.and.paleo) then
     2309            call physdem1("restartfi.nc",nsoilmx,ngrid,nlayer,nq, &
     2310                        ptimestep,ztime_restart,tsurf,                &
     2311                        tsoil,therm_inertia,emis,albedo,q2,qsurfpal,n2frac)
     2312            if (is_master) write(*,*)'PHYSIQ: writing restartfi at time =',ztime_restart
     2313      endif
    23082314
    23092315
     
    24202426      ! Diagnostics of optical thickness (dtau = dtau_gas + dtau_rayaer + dtau_cont).
    24212427      ! Warning this is exp(-dtau), I let you postproc with -log to have tau and k itself
    2422       ! VI
    2423       call write_output('dtauv_4656nm','Layer optical thickness attenuation in VI band','',int_dtauv(:,nlayer:1:-1,2))  ! 4.656 um (28 VIS Bands)
    2424       call write_output('dtauv_1181nm','Layer optical thickness attenuation in VI band','',int_dtauv(:,nlayer:1:-1,21)) ! 1.181 um (28 VIS Bands)
    2425       call write_output('dtauv_700nm','Layer optical thickness attenuation in VI band','',int_dtauv(:,nlayer:1:-1,24))  ! 0.700 um (28 VIS Bands)
    2426       call write_output('dtauv_185nm','Layer optical thickness attenuation in VI band','',int_dtauv(:,nlayer:1:-1,27))  ! 0.185 um (28 VIS Bands)
    2427       call write_output('dtauv_118nm','Layer optical thickness attenuation in VI band','',int_dtauv(:,nlayer:1:-1,28))  ! 0.118 um (28 VIS Bands)
    2428       ! IR
    2429       call write_output('dtaui_81250nm','Layer optical thickness attenuation in IR band','',int_dtaui(:,nlayer:1:-1,2)) ! 81.250 um (17 IR Bands)
    2430       call write_output('dtaui_3859nm','Layer optical thickness attenuation in IR band','',int_dtaui(:,nlayer:1:-1,16)) ! 3.859 um (17 IR Bands)
     2428      !! VI
     2429      !call write_output('dtauv_4656nm','Layer optical thickness attenuation in VI band','',int_dtauv(:,nlayer:1:-1,2))  ! 4.656 um (28 VIS Bands)
     2430      !call write_output('dtauv_1181nm','Layer optical thickness attenuation in VI band','',int_dtauv(:,nlayer:1:-1,21)) ! 1.181 um (28 VIS Bands)
     2431      !call write_output('dtauv_700nm','Layer optical thickness attenuation in VI band','',int_dtauv(:,nlayer:1:-1,24))  ! 0.700 um (28 VIS Bands)
     2432      !call write_output('dtauv_185nm','Layer optical thickness attenuation in VI band','',int_dtauv(:,nlayer:1:-1,27))  ! 0.185 um (28 VIS Bands)
     2433      !call write_output('dtauv_118nm','Layer optical thickness attenuation in VI band','',int_dtauv(:,nlayer:1:-1,28))  ! 0.118 um (28 VIS Bands)
     2434      !! IR
     2435      !call write_output('dtaui_81250nm','Layer optical thickness attenuation in IR band','',int_dtaui(:,nlayer:1:-1,2)) ! 81.250 um (17 IR Bands)
     2436      !call write_output('dtaui_3859nm','Layer optical thickness attenuation in IR band','',int_dtaui(:,nlayer:1:-1,16)) ! 3.859 um (17 IR Bands)
    24312437     
    24322438      !call write_output('dtaui_25250nm','Layer optical thickness attenuation in IR band','',int_dtaui(:,nlayer:1:-1,4)) ! 25.250 um (25 IR Bands)
     
    24352441      !call write_output('dtaui_15050nm','Layer optical thickness attenuation in IR band','',int_dtaui(:,nlayer:1:-1,10)) ! 15.050 um (25 IR Bands)
    24362442     
    2437       if (callmufi) then
     2443      !if (callmufi) then
    24382444         ! Aerosol optical thickness
    2439          call write_output('dtauv_aers_4656nm','Layer sph. aer. optical thickness attenuation in VI band','',int_dtauv_aer(:,nlayer:1:-1,2,1))
    2440          call write_output('dtauv_aerf_4656nm','Layer fra. aer. optical thickness attenuation in VI band','',int_dtauv_aer(:,nlayer:1:-1,2,2))
    2441          call write_output('dtauv_aers_1181nm','Layer sph. aer. optical thickness attenuation in VI band','',int_dtauv_aer(:,nlayer:1:-1,21,1))
    2442          call write_output('dtauv_aerf_1181nm','Layer fra. aer. optical thickness attenuation in VI band','',int_dtauv_aer(:,nlayer:1:-1,21,2))
    2443          call write_output('dtauv_aers_700nm','Layer sph. aer. optical thickness attenuation in VI band','',int_dtauv_aer(:,nlayer:1:-1,24,1))
    2444          call write_output('dtauv_aerf_700nm','Layer fra. aer. optical thickness attenuation in VI band','',int_dtauv_aer(:,nlayer:1:-1,24,2))
    2445          call write_output('dtauv_aers_185nm','Layer sph. aer. optical thickness attenuation in VI band','',int_dtauv_aer(:,nlayer:1:-1,27,1))
    2446          call write_output('dtauv_aerf_185nm','Layer fra. aer. optical thickness attenuation in VI band','',int_dtauv_aer(:,nlayer:1:-1,27,2))
    2447          call write_output('dtauv_aers_118nm','Layer sph. aer. optical thickness attenuation in VI band','',int_dtauv_aer(:,nlayer:1:-1,28,1))
    2448          call write_output('dtauv_aerf_118nm','Layer fra. aer. optical thickness attenuation in VI band','',int_dtauv_aer(:,nlayer:1:-1,28,2))
    2449          ! Aerosols single scattering albedo
    2450          call write_output('wbarv_aers_4656nm','Layer sph. aer. single scattering albedo in VI band','',int_wbarv_aer(:,nlayer:1:-1,2,1))
    2451          call write_output('wbarv_aerf_4656nm','Layer fra. aer. single scattering albedo in VI band','',int_wbarv_aer(:,nlayer:1:-1,2,2))
    2452          call write_output('wbarv_aers_1181nm','Layer sph. aer. single scattering albedo in VI band','',int_wbarv_aer(:,nlayer:1:-1,21,1))
    2453          call write_output('wbarv_aerf_1181nm','Layer fra. aer. single scattering albedo in VI band','',int_wbarv_aer(:,nlayer:1:-1,21,2))
    2454          call write_output('wbarv_aers_700nm','Layer sph. aer. single scattering albedo in VI band','',int_wbarv_aer(:,nlayer:1:-1,24,1))
    2455          call write_output('wbarv_aerf_700nm','Layer fra. aer. single scattering albedo in VI band','',int_wbarv_aer(:,nlayer:1:-1,24,2))
    2456          call write_output('wbarv_aers_185nm','Layer sph. aer. single scattering albedo in VI band','',int_wbarv_aer(:,nlayer:1:-1,27,1))
    2457          call write_output('wbarv_aerf_185nm','Layer fra. aer. single scattering albedo in VI band','',int_wbarv_aer(:,nlayer:1:-1,27,2))
    2458          call write_output('wbarv_aers_118nm','Layer sph. aer. single scattering albedo in VI band','',int_wbarv_aer(:,nlayer:1:-1,28,1))
    2459          call write_output('wbarv_aerf_118nm','Layer fra. aer. single scattering albedo in VI band','',int_wbarv_aer(:,nlayer:1:-1,28,2))
    2460       endif ! end callmufi
     2445         !call write_output('dtauv_aers_4656nm','Layer sph. aer. optical thickness attenuation in VI band','',int_dtauv_aer(:,nlayer:1:-1,2,1))
     2446         !call write_output('dtauv_aerf_4656nm','Layer fra. aer. optical thickness attenuation in VI band','',int_dtauv_aer(:,nlayer:1:-1,2,2))
     2447         !call write_output('dtauv_aers_1181nm','Layer sph. aer. optical thickness attenuation in VI band','',int_dtauv_aer(:,nlayer:1:-1,21,1))
     2448         !call write_output('dtauv_aerf_1181nm','Layer fra. aer. optical thickness attenuation in VI band','',int_dtauv_aer(:,nlayer:1:-1,21,2))
     2449         !call write_output('dtauv_aers_700nm','Layer sph. aer. optical thickness attenuation in VI band','',int_dtauv_aer(:,nlayer:1:-1,24,1))
     2450         !call write_output('dtauv_aerf_700nm','Layer fra. aer. optical thickness attenuation in VI band','',int_dtauv_aer(:,nlayer:1:-1,24,2))
     2451         !call write_output('dtauv_aers_185nm','Layer sph. aer. optical thickness attenuation in VI band','',int_dtauv_aer(:,nlayer:1:-1,27,1))
     2452         !call write_output('dtauv_aerf_185nm','Layer fra. aer. optical thickness attenuation in VI band','',int_dtauv_aer(:,nlayer:1:-1,27,2))
     2453         !call write_output('dtauv_aers_118nm','Layer sph. aer. optical thickness attenuation in VI band','',int_dtauv_aer(:,nlayer:1:-1,28,1))
     2454         !call write_output('dtauv_aerf_118nm','Layer fra. aer. optical thickness attenuation in VI band','',int_dtauv_aer(:,nlayer:1:-1,28,2))
     2455         !! Aerosols single scattering albedo
     2456         !call write_output('wbarv_aers_4656nm','Layer sph. aer. single scattering albedo in VI band','',int_wbarv_aer(:,nlayer:1:-1,2,1))
     2457         !call write_output('wbarv_aerf_4656nm','Layer fra. aer. single scattering albedo in VI band','',int_wbarv_aer(:,nlayer:1:-1,2,2))
     2458         !call write_output('wbarv_aers_1181nm','Layer sph. aer. single scattering albedo in VI band','',int_wbarv_aer(:,nlayer:1:-1,21,1))
     2459         !call write_output('wbarv_aerf_1181nm','Layer fra. aer. single scattering albedo in VI band','',int_wbarv_aer(:,nlayer:1:-1,21,2))
     2460         !call write_output('wbarv_aers_700nm','Layer sph. aer. single scattering albedo in VI band','',int_wbarv_aer(:,nlayer:1:-1,24,1))
     2461         !call write_output('wbarv_aerf_700nm','Layer fra. aer. single scattering albedo in VI band','',int_wbarv_aer(:,nlayer:1:-1,24,2))
     2462         !call write_output('wbarv_aers_185nm','Layer sph. aer. single scattering albedo in VI band','',int_wbarv_aer(:,nlayer:1:-1,27,1))
     2463         !call write_output('wbarv_aerf_185nm','Layer fra. aer. single scattering albedo in VI band','',int_wbarv_aer(:,nlayer:1:-1,27,2))
     2464         !call write_output('wbarv_aers_118nm','Layer sph. aer. single scattering albedo in VI band','',int_wbarv_aer(:,nlayer:1:-1,28,1))
     2465         !call write_output('wbarv_aerf_118nm','Layer fra. aer. single scattering albedo in VI band','',int_wbarv_aer(:,nlayer:1:-1,28,2))
     2466      !endif ! end callmufi
    24612467
    24622468      if (calllott) then
  • trunk/LMDZ.PLUTO/libf/phypluto/spreadglacier_paleo.F90

    r3585 r4040  
    1919!     ngrid                 Number of vertical columns
    2020!     pqsurf(ngrid,nq)      N2 ice tracer on surface kg/m2
    21 !     phisfi0               = phisfinew the geopotential of the bedrock
     21!     phisfi0               = phisfibed the geopotential of the bedrock
    2222!                             below the N2 ice
    2323!     ptsurf                surface temperature
  • trunk/LMDZ.PLUTO/libf/phypluto/surfdat_h.F90

    r3910 r4040  
    1313       real,allocatable,dimension(:) :: phisfi ! geopotential at ground level
    1414!$OMP THREADPRIVATE(phisfi)
     15       real,allocatable,dimension(:) :: phisfibed ! (paleo) geopotential of the bedrock (= phisfi-qsurf/1000*g)
     16!$OMP THREADPRIVATE(phisfibed)
     17       real,allocatable,dimension(:,:) :: qsurfyear ! (paleo) kg.m-2 averaged mass of ice lost/gained in the last year
     18!$OMP THREADPRIVATE(qsurfyear)
    1519       real,dimension(2) :: emisice ! ice emissivity; 1:Northern hemisphere 2:Southern hemisphere
    1620       real emissiv
Note: See TracChangeset for help on using the changeset viewer.