Ignore:
Timestamp:
Sep 11, 2023, 12:41:50 PM (15 months ago)
Author:
jbclement
Message:

Mars PEM:
Big cleaning of main program pem.F90 (indentation, declarations, comments, simplification of conditions/loops, etc).
JBC

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/LMDZ.COMMON/libf/evolution/pemetat0.F90

    r3019 r3039  
    1    SUBROUTINE pemetat0(filename,ngrid,nsoil_GCM,nsoil_PEM,nslope,timelen,timestep,TI_PEM,tsoil_PEM,ice_table, ice_table_thickness, &
     1SUBROUTINE pemetat0(filename,ngrid,nsoil_GCM,nsoil_PEM,nslope,timelen,timestep,TI_PEM,tsoil_PEM,ice_table, ice_table_thickness, &
    22                      tsurf_ave_yr1,tsurf_ave_yr2,q_co2,q_h2o,ps_inst,tsoil_inst,tend_h2oglaciers,tend_co2glaciers,co2ice,waterice, &
    33                      global_ave_pressure,watersurf_ave,watersoil_ave, m_co2_regolith_phys,deltam_co2_regolith_phys,  &
    44                      m_h2o_regolith_phys,deltam_h2o_regolith_phys, water_reservoir)
    55   
    6    use iostart_PEM, only:  open_startphy, close_startphy, get_field, get_var
    7    use comsoil_h_PEM, only: soil_pem,layer_PEM, mlayer_PEM,fluxgeo,inertiedat_PEM,water_reservoir_nom,depth_breccia,depth_bedrock,index_breccia,index_bedrock
    8    use comsoil_h, only:  volcapa,inertiedat
    9    use adsorption_mod, only : regolith_adsorption,adsorption_pem
    10    USE temps_mod_evol, ONLY: year_PEM
    11    USE ice_table_mod, only: computeice_table_equilibrium
    12    USE constants_marspem_mod,only: alpha_clap_h2o,beta_clap_h2o, &
    13                                    TI_breccia,TI_bedrock
    14    use soil_thermalproperties_mod, only: update_soil_thermalproperties
    15    use tracer_mod, only: mmol,igcm_h2o_vap ! tracer names and molar masses
     6  use iostart_PEM,                only: open_startphy, close_startphy, get_field, get_var
     7  use comsoil_h_PEM,              only: soil_pem, layer_PEM, mlayer_PEM, fluxgeo, inertiedat_PEM, water_reservoir_nom, depth_breccia, depth_bedrock, index_breccia, index_bedrock
     8  use comsoil_h,                  only: volcapa,inertiedat
     9  use adsorption_mod,             only: regolith_adsorption, adsorption_pem
     10  use ice_table_mod,              only: computeice_table_equilibrium
     11  use constants_marspem_mod,      only: alpha_clap_h2o, beta_clap_h2o, TI_breccia, TI_bedrock
     12  use soil_thermalproperties_mod, only: update_soil_thermalproperties
     13  use tracer_mod,                 only: mmol, igcm_h2o_vap ! tracer names and molar masses
    1614
    1715#ifndef CPP_STD
    18       USE comcstfi_h, only: r, mugaz
     16    use comcstfi_h,   only: r, mugaz
     17    use surfdat_h,    only: watercaptag
    1918#else
    20       USE comcstfi_mod, only: r, mugaz
     19    use comcstfi_mod, only: r, mugaz
    2120#endif
    2221
    23 
    24 #ifndef CPP_STD   
    25    use surfdat_h, only: watercaptag
    26 #endif
    27 
    2822  implicit none
    29  
     23
     24  include "callkeys.h"
     25
    3026  character(len=*), intent(in) :: filename                    ! name of the startfi_PEM.nc
    3127  integer,intent(in) :: ngrid                                 ! # of physical grid points
    … …  
    7369   real :: alph_tmp(ngrid,nsoil_PEM-1)                                ! Intermediate for tsoil computation []
    7470   real :: beta_tmp(ngrid,nsoil_PEM-1)                                ! Intermediate for tsoil computatio []                               
    75    real :: year_PEM_read                                              ! Year of the PEM previous run
    7671   LOGICAL :: startpem_file                                           ! boolean to check if we read the startfile or not
    7772#ifdef CPP_STD   
    … …  
    9792!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    9893
    99 
    10094!0. Check if the start_PEM exist.
    10195
    … …  
    110104   call open_startphy(filename)
    111105
    112    call get_var("Time",year_PEM_read,found)
    113    year_PEM=INT(year_PEM_read)
    114    if(.not.found) then
    115       write(*,*)'PEMetat0: failed loading year_PEM; take default=0'
    116    else
    117       write(*,*)'year_PEM of startpem=', year_PEM
    118    endif
    119 
    120     if(soil_pem) then
     106    if (soil_pem) then
    121107 
    122108!1. Thermal Inertia
    … …  
    161147ENDDO ! islope
    162148
    163 print *,'PEMETAT0: THERMAL INERTIA DONE'
     149write(*,*) 'PEMETAT0: THERMAL INERTIA DONE'
    164150
    165151! b. Special case for inertiedat, inertiedat_PEM
    … …  
    204190
    205191!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    206 
    207192!2. Soil Temperature
    208 
    209193DO islope=1,nslope
    210194  write(num,fmt='(i2.2)') islope
    … …  
    260244        enddo
    261245      enddo
    262 
    263    
    264246ENDDO ! islope
    265247
    266 print *,'PEMETAT0: SOIL TEMP  DONE'
    267 
    268 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    269 
     248write(*,*) 'PEMETAT0: SOIL TEMP  DONE'
     249
     250!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    270251!3. Ice Table
    271 
    272 
    273252  call get_field("ice_table",ice_table,found)
    274253   if(.not.found) then
    … …  
    282261   endif
    283262
    284 print *,'PEMETAT0: ICE TABLE  DONE'
    285 
    286 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    287 
     263write(*,*) 'PEMETAT0: ICE TABLE  DONE'
     264
     265!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    288266!4. CO2 & H2O Adsorption
    289 
    290267 if(adsorption_pem) then
    291268  DO islope=1,nslope
    … …  
    318295       deltam_h2o_regolith_phys(:) = 0. 
    319296    endif
    320 print *,'PEMETAT0: CO2 & H2O adsorption done  '
     297write(*,*) 'PEMETAT0: CO2 & H2O adsorption DONE'
    321298    endif
    322299endif ! soil_pem
    323300
    324 
    325 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    326 
     301!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    327302!. 5 water reservoir
    328 
    329303#ifndef CPP_STD   
    330304   call get_field("water_reservoir",water_reservoir,found)
    … …  
    342316#endif
    343317
    344 
    345318call close_startphy
    346319
    347320else !No startfi, let's build all by hand
    348321
    349     year_PEM=0
    350 
    351     if(soil_pem) then
     322    if (soil_pem) then
    352323
    353324!a) Thermal inertia
    … …  
    407378enddo
    408379
    409 
    410380!!! zone de transition
    411381delta = depth_bedrock
    … …  
    424394 enddo
    425395
    426 print *,'PEMETAT0: TI  DONE'
    427 
    428 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    429 
    430 
     396write(*,*) 'PEMETAT0: TI DONE'
     397
     398!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    431399!b) Soil temperature
    432400    do islope = 1,nslope
    … …  
    458426        enddo
    459427     enddo
    460  
    461       do isoil = nsoil_GCM+1,nsoil_PEM
    462         do ig = 1,ngrid
    463         watersoil_ave(ig,isoil,islope) = exp(beta_clap_h2o/tsoil_PEM(ig,isoil,islope) + alpha_clap_h2o)/tsoil_PEM(ig,isoil,islope)*mmol(igcm_h2o_vap)/(mugaz*r)
     428
     429        do isoil = nsoil_GCM+1,nsoil_PEM
     430          do ig = 1,ngrid
     431            watersoil_ave(ig,isoil,islope) = exp(beta_clap_h2o/tsoil_PEM(ig,isoil,islope) + alpha_clap_h2o)/tsoil_PEM(ig,isoil,islope)*mmol(igcm_h2o_vap)/(mugaz*r)
     432          enddo
    464433        enddo
    465       enddo
    466434enddo !islope
    467 print *,'PEMETAT0: TSOIL DONE  '
    468 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    469 !c) Ice table   
    470      call computeice_table_equilibrium(ngrid,nslope,nsoil_PEM,watercaptag,watersurf_ave,watersoil_ave,TI_PEM(:,1,:),ice_table,ice_table_thickness)
    471      call update_soil_thermalproperties(ngrid,nslope,nsoil_PEM,tend_h2oglaciers,waterice,global_ave_pressure,ice_table,ice_table_thickness,TI_PEM)
    472      do islope = 1,nslope
    473      call soil_pem_ini(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_ave_yr2(:,islope),tsoil_PEM(:,:,islope))
    474      enddo
    475    
    476 
    477 
    478 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    479 
    480 print *,'PEMETAT0: Ice table DONE  '
    481 
    482 
     435write(*,*) 'PEMETAT0: TSOIL DONE'
     436
     437!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     438!c) Ice table
     439       call computeice_table_equilibrium(ngrid,nslope,nsoil_PEM,watercaptag,watersurf_ave,watersoil_ave,TI_PEM(:,1,:),ice_table,ice_table_thickness)
     440       call update_soil_thermalproperties(ngrid,nslope,nsoil_PEM,tend_h2oglaciers,waterice,global_ave_pressure,ice_table,ice_table_thickness,TI_PEM)
     441       do islope = 1,nslope
     442         call soil_pem_ini(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_ave_yr2(:,islope),tsoil_PEM(:,:,islope))
     443       enddo
     444       write(*,*) 'PEMETAT0: Ice table DONE'
    483445
    484446!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    485447!d) Regolith adsorbed
    486 
    487448 if(adsorption_pem) then
    488449    m_co2_regolith_phys(:,:,:) = 0.
    … …  
    496457 endif
    497458
    498 print *,'PEMETAT0: CO2 adsorption done  '
    499 
     459write(*,*) 'PEMETAT0: CO2 adsorption DONE'
    500460endif !soil_pem
    501461
    502 
    503 
    504 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    505    
    506 
    507 
     462!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    508463!. e) water reservoir
    509 
    510464#ifndef CPP_STD   
    511 
    512465      do ig=1,ngrid
    513466        if(watercaptag(ig)) then
    … …  
    517470        endif
    518471      enddo
    519 
    520472#endif
    521473
    … …  
    527479      DO islope = 1,nslope
    528480        DO iloop = 1,nsoil_PEM
    529            if(isnan(tsoil_PEM(ig,iloop,islope))) then
    530          call abort_pem("PEM - pemetat0","NAN detected in Tsoil",1)
    531            endif
     481          if (isnan(tsoil_PEM(ig,iloop,islope))) call abort_pem("PEM - pemetat0","NaN detected in Tsoil",1)
    532482        ENDDO
    533483      ENDDO
    … …  
    535485endif!soil_pem
    536486
    537 
    538 
    539    END SUBROUTINE
     487END SUBROUTINE
Note: See TracChangeset for help on using the changeset viewer.