Ignore:
Timestamp:
Oct 23, 2023, 7:02:06 PM (15 months ago)
Author:
evos
Message:

Mars PCM

We added the to module vdifc the possibily of subsurface intercation, mostly to have the option of buried glaicer that can lose ice and create polar layers.
if no ice is on the surface the ice the atmosphere can directly interact with the subsurface ice.
if ice is on the surface it can interact with the subsurface ice.
For now it is not working with more then one subslope. I will add this feature soon
In order to use it there are two flags that needs to be set to true.
paleoclimate=.true.
lag_layer=.true.
in addition one need to give the following parametrs in startfi or run.def:
d_coef=X for the diffusion coeficent (4e-4 is the defualt)
<h2o_ice_depth=X> if its zero or lower there is no subsurface ice (default -1).
outputs: zdq_ssi, momentarily zdqsdif after the interaction with the SSIzdq_subtimestep,  momentarily flux after the interaction with the SSIzdq_ssi_frost momentarily zdqsdif after the interaction of the frost with the SSI
q1, atmospheric first layer water vapor quatity
qeq, SSI water vapor quantity

Eran Vos

Location:
trunk/LMDZ.MARS/libf/phymars
Files:
8 edited

Legend:

Unmodified
Added
Removed

  • trunk/LMDZ.MARS/libf/phymars/conf_phys.F

    r3094 r3098  
    3737      use aeropacity_mod, only: iddist, topdustref
    3838      USE mod_phys_lmdz_transfert_para, ONLY: bcast
    39       USE paleoclimate_mod,ONLY: paleoclimate,albedo_perenialco2
     39      USE paleoclimate_mod,ONLY: paleoclimate,albedo_perenialco2,
     40     &                           lag_layer
    4041      use microphys_h, only: mteta
    4142
    … …  
    338339! PALEOCLIMATE
    339340
    340          write(*,*)"Is it paleoclimate run?"
     341         write(*,*)"Using lag layer??"
     342         lag_layer=.false.
     343         call getin_p("lag_layer",lag_layer)
     344         write(*,*) " lag_layer = ", lag_layer
     345
     346        write(*,*)"Is it paleoclimate run?"
    341347         paleoclimate=.false. ! default value
    342348         call getin_p("paleoclimate",paleoclimate)
    343349         write(*,*)" paleoclimate = ",paleoclimate
     350
    344351
    345352         write(*,*)"Albedo for perenial CO2 ice?"

  • trunk/LMDZ.MARS/libf/phymars/dyn1d/init_testphys1d_mod.F90

    r3095 r3098  
    477477    call getin("subsurface_ice_depth",ice_depth)
    478478
     479
    479480    z0(1) = z0_default ! default value for roughness
    480481    write(*,*) 'Surface roughness length z0 (m)?'
    … …  
    716717        write(*,*) 'Prescribed atmospheric water vapor profile'
    717718        write(*,*) 'Unless it reaches saturation (maximal value)'
     719    else if (atm_wat_profile .eq. 2) then
     720        write(*,*) 'Prescribed atmospheric water vapor profile'
     721        write(*,*) 'like present day northern midlatitude'
     722        write(*,*) 'Unless it reaches saturation (maximal value)'
    718723    else
    719724        error stop 'Water vapor profile value not correct!'

  • trunk/LMDZ.MARS/libf/phymars/dyn1d/testphys1d.F90

    r3092 r3098  
    179179            endif
    180180        endif
    181     endif
     181    !EV profile
     182!            IF(atm_wat_profile.eq.2) THEN
     183!             DO ilayer=1,16
     184!              q(1,ilayer,igcm_h2o_vap)=min(zqsat(ilayer),6e-5)
     185!             ENDDO! ilayer=1,16
     186!             DO ilayer=17,22
     187!              q(1,ilayer,igcm_h2o_vap)=min(zqsat(ilayer),3e-5)
     188!             ENDDO! ilayer=17,22
     189!             DO ilayer=23,nlayer
     190!              q(1,ilayer,igcm_h2o_vap)=0
     191!             ENDDO! ilayer=23,nlayer
     192!             endif
     193     endif
    182194
    183195    ! Call physics

  • trunk/LMDZ.MARS/libf/phymars/paleoclimate_mod.F90

    r3007 r3098  
    1515!$OMP THREADPRIVATE(paleoclimate)
    1616
    17     real, save, allocatable :: lag_h2o_ice(:,:)  ! Thickness of the lag before H2O ice [m]
     17    real, save, allocatable :: h2o_ice_depth(:,:)  ! Thickness of the lag before H2O ice [m]
    1818    real, save, allocatable :: lag_co2_ice(:,:)  ! Thickness of the lag before CO2 ice [m]
    1919    real, save :: albedo_perenialco2             ! Albedo for perenial co2 ice [1]
    20 !$OMP THREADPRIVATE(lag_h2o_ice,lag_co2_ice,albedo_perenialco2)
     20    real, save, allocatable :: d_coef(:,:)  ! Diffusion coeficent
     21    LOGICAL,SAVE :: lag_layer ! does lag layer is present?
     22!$OMP THREADPRIVATE(h2o_ice_depth,lag_co2_ice,albedo_perenialco2)
    2123
    2224    CONTAINS
    … …  
    2931  integer,intent(in) :: nslope ! number of slope within a mesh
    3032
    31     allocate(lag_h2o_ice(ngrid,nslope))
     33    allocate(h2o_ice_depth(ngrid,nslope))
    3234    allocate(lag_co2_ice(ngrid,nslope))
     35    allocate(d_coef(ngrid,nslope))
    3336  end subroutine ini_paleoclimate_h
    3437
    … …  
    3639
    3740  implicit none
    38     if (allocated(lag_h2o_ice)) deallocate(lag_h2o_ice)
     41    if (allocated(d_coef)) deallocate(d_coef)
     42    if (allocated(h2o_ice_depth)) deallocate(h2o_ice_depth)
    3943    if (allocated(lag_co2_ice)) deallocate(lag_co2_ice)
    4044  end subroutine end_paleoclimate_h

  • trunk/LMDZ.MARS/libf/phymars/phyetat0_mod.F90

    r3037 r3098  
    2727  use comsoil_h, only: flux_geo
    2828  USE comslope_mod, ONLY: nslope, major_slope
    29   USE paleoclimate_mod, ONLY: paleoclimate, lag_h2o_ice, lag_co2_ice
     29  USE paleoclimate_mod, ONLY: paleoclimate, h2o_ice_depth, lag_co2_ice, d_coef
    3030  USE comcstfi_h, only: pi
    3131  use geometry_mod, only: latitude
    … …  
    744744  if (startphy_file) then
    745745! Depth of H2O lag
    746    call get_field("lag_h2o_ice",lag_h2o_ice,found,indextime)
    747    if (.not.found) then
    748      write(*,*) "phyetat0: Failed loading <lag_h2o_ice> : ", &
    749                           "<lag_h2o_ice> is set as -1 (no subsurface ice)"
    750      lag_h2o_ice(:,:) = -1.
    751    endif
     746   call get_field("h2o_ice_depth",h2o_ice_depth,found,indextime)
     747   if (.not.found) then
     748     write(*,*) "phyetat0: Failed loading <h2o_ice_depth> : ", &
     749                          "<h2o_ice_depth> is set as -1 (no subsurface ice)"
     750     h2o_ice_depth(:,:) = -1.
     751   endif
     752
     753! Diffuesion coeficent
     754   call get_field("d_coef",d_coef,found,indextime)
     755   if (.not.found) then
     756     write(*,*) "phyetat0: Failed loading <d_coef> : ", &
     757                          "<d_coef> is set as 4e-4 (defualt value)"
     758     d_coef(:,:) = 4e-4
     759   endif
     760
    752761
    753762!  Depth of CO2 lag
    … …  
    774783   endif ! notfound
    775784  else ! no startfiphyle
    776      lag_h2o_ice(:,:) = -1.
     785     h2o_ice_depth(:,:) = -1.
    777786     lag_co2_ice(:,:) = -1.
    778787     perenial_co2ice(:,:) = 0.
     788     d_coef(:,:)= 4.e-4
    779789     if (abs(latitude(ngrid)-(-pi/2.)).lt.1.e-5) then
    780790       DO islope = 1,nslope
    … …  
    786796  endif !startphy_file
    787797else
    788   lag_h2o_ice(:,:) = -1.
     798  h2o_ice_depth(:,:) = -1.
    789799  lag_co2_ice(:,:) = -1.
    790800  perenial_co2ice(:,:) = 0.
     801  d_coef(:,:)= 4.e-4
    791802endif !paleoclimate
    792803

  • trunk/LMDZ.MARS/libf/phymars/phyredem.F90

    r3051 r3098  
    2525  use time_phylmdz_mod, only: daysec
    2626  use comslope_mod, ONLY: nslope
    27   USE paleoclimate_mod, ONLY: paleoclimate, lag_h2o_ice, lag_co2_ice 
     27  USE paleoclimate_mod, ONLY: paleoclimate, h2o_ice_depth, lag_co2_ice 
    2828  implicit none
    2929 
    … …  
    156156  ! Paleoclimate outputs
    157157  if(paleoclimate) then
    158     call put_field("lag_h2o_ice","Depth of the H2O lags",lag_h2o_ice)
     158    call put_field("h2o_ice_depth","Depth to the fisrt H2O ice",h2o_ice_depth)
    159159    call put_field("lag_co2_ice","Depth of the CO2 lags",lag_co2_ice)
    160160  endif

  • trunk/LMDZ.MARS/libf/phymars/physiq_mod.F

    r3095 r3098  
    15141514     &        zdqdif,zdqsdif,wstar,zcdv,zcdh,hfmax_th,
    15151515     &        zcondicea_co2microp,sensibFlux,
    1516      &        dustliftday,local_time,watercap,dwatercap_dif)
     1516     &        dustliftday,local_time,watercap,dwatercap_dif,
     1517     &        tsoil,nsoilmx)
    15171518
    15181519          DO ig=1,ngrid

  • trunk/LMDZ.MARS/libf/phymars/vdifc_mod.F

    r3059 r3098  
    1313     $                pdqdif,pdqsdif,wstar,zcdv_true,zcdh_true,
    1414     $                hfmax,pcondicea_co2microp,sensibFlux,
    15      $                dustliftday,local_time,watercap, dwatercap_dif)
    16 
     15     $                dustliftday,local_time,watercap, dwatercap_dif,
     16     $                tsoil,nsoil)
    1717      use tracer_mod, only: noms, igcm_dust_mass, igcm_dust_number,
    1818     &                      igcm_dust_submicron, igcm_h2o_vap,
    … …  
    3232     &                      subslope_dist,major_slope,iflat
    3333      use microphys_h, only: To
     34      use paleoclimate_mod, only: d_coef,h2o_ice_depth,lag_layer
     35      use comsoil_h, only: layer, mlayer
    3436
    3537      IMPLICIT NONE
    … …  
    6163c   ----------
    6264
    63       INTEGER,INTENT(IN) :: ngrid,nlay
     65      INTEGER,INTENT(IN) :: ngrid,nlay,nsoil
    6466      REAL,INTENT(IN) :: ptimestep
    6567      REAL,INTENT(IN) :: pplay(ngrid,nlay),pplev(ngrid,nlay+1)
    … …  
    7577      REAL,INTENT(IN) :: pcapcal(ngrid,nslope)
    7678      REAL,INTENT(INOUT) :: pq2(ngrid,nlay+1)
     79      REAL,INTENT(IN) :: tsoil(ngrid,nsoil,nslope)
    7780
    7881c    Argument added for condensation:
    … …  
    102105      REAL :: pt(ngrid,nlay)
    103106 
    104       INTEGER ilev,ig,ilay,nlev,islope
     107      INTEGER ilev,ig,ilay,nlev,islope,ik,lice
    105108
    106109      REAL z4st,zdplanck(ngrid)
    … …  
    118121      REAL zcst1
    119122      REAL zu2(ngrid)
     123      REAL Tice(ngrid,nslope) ! subsurface temperature where ice is located.
     124      REAL qeq(ngrid,nslope) ! saturation water vapor in the subsurface
     125      REAL dist_up(ngrid,nslope) !distance from ice to layer above
     126      REAL dist_down(ngrid,nslope) !distance from ice to layer down
     127      REAL dist_sum(ngrid,nslope) ! sum of distance
     128      REAL zdqsdif_ssi(ngrid,nslope) !SSI flux
     129      REAL zdqsdif_ssi_frost(ngrid,nslope) !SSI-frost interaction
    120130
    121131      EXTERNAL SSUM,SCOPY
    … …  
    155165      REAL rho(ngrid) ! near surface air density
    156166      REAL qsat(ngrid)
     167      REAL qsat2(ngrid,nslope)
     168      REAL resist(ngrid,nslope) !subsurface ice flux reduction coef
    157169
    158170      REAL hdoflux(ngrid,nslope)       ! value of vapour flux of HDO
    … …  
    272284      zdqsdif(1:ngrid)=0
    273285      dwatercap_dif(1:ngrid,1:nslope)=0
    274 
     286!      h2o_ice_depth(1:ngrid,1:nslope)=5
    275287c    ** calcul de rho*dz et dt*rho/dz=dt*rho**2 g/dp
    276288c       avec rho=p/RT=p/ (R Theta) (p/ps)**kappa
    … …  
    914926c      Utilisation d'un sous pas de temps afin
    915927c      de decrire le flux de chaleur latente
    916 
     928         
    917929        do iq=1,nq 
    918930          if ((water).and.(iq.eq.igcm_h2o_vap)) then
    … …  
    920932          DO islope = 1,nslope
    921933           DO ig=1,ngrid
     934
     935             
     936
    922937             zqsurf(ig)=pqsurf(ig,igcm_h2o_ice,islope)/
    923938     &             cos(pi*def_slope_mean(islope)/180.)
    … …  
    943958            zq_tmp_vap(ig,:,:) =zq(ig,:,:)
    944959c           Debut du sous pas de temps
    945 
    946960            DO tsub=1,nsubtimestep(ig)
    947961
    … …  
    980994               if ((-zdqsdif(ig)*subtimestep)
    981995     &            .gt.(zqsurf(ig))) then
     996               
     997   
     998                 !EV subsurface ice
     999                IF(h2o_ice_depth(ig,islope) .gt. 0 .and. lag_layer)
     1000     &           then
     1001                zdqsdif(ig)=
     1002     &                        -zqsurf(ig)/subtimestep
     1003                zqsurf(ig)=0
     1004
     1005               DO ik=0,nsoil-2 ! go through all the layers to find the ice locations
     1006                IF((mlayer(ik).le.h2o_ice_depth(ig,islope)).and.
     1007     &           (mlayer(ik+1).gt.h2o_ice_depth(ig,islope))) THEN
     1008                 lice = ik+1
     1009                EXIT
     1010                ENDIF
     1011               ENDDO !of subsurface loop
     1012                IF (lice .gt. 1) then !calculate the distance from the layers
     1013                 dist_up(ig,islope)=(h2o_ice_depth(ig,islope)
     1014     &                              -mlayer(lice-1))
     1015                 dist_down(ig,islope)=(mlayer(lice)
     1016     &                                -h2o_ice_depth(ig,islope))
     1017                 dist_sum(ig,islope)=dist_up(ig,islope)
     1018     &                               +dist_down(ig,islope)
     1019                 Tice(ig,islope)=(dist_up(ig,islope) ! Linear interp to calculate the temp
     1020     &                          *tsoil(ig,lice-1,islope)
     1021     &           /dist_sum(ig,islope))+
     1022     &           (dist_down(ig,islope)*tsoil(ig,lice,islope)
     1023     &           /dist_sum(ig,islope))
     1024                        ELSE
     1025                                Tice(ig,islope)=tsoil(ig,1,islope)
     1026                        ENDIF
     1027              call watersat(1,Tice(ig,1),pplev(ig,1)
     1028     &                     ,qsat2(ig,1))
     1029              qeq(ig,1)=(ztsrf(ig)/Tice(ig,1))
     1030     &                        *qsat2(ig,1)
     1031              resist(ig,1)=(1+(h2o_ice_depth(ig,islope)*zcdh(ig)
     1032     &                           /d_coef(ig,islope)))
     1033              !write(*,*)'R=',resist(ig,islope)
     1034              !write(*,*)'zice=',h2o_ice_depth(ig,islope)
     1035              !write(*,*)'D=',d_coef(ig,islope)
     1036              !write(*,*)'zcdh=',zcdh(ig)
     1037              zb(ig,1)=zb(ig,1)/resist(ig,1) ! change zb to account subsurface ice
     1038              !vdifc algorithem !!!!!needs to change  reseist io to the mean
     1039              !beacuse of the slopes!!!!!
     1040              z1(ig)=1./(za(ig,nlay)+zb(ig,nlay))
     1041             zc(ig,nlay)=za(ig,nlay)*zq_tmp_vap(ig,nlay,iq)*z1(ig)
     1042             zd(ig,nlay)=zb(ig,nlay)*z1(ig)
     1043             DO ilay=nlay-1,2,-1
     1044                 z1(ig)=1./(za(ig,ilay)+zb(ig,ilay)+
     1045     $            zb(ig,ilay+1)*(1.-zd(ig,ilay+1)))
     1046                 zc(ig,ilay)=(za(ig,ilay)*zq_tmp_vap(ig,ilay,iq)+
     1047     $            zb(ig,ilay+1)*zc(ig,ilay+1))*z1(ig)
     1048                 zd(ig,ilay)=zb(ig,ilay)*z1(ig)
     1049             ENDDO
     1050             z1(ig)=1./(za(ig,1)+zb(ig,1)+
     1051     $          zb(ig,2)*(1.-zd(ig,2)))
     1052             zc(ig,1)=(za(ig,1)*zq_tmp_vap(ig,1,iq)+
     1053     $        zb(ig,2)*zc(ig,2)) * z1(ig)
     1054             zd(ig,1)=zb(ig,1)*z1(ig)
     1055             zq1temp(ig)=zc(ig,1)+ zd(ig,1)*qsat(ig)
     1056             zdqsdif_ssi(ig,1)=rho(ig)*dryness(ig)*zcdv(ig)   
     1057     &                      *(zq1temp(ig)-qeq(ig,islope))
     1058                      call write_output('zdq_ssi',
     1059     &                '','',zdqsdif_ssi(ig,1))
     1060                       call write_output('qeq',
     1061     &                '','',qeq(ig,1))
     1062                       call write_output('q1',
     1063     &                '','',zq1temp(ig))
     1064             !write(*,*)'qeq=',qeq(ig,1)
     1065             !write(*,*)'q1=',zq1temp(ig)
     1066             !write(*,*)'zdqsdif_ssi=',zdqsdif_ssi(ig,1)
     1067             !I should some all the interactions with the SSI accoridng
     1068             !to the statistics and the evaluate
     1069             
     1070             !write(*,*) "zdq_ssi*t", zdqsdif_ssi(ig)*subtimestep
     1071             if (zdqsdif_ssi(ig,1)<0) then
     1072               !zdqsdif(ig)=zdqsdif_ssi(ig,1)-(zqsurf(ig)/subtimestep)
     1073               zdqsdif(ig)=zdqsdif(ig)+zdqsdif_ssi(ig,1)
     1074                call write_output('zdq_zdqssi',
     1075     &                '','',zdqsdif(ig)+zdqsdif_ssi(ig,1))
     1076             endif
     1077              ELSE
    9821078c             pdqsdif > 0 : ice condensing
    9831079c             pdqsdif < 0 : ice subliming
    984 c             write(*,*) "subliming more than available frost:  qsurf!"
     1080             write(*,*) "subliming more than available frost:  qsurf!"
    9851081                  zdqsdif(ig)=
    9861082     &                        -zqsurf(ig)/subtimestep
    … …  
    9911087     $           (-zdqsdif(ig)) *subtimestep) *z1(ig)
    9921088                 zq1temp(ig)=zc(ig,1)
     1089                 ENDIF !if h2o_ice_depth>0 and lag_layer
    9931090               endif  !if .not.watercaptag(ig)
    9941091             endif ! if sublim more than surface
     1092
     1093              ! subsurface ice <--> frost interaction !EV
     1094             IF(h2o_ice_depth(ig,islope) .gt. 4e-4 .and. lag_layer
     1095     &        .and. zqsurf(ig) .gt. 0) then
     1096               zdqsdif_ssi_frost(ig,1)=(d_coef(ig,1)
     1097     &                         /h2o_ice_depth(ig,1))
     1098     &          *rho(ig)*dryness(ig)*(qsat(ig)-qeq(ig,1))
     1099               !needs to change to the mean of eq
     1100               zdqsdif(ig)=zdqsdif(ig)-zdqsdif_ssi_frost(ig,1)
     1101                !!!!! zdsqdif_ssi_frosst need to be changed to an
     1102                !average
     1103              ELSEIF (h2o_ice_depth(ig,islope) .gt. 4e-4 .and. lag_layer
     1104     &        .and. watercaptag(ig)) then
     1105              zdqsdif_ssi_frost(ig,1)=(d_coef(ig,1)
     1106     &                         /h2o_ice_depth(ig,1))
     1107     &          *rho(ig)*dryness(ig)*(qsat(ig)-qeq(ig,1))
     1108               zdqsdif(ig)=zdqsdif(ig)-zdqsdif_ssi_frost(ig,1)
     1109               !needs to change to the mean of eq
     1110              ENDIF
     1111              call write_output('zdqsdif_ssi_frost',
     1112     &                '','',zdqsdif_ssi_frost(ig,1))
     1113              call write_output('subtimestep',
     1114     &                '','',subtimestep)
     1115             ! ENDDO !subsurface ice subslope
     1116
     1117
    9951118
    9961119c             Starting upward calculations for water :
    9971120c             Actualisation de h2o_vap dans le premier niveau
    9981121             zq_tmp_vap(ig,1,igcm_h2o_vap)=zq1temp(ig)
    999 
    10001122c    Take into account the H2O latent heat impact on the surface temperature
    10011123              if (latentheat_surfwater) then
    … …  
    10091131     &                              *zq_tmp_vap(ig,ilay-1,iq)
    10101132              ENDDO
    1011 
    10121133c             Subtimestep water budget :
    10131134
    … …  
    10161137              zqsurf(ig)= zqsurf(ig)+(
    10171138     &                       zdqsdif(ig))*subtimestep
    1018 
     1139              if (zqsurf(ig)<0) then
     1140                      zqsurf(ig)=0
     1141              endif
    10191142c             Monitoring instantaneous latent heat flux in W.m-2 :
    10201143              zsurf_h2o_lh(ig,islope) = zsurf_h2o_lh(ig,islope)+
    … …  
    10331156c             Fin du sous pas de temps
    10341157            ENDDO ! tsub=1,nsubtimestep
    1035 
    10361158c             Integration of subtimestep temp and water budget :
    10371159c             (btw could also compute the post timestep temp and ice
    … …  
    10441166     &                       cos(pi*def_slope_mean(islope)/180.))
    10451167     &                       /ptimestep
    1046 
    10471168c if subliming more than qsurf(ice) and on watercaptag, water
    10481169c sublimates from watercap reservoir (dwatercap_dif is <0)
    … …  
    10631184           ENDDO ! of DO ig=1,ngrid
    10641185          ENDDO ! islope
    1065 
    10661186c Some grid box averages: interface with the atmosphere
    10671187       DO ig = 1,ngrid
    … …  
    10751195         ENDDO
    10761196       ENDDO
    1077 
    10781197! Recompute values in kg/m^2 slopped
    10791198        DO ig = 1,ngrid
    … …  
    11671286     &                          "Ground ice latent heat flux",
    11681287     &                               "W.m-2",surf_h2o_lh(:,iflat))
    1169 
     1288                       call write_output('zdq_subtimestep',
     1289     &                '','',zdqsdif(:)*subtimestep)
     1290                       call write_output('zdq_end',
     1291     &                '','',zdqsdif(:))
     1292                       call write_output('vdifc_subtimestep',
     1293     &                '','',subtimestep)
    11701294C       Diagnostic output for HDO
    11711295!        if (hdo) then
Note: See TracChangeset for help on using the changeset viewer.