Changeset 3483


Ignore:
Timestamp:
Oct 24, 2024, 1:10:51 PM (6 weeks ago)
Author:
afalco
Message:

Pluto PCM: Added thermal inertia input when startphy_file == false.
Renamed thermal inertia variables to therm_inertia.
AF

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

Legend:

Unmodified
Added
Removed

  • trunk/LMDZ.PLUTO/libf/phypluto/phyetat0_mod.F90

    r3275 r3483  
    88                     ngrid,nlayer,fichnom,tab0,Lmodif,nsoil,nq, &
    99                     day_ini,time,tsurf,tsoil, &
    10                      emis,q2,qsurf)
     10                     emis,q2,qsurf,therm_inertia)
    1111                    !  ,cloudfrac,totcloudfrac,hice, &
    1212                    !  rnat,pctsrf_sic)
    … …  
    1818  use write_field_phy, only: Writefield_phy
    1919!!
     20  use comsoil_h, only: nsoilmx
    2021  use tabfi_mod, only: tabfi
    2122  USE tracer_h, ONLY: noms
    … …  
    4647  real,intent(out) :: tsoil(ngrid,nsoil) ! soil temperature
    4748  real,intent(out) :: emis(ngrid) ! surface emissivity
    48   real,intent(out) :: q2(ngrid,nlayer+1) ! 
     49  real,intent(out) :: q2(ngrid,nlayer+1) !
    4950  real,intent(out) :: qsurf(ngrid,nq) ! tracers on surface
     51  real,intent(out) :: therm_inertia(ngrid,nsoilmx) ! thermal inertia
    5052! real n2ice(ngrid) ! n2 ice cover
    5153
    … …  
    6264      INTEGER nid, nvarid
    6365      INTEGER ierr, i, nsrf
    64 !      integer isoil 
     66!      integer isoil
    6567!      INTEGER length
    6668!      PARAMETER (length=100)
    … …  
    6971      CHARACTER*1 yes
    7072!
    71       REAL p_rad,p_omeg,p_g,p_cpp,p_mugaz,p_daysec
     73      REAL p_rad,p_omeg,p_g,p_cpp,p_mugaz,p_daysec,TI
    7274      INTEGER nqold
    7375
    … …  
    7678      integer :: count
    7779      character(len=30) :: txt ! to store some text
    78      
     80
    7981      INTEGER :: indextime=1 ! index of selected time, default value=1
    8082      logical :: found
    81      
     83
    8284      character(len=8) :: modname="phyetat0"
    8385
    … …  
    265267    call soil_settings(nid_start,ngrid,nsoil,tsurf,tsoil,indextime)
    266268  endif
     269else
     270  call getin_p("therm_inertia",TI)
     271  therm_inertia(:,:) = TI
     272    !AF24
    267273endif ! of if (startphy_file)
    268274

  • trunk/LMDZ.PLUTO/libf/phypluto/physiq_mod.F90

    r3455 r3483  
    2828      use comdiurn_h, only: coslat, sinlat, coslon, sinlon
    2929      use comsaison_h, only: mu0, fract, dist_star, declin, right_ascen
    30       use comsoil_h, only: nsoilmx, layer, mlayer, inertiedat
     30      use comsoil_h, only: nsoilmx, layer, mlayer, inertiedat, volcapa
    3131      use geometry_mod, only: latitude, longitude, cell_area
    3232      USE comgeomfi_h, only: totarea, totarea_planet
    … …  
    263263      REAL,save  :: tcond1p4Pa
    264264      DATA tcond1p4Pa/38/
    265       real :: tidat(ngrid,nsoilmx)    ! thermal inertia soil
    266265
    267266
    … …  
    523522#endif
    524523
     524      ! local variables for skin depth check
     525      real :: inertia_min,inertia_max
     526      real :: diurnal_skin ! diurnal skin depth (m)
     527      real :: annual_skin ! anuual skin depth (m)
     528
     529
    525530      ! flags to trigger extra sanity checks
    526531      logical, save ::  check_physics_inputs=.false.
    … …  
    585590         call phyetat0(startphy_file,                                 &
    586591                       ngrid,nlayer,"startfi.nc",0,0,nsoilmx,nq,      &
    587                        day_ini,time_phys,tsurf,tsoil,emis,q2,qsurf)
     592                       day_ini,time_phys,tsurf,tsoil,emis,q2,qsurf,inertiedat)
    588593
    589594#else
    … …  
    9991004      end if ! if ngrid ne 1
    10001005
    1001       call surfprop(ngrid,nq,fract,qsurf,tsurf,tidat,   &
    1002       capcal,adjust,dist_star,albedo,emis,flusurfold,ptimestep,   &
    1003       zls)
     1006      call surfprop(ngrid,nq,fract,qsurf,tsurf,        &
     1007      capcal,adjust,dist_star,flusurfold,ptimestep,zls,&
     1008      albedo,emis,inertiedat)
    10041009      ! do i=2,ngrid
    10051010      !    albedo(i,:) = albedo(1,:)
    10061011      ! enddo
     1012
     1013      if (callsoil) then
     1014         ! AF24 Originally in soil.F, but inertiedat is modified by surfprop
     1015         ! Additional checks: is the vertical discretization sufficient
     1016         ! to resolve diurnal and annual waves?
     1017         do ig=1,ngrid
     1018            ! extreme inertia for this column
     1019            inertia_min=minval(inertiedat(ig,:))
     1020            inertia_max=maxval(inertiedat(ig,:))
     1021            ! diurnal and annual skin depth
     1022            diurnal_skin=(inertia_min/volcapa)*sqrt(daysec/pi)
     1023            annual_skin=(inertia_max/volcapa)*sqrt(year_day*daysec/pi)
     1024            if (0.5*diurnal_skin<layer(1)) then
     1025            ! one should have the fist layer be at least half of diurnal skin depth
     1026            write(*,*) "soil Error: grid point ig=",ig
     1027            write(*,*) "            longitude=",longitude(ig)*(180./pi)
     1028            write(*,*) "             latitude=",latitude(ig)*(180./pi)
     1029            write(*,*) "  first soil layer depth ",layer(1)
     1030            write(*,*) "  not small enough for a diurnal skin depth of ", &
     1031                        diurnal_skin
     1032            write(*,*) " change soil layer distribution (comsoil_h.F90)"
     1033            call abort_physic("soil","change soil layer distribution (comsoil_h.F90)",1)
     1034            endif
     1035            if (2.*annual_skin>layer(nsoilmx)) then
     1036            ! one should have the full soil be at least twice the diurnal skin depth
     1037            write(*,*) "soil Error: grid point ig=",ig
     1038            write(*,*) "            longitude=",longitude(ig)*(180./pi)
     1039            write(*,*) "             latitude=",latitude(ig)*(180./pi)
     1040            write(*,*) "  total soil layer depth ",layer(nsoilmx)
     1041            write(*,*) "  not large enough for an annual skin depth of ", &
     1042                        annual_skin
     1043            write(*,*) " change soil layer distribution (comsoil_h.F90)"
     1044            call abort_physic("soil","change soil layer distribution (comsoil_h.F90)",1)
     1045            endif
     1046         enddo ! of do ig=1,ngrid
     1047
     1048      end if ! callsoil
    10071049
    10081050      if (callrad) then
    … …  
    17951837      ! tidat_out(:,:)=0.
    17961838      ! DO l=1,min(nlayermx,nsoilmx)
    1797       !    tidat_out(:,l)=tidat(:,l)
     1839      !    tidat_out(:,l)=inertiedat(:,l)
    17981840      ! ENDDO
    17991841
    … …  
    22172259            !      ptimestep,pdaypal, &
    22182260            !      ztime_fin,tsurf,tsoil,emis,q2,qsurfpal, &
    2219             !      cell_area,albedodat,tidat,zmea,zstd,zsig, &
     2261            !      cell_area,albedodat,inertiedat,zmea,zstd,zsig, &
    22202262            !      zgam,zthe,oblipal,eccpal,tpalnew,adjustnew,phisfipal,  &
    22212263            !      peri_daypal)

  • trunk/LMDZ.PLUTO/libf/phypluto/soil.F90

    r3482 r3483  
    1       subroutine soil(ngrid,nsoil,firstcall,lastcall,
    2      &          therm_i,
    3      &          timestep,tsurf,tsoil,
    4      &          capcal,fluxgrd)
     1      subroutine soil(ngrid,nsoil,firstcall,lastcall,   &
     2               therm_i, &
     3               timestep,tsurf,tsoil,    &
     4               capcal,fluxgrd)
    55
    66      use comsoil_h, only: layer, mlayer, volcapa, inertiedat
    … …  
    2121!-----------------------------------------------------------------------
    2222
    23 c-----------------------------------------------------------------------
     23!-----------------------------------------------------------------------
    2424!  arguments
    2525!  ---------
    … …  
    4949! local variables:
    5050      integer ig,ik
    51       real :: inertia_min,inertia_max
    52       real :: diurnal_skin ! diurnal skin depth (m)
    53       real :: annual_skin ! anuual skin depth (m)
    5451
    5552! 0. Initialisations and preprocessing step
    … …  
    7673      do ig=1,ngrid
    7774        do ik=1,nsoil-1
    78       thermdiff(ig,ik)=((layer(ik)-mlayer(ik-1))*mthermdiff(ig,ik)
    79      &                +(mlayer(ik)-layer(ik))*mthermdiff(ig,ik-1))
    80      &                    /(mlayer(ik)-mlayer(ik-1))
     75      thermdiff(ig,ik)=((layer(ik)-mlayer(ik-1))*mthermdiff(ig,ik)  &
     76                     +(mlayer(ik)-layer(ik))*mthermdiff(ig,ik-1))   &
     77                         /(mlayer(ik)-mlayer(ik-1))
    8178!       write(*,*),'soil: ik: ',ik,' thermdiff:',thermdiff(ig,ik)
    8279        enddo
    … …  
    9188        ! q_{k+1/2}
    9289        do ik=1,nsoil-1
    93           coefq(ik)=volcapa*(layer(ik+1)-layer(ik))
    94      &                 /timestep
     90          coefq(ik)=volcapa*(layer(ik+1)-layer(ik)) &
     91                      /timestep
    9592        enddo
    9693
    … …  
    10299
    103100        ! alph_{N-1}
    104         alph(ig,nsoil-1)=coefd(ig,nsoil-1)/
    105      &                  (coefq(nsoil-1)+coefd(ig,nsoil-1))
     101        alph(ig,nsoil-1)=coefd(ig,nsoil-1)/ &
     102                       (coefq(nsoil-1)+coefd(ig,nsoil-1))
    106103        ! alph_k
    107104        do ik=nsoil-2,1,-1
    108           alph(ig,ik)=coefd(ig,ik)/(coefq(ik)+coefd(ig,ik+1)*
    109      &                              (1.-alph(ig,ik+1))+coefd(ig,ik))
     105          alph(ig,ik)=coefd(ig,ik)/(coefq(ik)+coefd(ig,ik+1)*   &
     106                                   (1.-alph(ig,ik+1))+coefd(ig,ik))
    110107        enddo
    111108
    112109        ! capcal
    113110! Cstar
    114         capcal(ig)=volcapa*layer(1)+
    115      &              (thermdiff(ig,1)/(mlayer(1)-mlayer(0)))*
    116      &              (timestep*(1.-alph(ig,1)))
     111        capcal(ig)=volcapa*layer(1)+    &
     112                   (thermdiff(ig,1)/(mlayer(1)-mlayer(0)))* &
     113                   (timestep*(1.-alph(ig,1)))
    117114! Cs
    118         capcal(ig)=capcal(ig)/(1.+mu*(1.0-alph(ig,1))*
    119      &                         thermdiff(ig,1)/mthermdiff(ig,0))
     115        capcal(ig)=capcal(ig)/(1.+mu*(1.0-alph(ig,1))*  &
     116                              thermdiff(ig,1)/mthermdiff(ig,0))
    120117      !write(*,*)'soil: ig=',ig,' capcal(ig)=',capcal(ig)
    121       enddo ! of do ig=1,ngrid
    122 
    123       ! Additional checks: is the vertical discretization sufficient
    124       ! to resolve diurnal and annual waves?
    125       do ig=1,ngrid
    126         ! extreme inertia for this column
    127         inertia_min=minval(inertiedat(ig,:))
    128         inertia_max=maxval(inertiedat(ig,:))
    129         ! diurnal and annual skin depth
    130         diurnal_skin=(inertia_min/volcapa)*sqrt(daysec/pi)
    131         annual_skin=(inertia_max/volcapa)*sqrt(year_day*daysec/pi)
    132         if (0.5*diurnal_skin<layer(1)) then
    133         ! one should have the fist layer be at least half of diurnal skin depth
    134           write(*,*) "soil Error: grid point ig=",ig
    135           write(*,*) "            longitude=",longitude(ig)*(180./pi)
    136           write(*,*) "             latitude=",latitude(ig)*(180./pi)
    137           write(*,*) "  first soil layer depth ",layer(1)
    138           write(*,*) "  not small enough for a diurnal skin depth of ",
    139      &                diurnal_skin
    140           write(*,*) " change soil layer distribution (comsoil_h.F90)"
    141           stop
    142         endif
    143         if (2.*annual_skin>layer(nsoil)) then
    144         ! one should have the full soil be at least twice the diurnal skin depth
    145           write(*,*) "soil Error: grid point ig=",ig
    146           write(*,*) "            longitude=",longitude(ig)*(180./pi)
    147           write(*,*) "             latitude=",latitude(ig)*(180./pi)
    148           write(*,*) "  total soil layer depth ",layer(nsoil)
    149           write(*,*) "  not large enough for an annual skin depth of ",
    150      &                annual_skin
    151           write(*,*) " change soil layer distribution (comsoil_h.F90)"
    152           stop
    153         endif
    154118      enddo ! of do ig=1,ngrid
    155119
    … …  
    160124! First layer:
    161125      do ig=1,ngrid
    162         tsoil(ig,1)=(tsurf(ig)+mu*beta(ig,1)*
    163      &                         thermdiff(ig,1)/mthermdiff(ig,0))/
    164      &              (1.+mu*(1.0-alph(ig,1))*
    165      &               thermdiff(ig,1)/mthermdiff(ig,0))
     126        tsoil(ig,1)=(tsurf(ig)+mu*beta(ig,1)*   &
     127                              thermdiff(ig,1)/mthermdiff(ig,0))/    &
     128                   (1.+mu*(1.0-alph(ig,1))* &
     129                    thermdiff(ig,1)/mthermdiff(ig,0))
    166130      enddo
    167131! Other layers:
    … …  
    177141! Bottom layer, beta_{N-1}
    178142      do ig=1,ngrid
    179         beta(ig,nsoil-1)=coefq(nsoil-1)*tsoil(ig,nsoil)
    180      &                   /(coefq(nsoil-1)+coefd(ig,nsoil-1))
     143        beta(ig,nsoil-1)=coefq(nsoil-1)*tsoil(ig,nsoil) &
     144                        /(coefq(nsoil-1)+coefd(ig,nsoil-1))
    181145      enddo
    182146! Other layers
    183147      do ik=nsoil-2,1,-1
    184148        do ig=1,ngrid
    185           beta(ig,ik)=(coefq(ik)*tsoil(ig,ik+1)+
    186      &                 coefd(ig,ik+1)*beta(ig,ik+1))/
    187      &                 (coefq(ik)+coefd(ig,ik+1)*(1.0-alph(ig,ik+1))
    188      &                  +coefd(ig,ik))
     149          beta(ig,ik)=(coefq(ik)*tsoil(ig,ik+1)+    &
     150                      coefd(ig,ik+1)*beta(ig,ik+1))/    &
     151                      (coefq(ik)+coefd(ig,ik+1)*(1.0-alph(ig,ik+1)) &
     152                       +coefd(ig,ik))
    189153        enddo
    190154      enddo
    … …  
    206170!         print*,'tsoil=',tsoil(ig,1)
    207171
    208         fluxgrd(ig)=(thermdiff(ig,1)/(mlayer(1)-mlayer(0)))*
    209      &              (beta(ig,1)+(alph(ig,1)-1.0)*tsoil(ig,1))
     172        fluxgrd(ig)=(thermdiff(ig,1)/(mlayer(1)-mlayer(0)))*    &
     173                   (beta(ig,1)+(alph(ig,1)-1.0)*tsoil(ig,1))
    210174
    211175!        mu=mlayer(ig,0)/(mlayer(ig,1)-mlayer(ig,0))
    … …  
    213177!     &                         thermdiff(ig,1)/mthermdiff(ig,0))
    214178! Fs
    215         fluxgrd(ig)=fluxgrd(ig)+(capcal(ig)/timestep)*
    216      &              (tsoil(ig,1)*(1.+mu*(1.0-alph(ig,1))*
    217      &                         thermdiff(ig,1)/mthermdiff(ig,0))
    218      &               -tsurf(ig)-mu*beta(ig,1)*
    219      &                          thermdiff(ig,1)/mthermdiff(ig,0))
     179        fluxgrd(ig)=fluxgrd(ig)+(capcal(ig)/timestep)*  &
     180                   (tsoil(ig,1)*(1.+mu*(1.0-alph(ig,1))*    &
     181                              thermdiff(ig,1)/mthermdiff(ig,0)) &
     182                    -tsurf(ig)-mu*beta(ig,1)*   &
     183                               thermdiff(ig,1)/mthermdiff(ig,0))
    220184      enddo
    221185

  • trunk/LMDZ.PLUTO/libf/phypluto/surfprop.F90

    r3377 r3483  
    1       subroutine surfprop(ngrid,nq,pfra,qsurf,tsurface,     &
    2                      tidat,capcal,adjust,dist,albedo,emis,fluold, &
    3                      ptimestep,zls)
     1      subroutine surfprop(ngrid,nq,pfra,qsurf,tsurface,       &
     2                     capcal,adjust,dist,fluold,ptimestep,zls, &
     3                     albedo,emis,therm_inertia)
    44
    55    !   use comgeomfi_h, only:
    … …  
    4343!     albedo(ngrid)
    4444!     emis(ngrid)
     45!     therm_inertia(ngrid,nsoilmx)
    4546!
    4647!     Both
    … …  
    5657!     Arguments
    5758
     59
    5860      INTEGER ngrid, nq
     61      REAL,INTENT(IN) :: pfra(ngrid)
    5962      REAL,INTENT(IN) :: qsurf(ngrid,nq)
     63      REAL,INTENT(IN) :: tsurface(ngrid)
    6064      REAL,INTENT(IN) :: fluold(ngrid,nq)
    6165      REAL,INTENT(IN) :: ptimestep
    6266      REAL,INTENT(IN) :: zls
    63       REAL,INTENT(IN) :: tsurface(ngrid)
    6467      REAL,INTENT(IN) :: capcal(ngrid)
    6568      REAL,INTENT(IN) :: adjust
    … …  
    6770      REAL,INTENT(OUT) :: albedo(ngrid,L_NSPECTV)
    6871      REAL,INTENT(OUT) :: emis(ngrid)
    69       REAL,INTENT(OUT) :: tidat(ngrid,nsoilmx)
    70       REAL,INTENT(IN) :: pfra(ngrid)
     72      REAL,INTENT(OUT) :: therm_inertia(ngrid,nsoilmx) ! therm_inertia = inertiedat
    7173!-----------------------------------------------------------------------
    7274!     Local variables
    … …  
    448450           do isoil=0,nsoilmx-1
    449451              if (mlayer(isoil).le.emin) then ! diurnal TI
    450                    tidat(ig,isoil+1)=tid
     452                   therm_inertia(ig,isoil+1)=tid
    451453              else if (isoil.gt.0.and.(mlayer(isoil).gt.emin).and.(mlayer(isoil-1).lt.emin)) then ! still diurnal TI
    452                    tidat(ig,isoil+1)=tid
     454                   therm_inertia(ig,isoil+1)=tid
    453455              else if ((mlayer(isoil).gt.emin).and.(mlayer(isoil).le.emax)) then ! TI N2
    454                    tidat(ig,isoil+1)=ITN2
     456                   therm_inertia(ig,isoil+1)=ITN2
    455457              else if ((mlayer(isoil).gt.emin).and.(mlayer(isoil).le.emax2)) then
    456                    tidat(ig,isoil+1)=ITCH4
     458                   therm_inertia(ig,isoil+1)=ITCH4
    457459              else
    458                    tidat(ig,isoil+1)=ITH2O
     460                   therm_inertia(ig,isoil+1)=ITH2O
    459461              endif
    460462
    … …  
    465467
    466468        DO ig=1,ngrid
    467            tidat(ig,:)=inertiedat(ig,:)
     469           therm_inertia(ig,:)=inertiedat(ig,:)
    468470        ENDDO
    469471
Note: See TracChangeset for help on using the changeset viewer.