Index: trunk/LMDZ.COMMON/libf/evolution/changelog.txt =================================================================== --- trunk/LMDZ.COMMON/libf/evolution/changelog.txt (revision 3063) +++ trunk/LMDZ.COMMON/libf/evolution/changelog.txt (revision 3065) @@ -77,5 +77,6 @@ == 26/09/2023 == JBC Minor changes concerning the form of the code in the PEM. +Addition of a file "changelog.txt" in LMDZ.COMMON/libf/evolution/ specific to the PEM rather than using the one for Mars. Completed with changesets since 01/07/2023. -== 28/09/2023 == JBC -Addition of a file "changelog.txt" in LMDZ.COMMON/libf/evolution/ specific to the PEM rather than using the one for Mars. Completed with changesets since 01/07/2023. +== 02/10/2023 == JBC +Initialization of the PEM in 1D through the subroutine "init_testphys1d_mod.F90" + Some adaptations of the Mars PCM in 1D + Update of "launch_pem.sh" in deftank. Index: trunk/LMDZ.COMMON/libf/evolution/init_pem1d_mod.F90 =================================================================== --- trunk/LMDZ.COMMON/libf/evolution/init_pem1d_mod.F90 (revision 3063) +++ (revision ) @@ -1,581 +1,0 @@ -module init_pem1d_mod - -implicit none - -contains - -#ifdef CPP_1D - -subroutine init_pem1d(llm_in,nqtot_in,u,v,temp,q,psurf,time,ap_out,bp_out) - -!-------------- init_pem1d ------------- -! -! This function is a copy of the test_phys_1d.F90 initialisation part. -! It is meant to initialize all the variables in modules and for the ouput -! that can't be initialzed via the startfi file. -! -! Basically it reads run.def and save variables values -! -!-------------- ------------- - - USE ioipsl_getincom, only: getin - use comcstfi_h, only: pi, rad, omeg, g, mugaz, rcp, r, cpp - use time_phylmdz_mod, only: daysec, dtphys, day_step, ecritphy, iphysiq - use planete_h, only: year_day, periheli, aphelie, peri_day, obliquit, emin_turb, lmixmin - use surfdat_h, only: albedodat, z0_default, emissiv, emisice,& - albedice, iceradius, dtemisice, z0,& - zmea, zstd, zsig, zgam, zthe, phisfi,& - watercaptag, watercap, hmons, summit, base,& - tsurf, emis,qsurf - use infotrac, only: nqtot, tname, nqperes,nqfils - use regular_lonlat_mod, only: init_regular_lonlat - use mod_grid_phy_lmdz, only : regular_lonlat - USE comvert_mod, ONLY: ap,bp,aps,bps,pa,preff,sig,& - presnivs,pseudoalt,scaleheight - use dimradmars_mod, only: tauvis,totcloudfrac - use mod_interface_dyn_phys, only: init_interface_dyn_phys - use geometry_mod, only: init_geometry - use dimphy, only : init_dimphy - USE phys_state_var_init_mod, ONLY: phys_state_var_init - use comgeomfi_h, only: sinlat, ini_fillgeom - use slope_mod, only: theta_sl, psi_sl - use comslope_mod, only: def_slope,subslope_dist,def_slope_mean - use dimradmars_mod, only: tauvis,totcloudfrac - use dust_param_mod, only: tauscaling - use tracer_mod, only: igcm_h2o_vap,igcm_h2o_ice,igcm_co2,noms - USE logic_mod, ONLY: hybrid - USE vertical_layers_mod, ONLY: init_vertical_layers - use comsoil_h, only: volcapa, layer, mlayer, inertiedat, & - inertiesoil,nsoilmx,flux_geo - USE read_profile_mod, ONLY: read_profile - use inichim_newstart_mod, only: inichim_newstart - use physics_distribution_mod, only: init_physics_distribution - use iostart, only: open_startphy,get_var, close_startphy -#ifdef CPP_XIOS - use mod_const_mpi, only: COMM_LMDZ -#endif - include "dimensions.h" - integer, parameter :: ngrid = 1 !(2+(jjm-1)*iim - 1/jjm) - integer, parameter :: nlayer = llm - include "callkeys.h" - include "netcdf.inc" - -!-------------- INPUT VARIABLES ------------- - - integer, intent(in) :: llm_in,nqtot_in - -!-------------- OUTPUT VARIABLES ------------- - - real, intent(out) :: u(nlayer), v(nlayer) ! zonal, meridional wind - real, intent(out) :: temp(nlayer) ! temperature at the middle of the layers - real, intent(out) :: q(2,llm_in,nqtot_in) ! tracer mixing ratio (e.g. kg/kg) - real, intent(out) :: psurf(2) - real, intent(out) :: time ! time (0=1!" - stop - endif - endif - ! allocate arrays: - allocate(tname(nq)) - allocate(qsurf(1,1,nq)) - allocate(tnom_transp(nq)) - - ! read tracer names from file traceur.def - do iq = 1,nq - read(90,'(80a)',iostat = ierr) line ! store the line from traceur.def - if (ierr /= 0) then - write(*,*) 'init_pem1d: error reading tracer names...' - stop - endif - ! if format is tnom_0, tnom_transp (isotopes) - read(line,*,iostat = ierr0) tname(iq),tnom_transp(iq) - if (ierr0 /= 0) then - read(line,*) tname(iq) - tnom_transp(iq)='air' - endif - enddo - close(90) - - ! Isotopes: as in the 3D case we have to determine father/son relations for isotopes and carrying fluid - ALLOCATE(nqfils(nqtot)) - nqperes=0 - nqfils(:)=0 - DO iq=1,nqtot - if (tnom_transp(iq) == 'air') then - ! ceci est un traceur père - WRITE(*,*) 'Le traceur',iq,', appele ',trim(tname(iq)),', est un pere' - nqperes=nqperes+1 - else !if (tnom_transp(iq) == 'air') then - ! ceci est un fils. Qui est son père? - WRITE(*,*) 'Le traceur',iq,', appele ',trim(tname(iq)),', est un fils' - continu=.true. - ipere=1 - do while (continu) - if (tnom_transp(iq) .eq. tname(ipere)) then - ! Son père est ipere - WRITE(*,*) 'Le traceur',iq,'appele ',trim(tname(iq)),' est le fils de ',ipere,'appele ',trim(tname(ipere)) - nqfils(ipere)=nqfils(ipere)+1 - continu=.false. - else !if (tnom_transp(iq) == tnom_0(ipere)) then - ipere=ipere+1 - if (ipere.gt.nqtot) then - WRITE(*,*) 'Le traceur',iq,'appele ',trim(tname(iq)),', est orpelin.' - CALL abort_gcm('infotrac_init','Un traceur est orphelin',1) - endif !if (ipere.gt.nqtot) then - endif !if (tnom_transp(iq) == tnom_0(ipere)) then - enddo !do while (continu) - endif !if (tnom_transp(iq) == 'air') then - enddo !DO iq=1,nqtot - WRITE(*,*) 'nqperes=',nqperes - WRITE(*,*) 'nqfils=',nqfils - - ! Initialize tracers here: - write(*,*) "init_pem1d: initializing tracers" - do iq = 1,nq - open(3,file = 'start1D.txt',status = "old",action = "read") - read(3,*) header, qsurf(1,1,iq), (q(1,ilayer,iq), ilayer = 1,nlayer) - if (tname(iq) /= trim(header)) then - write(*,*) 'Tracer names not compatible for initialization with "start1D.txt"!' - stop - endif - enddo - q(2,:,:) = q(1,:,:) - -#ifdef CPP_XIOS - call init_physics_distribution(regular_lonlat,4,1,1,1,nlayer,COMM_LMDZ) -#else - call init_physics_distribution(regular_lonlat,4,1,1,1,nlayer,1) -#endif - -! Discretization (Definition of grid and time steps) -! -------------- - nlevel=nlayer+1 - nsoil=nsoilmx - - day_step=48 ! default value for day_step - write(*,*)'Number of time steps per sol ?' - call getin("day_step",day_step) - write(*,*) " day_step = ",day_step - - ecritphy=day_step ! default value for ecritphy, output every time step - - ndt=10 ! default value for ndt - write(*,*)'Number of sols to run ?' - call getin("ndt",ndt) - write(*,*) " ndt = ",ndt - - dayn=day0+ndt - ndt=ndt*day_step - dtphys=daysec/day_step - -! Imposed surface pressure -! ------------------------------------ - psurf(1) = 610. ! Default value for psurf - write(*,*) 'Surface pressure (Pa)?' - if (.not. therestart1D) then - call getin("psurf",psurf) - else - read(3,*) header, psurf(1) - endif - write(*,*) " psurf = ",psurf(1) - psurf(2) = psurf(1) - -! Reference pressures - pa=20. ! transition pressure (for hybrid coord.) - preff=610. ! reference surface pressure - -! Aerosol properties -! -------------------------------- - tauvis=0.2 ! default value for tauvis (dust opacity) - write(*,'("Reference dust opacity at ",f4.0," Pa ?")')odpref - call getin("tauvis",tauvis) - write(*,*) " tauvis = ",tauvis - -! latitude/longitude -! ------------------ - latitude(1)=0 ! default value for latitude - write(*,*)'latitude (in degrees) ?' - call getin("latitude",latitude(1)) - write(*,*) " latitude = ",latitude - latitude=latitude*pi/180.E+0 - longitude=0.E+0 - longitude=longitude*pi/180.E+0 - -! some initializations (some of which have already been -! done above!) and loads parameters set in callphys.def -! and allocates some arrays -! Mars possible matter with dtphys in input and include!!! -! Initializations below should mimick what is done in iniphysiq for 3D GCM - call init_interface_dyn_phys - call init_regular_lonlat(1,1,longitude,latitude,& - (/0.,0./),(/0.,0./)) - call init_geometry(1,longitude,latitude,& - (/0.,0.,0.,0./),(/0.,0.,0.,0./),& - cell_area) -! Ehouarn: init_vertial_layers called later (because disvert not called yet) -! call init_vertical_layers(nlayer,preff,scaleheight, -! & ap,bp,aps,bps,presnivs,pseudoalt) - call init_dimphy(1,nlayer) ! Initialize dimphy module - call phys_state_var_init(1,llm,nq,tname,& - day0,dayn,time,& - daysec,dtphys,& - rad,g,r,cpp,& - nqperes,nqfils)! MVals: variables isotopes - call ini_fillgeom(1,latitude,longitude,(/1.0/)) - call conf_phys(1,llm,nq) - - ! in 1D model physics are called every time step - ! ovverride iphysiq value that has been set by conf_phys - if (iphysiq/=1) then - write(*,*) "init_pem1d: setting iphysiq=1" - iphysiq=1 - endif - -! Initialize local slope parameters (only matters if "callslope" -! is .true. in callphys.def) - ! slope inclination angle (deg) 0: horizontal, 90: vertical - theta_sl(1)=0.0 ! default: no inclination - call getin("slope_inclination",theta_sl(1)) - ! slope orientation (deg) - ! 0 == Northward, 90 == Eastward, 180 == Southward, 270 == Westward - psi_sl(1)=0.0 ! default value - call getin("slope_orientation",psi_sl(1)) - - ! sub-slopes parameters (assuming no sub-slopes distribution for now). - def_slope(1)=-90 ! minimum slope angle - def_slope(2)=90 ! maximum slope angle - subslope_dist(1,1)=1 ! fraction of subslopes in mesh -! -! for the gravity wave scheme -! --------------------------------- - zmea(1)=0.E+0 - zstd(1)=0.E+0 - zsig(1)=0.E+0 - zgam(1)=0.E+0 - zthe(1)=0.E+0 -! -! for the slope wind scheme -! --------------------------------- -! - hmons(1)=0.E+0 - write(*,*)'hmons is initialized to ',hmons(1) - summit(1)=0.E+0 - write(*,*)'summit is initialized to ',summit(1) - base(1)=0.E+0 -! -! Default values initializing the coefficients calculated later -! --------------------------------- - tauscaling(1)=1. ! calculated in aeropacity_mod.F - totcloudfrac(1)=1. ! calculated in watercloud_mod.F - -! Specific initializations for "physiq" -! ------------------------------------- -! surface geopotential is not used (or useful) since in 1D -! everything is controled by surface pressure - phisfi(1)=0.E+0 - -! Initialization to take into account prescribed winds -! ------------------------------------------------------ - ptif=2.E+0*omeg*sinlat(1) - -! geostrophic wind - gru=10. ! default value for gru - write(*,*)'zonal eastward component of the geostrophic wind (m/s) ?' - call getin("u",gru) - write(*,*) " u = ",gru - grv=0. !default value for grv - write(*,*)'meridional northward component of the geostrophic',& - ' wind (m/s) ?' - call getin("v",grv) - write(*,*) " v = ",grv - -! Initialize winds for first time step - read(3,*) header, (u(ilayer), ilayer = 1,nlayer) - read(3,*) header, (v(ilayer), ilayer = 1,nlayer) - w = 0. ! default: no vertical wind - -! Initialize turbulente kinetic energy - DO ilevel=1,nlevel - q2(ilevel)=0.E+0 - ENDDO - -! CO2 ice on the surface -! ------------------- - ! get the index of co2 tracer (not known at this stage) - igcm_co2=0 - do iq=1,nq - if (trim(tname(iq))=="co2") then - igcm_co2=iq - endif - enddo - if (igcm_co2==0) then - write(*,*) "init_pem1d error, missing co2 tracer!" - stop - endif - -! Compute pressures and altitudes of atmospheric levels -! ---------------------------------------------------------------- -! Vertical Coordinates -! """""""""""""""""""" - hybrid=.true. - write(*,*)'Hybrid coordinates ?' - call getin("hybrid",hybrid) - write(*,*) " hybrid = ", hybrid - - CALL disvert_noterre - ! now that disvert has been called, initialize module vertical_layers_mod - call init_vertical_layers(nlayer,preff,scaleheight,& - ap,bp,aps,bps,presnivs,pseudoalt) - - DO ilevel=1,nlevel - plev(ilevel)=ap(ilevel)+psurf(1)*bp(ilevel) - ENDDO - - DO ilayer=1,nlayer - play(ilayer)=aps(ilayer)+psurf(1)*bps(ilayer) - ENDDO - - DO ilayer=1,nlayer - zlay(ilayer)=-200.E+0 *r*log(play(ilayer)/plev(1))/g - ENDDO - - -! Initialize temperature profile -! -------------------------------------- - pks=psurf(1)**rcp - -! altitude in km in profile: divide zlay by 1000 - tmp1(0)=0.E+0 - DO ilayer=1,nlayer - tmp1(ilayer)=zlay(ilayer)/1000.E+0 - ENDDO - - call profile(nlayer+1,tmp1,tmp2) - - read(3,*) header, tsurf, (temp(ilayer), ilayer = 1,nlayer) - close(3) - -! Initialize soil properties and temperature -! ------------------------------------------ - volcapa=1.e6 ! volumetric heat capacity - - flux_geo_tmp=0. - call getin("flux_geo",flux_geo_tmp) - flux_geo(:,:) = flux_geo_tmp - -! Initialize depths -! ----------------- - do isoil=0,nsoil-1 - mlayer(isoil)=2.e-4*(2.**(isoil-0.5)) ! mid-layer depth - enddo - do isoil=1,nsoil - layer(isoil)=2.e-4*(2.**(isoil-1)) ! layer depth - enddo - -! Initialize traceurs -! --------------------------- - if (photochem.or.callthermos) then - write(*,*) 'Initializing chemical species' - ! flagthermo=0: initialize over all atmospheric layers - flagthermo=0 - ! check if "h2o_vap" has already been initialized - ! (it has been if there is a "profile_h2o_vap" file around) - inquire(file="profile_h2o_vap",exist=present) - if (present) then - flagh2o=0 ! 0: do not initialize h2o_vap - else - flagh2o=1 ! 1: initialize h2o_vap in inichim_newstart - endif - - ! hack to accomodate that inichim_newstart assumes that - ! q & psurf arrays are on the dynamics scalar grid - allocate(qdyn(2,1,llm,nq),psdyn(2,1)) - qdyn(1,1,1:llm,1:nq)=q(1,1:llm,1:nq) - psdyn(1:2,1)=psurf(1) - call inichim_newstart(ngrid,nq,qdyn,qsurf,psdyn,flagh2o,flagthermo) - q(1,1:llm,1:nq)=qdyn(1,1,1:llm,1:nq) - endif - -! Check if the surface is a water ice reservoir -! -------------------------------------------------- - watercaptag(1)=.false. ! Default: no water ice reservoir - write(*,*)'Water ice cap on ground ?' - call getin("watercaptag",watercaptag) - write(*,*) " watercaptag = ",watercaptag - -! Check if the atmospheric water profile is specified -! --------------------------------------------------- - ! Adding an option to force atmospheric water values JN - atm_wat_profile = -1. ! Default: free atm wat profile - if (water) then - write(*,*)'Force atmospheric water vapor profile?' - call getin('atm_wat_profile',atm_wat_profile) - write(*,*) 'atm_wat_profile = ', atm_wat_profile - if (abs(atm_wat_profile + 1.) < 1.e-15) then ! if == -1. - write(*,*) 'Free atmospheric water vapor profile' - write(*,*) 'Total water is conserved in the column' - else if (abs(atm_wat_profile) < 1.e-15) then ! if == 0. - write(*,*) 'Dry atmospheric water vapor profile' - else if (0. < atm_wat_profile .and. atm_wat_profile <= 1.) then - write(*,*) 'Prescribed atmospheric water vapor profile' - write(*,*) 'Unless it reaches saturation (maximal value)' - else - write(*,*) 'Water vapor profile value not correct!' - stop - endif - endif - -! Check if the atmospheric water profile relaxation is specified -! -------------------------------------------------------------- - ! Adding an option to relax atmospheric water values JBC - atm_wat_tau = -1. ! Default: no time relaxation - if (water) then - write(*,*) 'Relax atmospheric water vapor profile?' - call getin('atm_wat_tau',atm_wat_tau) - write(*,*) 'atm_wat_tau = ', atm_wat_tau - if (atm_wat_tau < 0.) then - write(*,*) 'Atmospheric water vapor profile is not relaxed.' - else - if (0. <= atm_wat_profile .and. atm_wat_profile <= 1.) then - write(*,*) 'Relaxed atmospheric water vapor profile towards ', atm_wat_profile - write(*,*) 'Unless it reaches saturation (maximal value)' - else - write(*,*) 'Reference atmospheric water vapor profile not known!' - write(*,*) 'Please, specify atm_wat_profile' - stop - endif - endif - endif - - ap_out = ap - bp_out = bp - -end subroutine init_pem1d - -#endif - -end module init_pem1d_mod Index: trunk/LMDZ.COMMON/libf/evolution/pem.F90 =================================================================== --- trunk/LMDZ.COMMON/libf/evolution/pem.F90 (revision 3063) +++ trunk/LMDZ.COMMON/libf/evolution/pem.F90 (revision 3065) @@ -96,5 +96,6 @@ use physics_distribution_mod, only: init_physics_distribution use mod_grid_phy_lmdz, only: regular_lonlat - use init_pem1d_mod, only: init_pem1d + use init_testphys1d_mod, only: init_testphys1d + use comvert_mod, only: ap, bp #endif @@ -111,12 +112,12 @@ ! Same variable names as in the GCM -integer :: ngrid ! Number of physical grid points -integer :: nlayer ! Number of vertical layer -integer :: nq ! Number of tracer -integer :: day_ini ! First day of the simulation -real :: pday ! Physical day -real :: time_phys ! Same as GCM -real :: ptimestep ! Same as GCM -real :: ztime_fin ! Same as GCM +integer, parameter :: nlayer = llm ! Number of vertical layer +integer :: ngrid ! Number of physical grid points +integer :: nq ! Number of tracer +integer :: day_ini ! First day of the simulation +real :: pday ! Physical day +real :: time_phys ! Same as GCM +real :: ptimestep ! Same as GCM +real :: ztime_fin ! Same as GCM ! Variables to read start.nc @@ -124,18 +125,19 @@ ! Dynamic variables -real :: vcov(ip1jm,llm), ucov(ip1jmp1,llm) ! vents covariants -real :: teta(ip1jmp1,llm) ! temperature potentielle -real, dimension(:,:,:), allocatable :: q ! champs advectes -real :: ps(ip1jmp1) ! pression au sol -real, dimension(:), allocatable :: ps_start_GCM !(ngrid) pression au sol -real, dimension(:,:), allocatable :: ps_timeseries !(ngrid x timelen) ! pression au sol instantannées -real :: masse(ip1jmp1,llm) ! masse d'air -real :: phis(ip1jmp1) ! geopotentiel au sol +real, dimension(ip1jm,llm) :: vcov ! vents covariants +real, dimension(ip1jmp1,llm) :: ucov ! vents covariants +real, dimension(ip1jmp1,llm) :: teta ! temperature potentielle +real, dimension(:,:,:), allocatable :: q ! champs advectes +real, dimension(ip1jmp1) :: ps ! pression au sol +real, dimension(:), allocatable :: ps_start_GCM ! (ngrid) pression au sol +real, dimension(:,:), allocatable :: ps_timeseries ! (ngrid x timelen) ! pression au sol instantannées +real, dimension(ip1jmp1,llm) :: masse ! masse d'air +real, dimension(ip1jmp1) :: phis ! geopotentiel au sol real :: time_0 ! Variables to read starfi.nc -character (len = *), parameter :: FILE_NAME = "startfi_evol.nc" !Name of the file used for initialsing the PEM +character (len = *), parameter :: FILE_NAME = "startfi_evol.nc" ! Name of the file used for initialsing the PEM character*2 str2 -integer :: ncid, varid,status ! Variable for handling opening of files +integer :: ncid, varid, status ! Variable for handling opening of files integer :: phydimid, subdimid, nlayerdimid, nqdimid ! Variable ID for Netcdf files integer :: lonvarid, latvarid, areavarid, sdvarid ! Variable ID for Netcdf files @@ -143,75 +145,74 @@ ! Variables to read starfi.nc and write restartfi.nc -real, dimension(:), allocatable :: longitude ! Longitude read in FILE_NAME and written in restartfi -real, dimension(:), allocatable :: latitude ! Latitude read in FILE_NAME and written in restartfi -real, dimension(:), allocatable :: ap ! Coefficient ap read in FILE_NAME_start and written in restart -real, dimension(:), allocatable :: bp ! Coefficient bp read in FILE_NAME_start and written in restart -real, dimension(:), allocatable :: cell_area ! Cell_area read in FILE_NAME and written in restartfi +real, dimension(:), allocatable :: longitude ! Longitude read in FILE_NAME and written in restartfi +real, dimension(:), allocatable :: latitude ! Latitude read in FILE_NAME and written in restartfi +real, dimension(:), allocatable :: cell_area ! Cell_area read in FILE_NAME and written in restartfi real :: Total_surface ! Total surface of the planet ! Variables for h2o_ice evolution -real :: ini_surf_h2o ! Initial surface of sublimating h2o ice -real :: ini_surf_co2 ! Initial surface of sublimating co2 ice -real :: global_ave_press_GCM ! constant: global average pressure retrieved in the GCM [Pa] -real :: global_ave_press_old ! constant: Global average pressure of initial/previous time step -real :: global_ave_press_new ! constant: Global average pressure of current time step -real, dimension(:,:), allocatable :: zplev_new ! Physical x Atmospheric field : mass of the atmospheric layers in the pem at current time step [kg/m^2] -real, dimension(:,:), allocatable :: zplev_gcm ! same but retrieved from the gcm [kg/m^2] -real, dimension(:,:,:), allocatable :: zplev_new_timeseries ! Physical x Atmospheric x Time: same as zplev_new, but in times series [kg/m ^2] -real, dimension(:,:,:), allocatable :: zplev_old_timeseries ! same but with the time series, for oldest time step -logical :: STOPPING_water ! Logical : is the criterion (% of change in the surface of sublimating water ice) reached? -logical :: STOPPING_1_water ! Logical : is there still water ice to sublimate? -logical :: STOPPING_co2 ! Logical : is the criterion (% of change in the surface of sublimating water ice) reached? -logical :: STOPPING_pressure ! Logical : is the criterion (% of change in the surface pressure) reached? -integer :: criterion_stop ! which criterion is reached ? 1= h2o ice surf, 2 = co2 ice surf, 3 = ps, 4 = orb param -real, save :: A , B, mmean ! Molar mass: intermediate A, B for computations of the mean molar mass of the layer [mol/kg] -real, allocatable :: vmr_co2_gcm(:,:) ! Physics x Times co2 volume mixing ratio retrieve from the gcm [m^3/m^3] -real, allocatable :: vmr_co2_pem_phys(:,:) ! Physics x Times co2 volume mixing ratio used in the PEM -real, allocatable :: q_co2_PEM_phys(:,:) ! Physics x Times co2 mass mixing ratio in the first layer computed in the PEM, first value comes from GCM [kg/kg] -real, allocatable :: q_h2o_PEM_phys(:,:) ! Physics x Times: h2o mass mixing ratio computed in the PEM, first value comes from GCM [kg/kg] -integer :: timelen ! # time samples -real :: ave ! intermediate varibale to compute average -real, allocatable :: p(:,:) ! Physics x Atmosphere: pressure to recompute and write in restart (ngrid,llmp1) -real :: extra_mass ! Intermediate variables Extra mass of a tracer if it is greater than 1 +real :: ini_surf_h2o ! Initial surface of sublimating h2o ice +real :: ini_surf_co2 ! Initial surface of sublimating co2 ice +real :: global_ave_press_GCM ! constant: global average pressure retrieved in the GCM [Pa] +real :: global_ave_press_old ! constant: Global average pressure of initial/previous time step +real :: global_ave_press_new ! constant: Global average pressure of current time step +real, dimension(:,:), allocatable :: zplev_new ! Physical x Atmospheric field : mass of the atmospheric layers in the pem at current time step [kg/m^2] +real, dimension(:,:), allocatable :: zplev_gcm ! same but retrieved from the gcm [kg/m^2] +real, dimension(:,:,:), allocatable :: zplev_new_timeseries ! Physical x Atmospheric x Time: same as zplev_new, but in times series [kg/m ^2] +real, dimension(:,:,:), allocatable :: zplev_old_timeseries ! same but with the time series, for oldest time step +logical :: STOPPING_water ! Logical: is the criterion (% of change in the surface of sublimating water ice) reached? +logical :: STOPPING_1_water ! Logical: is there still water ice to sublimate? +logical :: STOPPING_co2 ! Logical: is the criterion (% of change in the surface of sublimating water ice) reached? +logical :: STOPPING_pressure ! Logical: is the criterion (% of change in the surface pressure) reached? +integer :: criterion_stop ! which criterion is reached ? 1= h2o ice surf, 2 = co2 ice surf, 3 = ps, 4 = orb param +real, save :: A , B, mmean ! Molar mass: intermediate A, B for computations of the mean molar mass of the layer [mol/kg] +real, dimension(:,:), allocatable :: vmr_co2_gcm ! Physics x Times co2 volume mixing ratio retrieve from the gcm [m^3/m^3] +real, dimension(:,:), allocatable :: vmr_co2_pem_phys ! Physics x Times co2 volume mixing ratio used in the PEM +real, dimension(:,:), allocatable :: q_co2_PEM_phys ! Physics x Times co2 mass mixing ratio in the first layer computed in the PEM, first value comes from GCM [kg/kg] +real, dimension(:,:), allocatable :: q_h2o_PEM_phys ! Physics x Times: h2o mass mixing ratio computed in the PEM, first value comes from GCM [kg/kg] +integer :: timelen ! # time samples +real :: ave ! intermediate varibale to compute average +real, dimension(:,:), allocatable :: p ! Physics x Atmosphere: pressure to recompute and write in restart (ngrid,llmp1) +real :: extra_mass ! Intermediate variables Extra mass of a tracer if it is greater than 1 ! Variables for slopes -real, dimension(:,:), allocatable :: min_co2_ice_1 ! ngrid field : minimum of co2 ice at each point for the first year [kg/m^2] -real, dimension(:,:), allocatable :: min_co2_ice_2 ! ngrid field : minimum of co2 ice at each point for the second year [kg/m^2] -real, dimension(:,:), allocatable :: min_h2o_ice_1 ! ngrid field : minimum of water ice at each point for the first year [kg/m^2] -real, dimension(:,:), allocatable :: min_h2o_ice_2 ! ngrid field : minimum of water ice at each point for the second year [kg/m^2] -real, dimension(:,:,:), allocatable :: co2_ice_GCM ! Physics x NSLOPE x Times field : co2 ice given by the GCM [kg/m^2] -real, dimension(:,:), allocatable :: initial_co2_ice_sublim ! physical point field : Logical array indicating sublimating point of co2 ice -real, dimension(:,:), allocatable :: initial_h2o_ice ! physical point field : Logical array indicating if there is water ice at initial state -real, dimension(:,:), allocatable :: initial_co2_ice ! physical point field : Logical array indicating if there is co2 ice at initial state -real, dimension(:,:), allocatable :: tendencies_co2_ice ! physical point xslope field : Tendency of evolution of perenial co2 ice over a year +real, dimension(:,:), allocatable :: min_co2_ice_1 ! ngrid field: minimum of co2 ice at each point for the first year [kg/m^2] +real, dimension(:,:), allocatable :: min_co2_ice_2 ! ngrid field: minimum of co2 ice at each point for the second year [kg/m^2] +real, dimension(:,:), allocatable :: min_h2o_ice_1 ! ngrid field: minimum of water ice at each point for the first year [kg/m^2] +real, dimension(:,:), allocatable :: min_h2o_ice_2 ! ngrid field: minimum of water ice at each point for the second year [kg/m^2] +real, dimension(:,:,:), allocatable :: co2_ice_GCM ! Physics x NSLOPE x Times field: co2 ice given by the GCM [kg/m^2] +real, dimension(:,:), allocatable :: initial_co2_ice_sublim ! physical point field: Logical array indicating sublimating point of co2 ice +real, dimension(:,:), allocatable :: initial_h2o_ice ! physical point field: Logical array indicating if there is water ice at initial state +real, dimension(:,:), allocatable :: initial_co2_ice ! physical point field: Logical array indicating if there is co2 ice at initial state +real, dimension(:,:), allocatable :: tendencies_co2_ice ! physical point x slope field: Tendency of evolution of perenial co2 ice over a year real, dimension(:,:), allocatable :: tendencies_co2_ice_ini ! physical point x slope field x nslope: Tendency of evolution of perenial co2 ice over a year in the GCM -real, dimension(:,:), allocatable :: tendencies_h2o_ice ! physical pointx slope field : Tendency of evolution of perenial h2o ice -real, dimension(:,:), allocatable :: flag_co2flow(:,:) ! (ngrid,nslope): Flag where there is a CO2 glacier flow -real, dimension(:), allocatable :: flag_co2flow_mesh(:) ! (ngrid) : Flag where there is a CO2 glacier flow -real, dimension(:,:), allocatable :: flag_h2oflow(:,:) ! (ngrid,nslope): Flag where there is a H2O glacier flow -real, dimension(:), allocatable :: flag_h2oflow_mesh(:) ! (ngrid) : Flag where there is a H2O glacier flow +real, dimension(:,:), allocatable :: tendencies_h2o_ice ! physical point x slope field: Tendency of evolution of perenial h2o ice +real, dimension(:,:), allocatable :: flag_co2flow ! (ngrid,nslope): Flag where there is a CO2 glacier flow +real, dimension(:), allocatable :: flag_co2flow_mesh ! (ngrid) : Flag where there is a CO2 glacier flow +real, dimension(:,:), allocatable :: flag_h2oflow ! (ngrid,nslope): Flag where there is a H2O glacier flow +real, dimension(:), allocatable :: flag_h2oflow_mesh ! (ngrid) : Flag where there is a H2O glacier flow ! Variables for surface and soil -real, allocatable :: tsurf_ave(:,:) ! Physic x SLOPE field : Averaged Surface Temperature [K] -real, allocatable :: tsoil_ave(:,:,:) ! Physic x SOIL x SLOPE field : Averaged Soil Temperature [K] -real, allocatable :: tsurf_GCM_timeseries(:,:,:) ! ngrid x SLOPE XTULES field : Surface Temperature in timeseries [K] -real, allocatable :: tsoil_phys_PEM_timeseries(:,:,:,:) ! IG x SLOPE XTULES field : NOn averaged Soil Temperature [K] -real, allocatable :: tsoil_GCM_timeseries(:,:,:,:) ! IG x SLOPE XTULES field : NOn averaged Soil Temperature [K] -real, allocatable :: tsurf_ave_yr1(:,:) ! Physic x SLOPE field : Averaged Surface Temperature of first call of the gcm [K] -real, allocatable :: TI_locslope(:,:) ! Physic x Soil: Intermediate thermal inertia to compute Tsoil [SI] -real, allocatable :: Tsoil_locslope(:,:) ! Physic x Soil: intermediate when computing Tsoil [K] -real, allocatable :: Tsurf_locslope(:) ! Physic x Soil: Intermediate surface temperature to compute Tsoil [K] -real, allocatable :: watersoil_density_timeseries(:,:,:,:) ! Physic x Soil x Slope x Times water soil density, time series [kg /m^3] -real, allocatable :: watersurf_density_ave(:,:) ! Physic x Slope, water surface density, yearly averaged [kg/m^3] -real, allocatable :: watersoil_density_PEM_timeseries(:,:,:,:) ! Physic x Soil x Slope x Times, water soil density, time series [kg/m^3] -real, allocatable :: watersoil_density_PEM_ave(:,:,:) ! Physic x Soil x SLopes, water soil density, yearly averaged [kg/m^3] -real, allocatable :: Tsurfave_before_saved(:,:) ! Surface temperature saved from previous time step [K] -real, allocatable :: delta_co2_adsorbded(:) ! Physics: quantity of CO2 that is exchanged because of adsorption / desorption [kg/m^2] -real, allocatable :: delta_h2o_adsorbded(:) ! Physics: quantity of H2O that is exchanged because of adsorption / desorption [kg/m^2] -real :: totmassco2_adsorbded ! Total mass of CO2 that is exchanged because of adsorption / desoprtion over the planets [kg] -real :: totmassh2o_adsorbded ! Total mass of H2O that is exchanged because of adsorption / desoprtion over the planets [kg] -logical :: bool_sublim ! logical to check if there is sublimation or not -real,allocatable :: porefillingice_thickness_prev_iter(:,:) ! ngrid x nslope: Thickness of the ice table at the previous iteration [m] -real,allocatable :: delta_h2o_icetablesublim(:) ! ngrid x Total mass of the H2O that has sublimated / condenses from the ice table [kg] +real, dimension(:,:), allocatable :: tsurf_ave ! Physic x SLOPE field : Averaged Surface Temperature [K] +real, dimension(:,:,:), allocatable :: tsoil_ave ! Physic x SOIL x SLOPE field : Averaged Soil Temperature [K] +real, dimension(:,:,:), allocatable :: tsurf_GCM_timeseries ! ngrid x SLOPE XTULES field : Surface Temperature in timeseries [K] +real, dimension(:,:,:,:), allocatable :: tsoil_phys_PEM_timeseries ! IG x SLOPE XTULES field : NOn averaged Soil Temperature [K] +real, dimension(:,:,:,:), allocatable :: tsoil_GCM_timeseries ! IG x SLOPE XTULES field : NOn averaged Soil Temperature [K] +real, dimension(:,:), allocatable :: tsurf_ave_yr1 ! Physic x SLOPE field : Averaged Surface Temperature of first call of the gcm [K] +real, dimension(:,:), allocatable :: TI_locslope ! Physic x Soil: Intermediate thermal inertia to compute Tsoil [SI] +real, dimension(:,:), allocatable :: Tsoil_locslope ! Physic x Soil: intermediate when computing Tsoil [K] +real, dimension(:), allocatable :: Tsurf_locslope ! Physic x Soil: Intermediate surface temperature to compute Tsoil [K] +real, dimension(:,:,:,:), allocatable :: watersoil_density_timeseries ! Physic x Soil x Slope x Times water soil density, time series [kg /m^3] +real, dimension(:,:), allocatable :: watersurf_density_ave ! Physic x Slope, water surface density, yearly averaged [kg/m^3] +real, dimension(:,:,:,:), allocatable :: watersoil_density_PEM_timeseries ! Physic x Soil x Slope x Times, water soil density, time series [kg/m^3] +real, dimension(:,:,:), allocatable :: watersoil_density_PEM_ave ! Physic x Soil x SLopes, water soil density, yearly averaged [kg/m^3] +real, dimension(:,:), allocatable :: Tsurfave_before_saved ! Surface temperature saved from previous time step [K] +real, dimension(:), allocatable :: delta_co2_adsorbded ! Physics: quantity of CO2 that is exchanged because of adsorption / desorption [kg/m^2] +real, dimension(:), allocatable :: delta_h2o_adsorbded ! Physics: quantity of H2O that is exchanged because of adsorption / desorption [kg/m^2] +real :: totmassco2_adsorbded ! Total mass of CO2 that is exchanged because of adsorption / desoprtion over the planets [kg] +real :: totmassh2o_adsorbded ! Total mass of H2O that is exchanged because of adsorption / desoprtion over the planets [kg] +logical :: bool_sublim ! logical to check if there is sublimation or not +real, dimension(:,:), allocatable :: porefillingice_thickness_prev_iter ! ngrid x nslope: Thickness of the ice table at the previous iteration [m] +real, dimension(:), allocatable :: delta_h2o_icetablesublim(:) ! ngrid x Total mass of the H2O that has sublimated / condenses from the ice table [kg] + ! Some variables for the PEM run real, parameter :: year_step = 1 ! timestep for the pem @@ -221,31 +222,46 @@ integer :: n_myear ! Maximum number of Martian years of the chained simulations real :: timestep ! timestep [s] -real :: watercap_sum ! total mass of water cap [kg/m^2] -real :: water_sum ! total mass of water in the mesh [kg/m^2] +real :: watercap_sum ! total mass of water cap [kg/m^2] +real :: water_sum ! total mass of water in the mesh [kg/m^2] #ifdef CPP_STD - real :: frost_albedo_threshold = 0.05 ! frost albedo threeshold to convert fresh frost to old ice - real :: albedo_h2o_frost ! albedo of h2o frost - real, allocatable :: tsurf_read_generic(:) ! Temporary variable to do the subslope transfert dimensiion when reading form generic - real, allocatable :: qsurf_read_generic(:,:) ! Temporary variable to do the subslope transfert dimensiion when reading form generic - real, allocatable :: tsoil_read_generic(:,:) ! Temporary variable to do the subslope transfert dimensiion when reading form generic - real, allocatable :: emis_read_generic(:) ! Temporary variable to do the subslope transfert dimensiion when reading form generic - real, allocatable :: albedo_read_generic(:,:) ! Temporary variable to do the subslope transfert dimensiion when reading form generic - real, allocatable :: tsurf(:,:) ! Subslope variable, only needed in the GENERIC case - real, allocatable :: qsurf(:,:,:) ! Subslope variable, only needed in the GENERIC case - real, allocatable :: tsoil(:,:,:) ! Subslope variable, only needed in the GENERIC case - real, allocatable :: emis(:,:) ! Subslope variable, only needed in the GENERIC case - real, allocatable :: watercap(:,:) ! Subslope variable, only needed in the GENERIC case =0 no watercap in generic model - logical, allocatable :: WATERCAPTAG(:) ! Subslope variable, only needed in the GENERIC case =false no watercaptag in generic model - real, allocatable :: albedo(:,:,:) ! Subslope variable, only needed in the GENERIC case - real, allocatable :: inertiesoil(:,:,:) ! Subslope variable, only needed in the GENERIC case + real :: frost_albedo_threshold = 0.05 ! frost albedo threeshold to convert fresh frost to old ice + real :: albedo_h2o_frost ! albedo of h2o frost + real, dimension(:), allocatable :: tsurf_read_generic ! Temporary variable to do the subslope transfert dimensiion when reading form generic + real, dimension(:,:), allocatable :: qsurf_read_generic ! Temporary variable to do the subslope transfert dimensiion when reading form generic + real, dimension(:,:), allocatable :: tsoil_read_generic ! Temporary variable to do the subslope transfert dimensiion when reading form generic + real, dimension(:), allocatable :: emis_read_generic ! Temporary variable to do the subslope transfert dimensiion when reading form generic + real, dimension(:,:), allocatable :: albedo_read_generic ! Temporary variable to do the subslope transfert dimensiion when reading form generic + real, dimension(:,:), allocatable :: tsurf ! Subslope variable, only needed in the GENERIC case + real, dimension(:,:,:), allocatable :: qsurf ! Subslope variable, only needed in the GENERIC case + real, dimension(:,:,:), allocatable :: tsoil ! Subslope variable, only needed in the GENERIC case + real, dimension(:,:), allocatable :: emis ! Subslope variable, only needed in the GENERIC case + real, dimension(:,:), allocatable :: watercap ! Subslope variable, only needed in the GENERIC case =0 no watercap in generic model + logical, dimension(:), allocatable :: WATERCAPTAG ! Subslope variable, only needed in the GENERIC case =false no watercaptag in generic model + real, dimension(:,:,:), allocatable :: albedo ! Subslope variable, only needed in the GENERIC case + real, dimension(:,:,:), allocatable :: inertiesoil ! Subslope variable, only needed in the GENERIC case #endif #ifdef CPP_1D - integer :: nsplit_phys - integer, parameter :: jjm_value = jjm - 1 - integer :: ierr + integer :: nsplit_phys + integer, parameter :: jjm_value = jjm - 1 + integer :: ierr + + ! Dummy variables to use the subroutine 'init_testphys1d' + real, dimension(:,:), allocatable :: q_tmp ! Temporary tracer mixing ratio (e.g. kg/kg) + logical :: startfiles_1D, therestart1D, therestartfi + integer :: ndt, day0 + real :: ptif, pks, day, gru, grv, atm_wat_profile, atm_wat_tau + real, dimension(:), allocatable :: zqsat, qsurf_tmp + real, dimension(:,:), allocatable :: dq, dqdyn + real, dimension(nlayer) :: play, w + real, dimension(nlayer + 1) :: plev, q2_tmp + real, dimension(nsoilmx) :: tsoil_tmp + real, dimension(1) :: tsurf_tmp, emis_tmp + real, dimension(1,1) :: albedo_tmp #else - integer, parameter :: jjm_value = jjm + integer, parameter :: jjm_value = jjm + real, dimension(:), allocatable :: ap ! Coefficient ap read in FILE_NAME_start and written in restart + real, dimension(:), allocatable :: bp ! Coefficient bp read in FILE_NAME_start and written in restart #endif @@ -262,14 +278,10 @@ #endif +! Some constants day_ini = 0 ! test time_phys = 0. ! test - -! Some constants ngrid = ngridmx -nlayer = llm - A = (1/m_co2 - 1/m_noco2) B = 1/m_noco2 - year_day = 669 daysec = 88775. @@ -287,36 +299,14 @@ nq = nqtot allocate(q(ip1jmp1,llm,nqtot)) + allocate(longitude(ngrid),latitude(ngrid),cell_area(ngrid)) #else - ! load tracer names from file 'traceur.def' - open(90,file = 'traceur.def',status = 'old',form = 'formatted',iostat = ierr) - if (ierr /= 0) then - write(*,*) 'Cannot find required file "traceur.def"' - write(*,*) ' If you want to run with tracers, I need it' - write(*,*) ' ... might as well stop here ...' - stop - else - write(*,*) "pem1d: Reading file traceur.def" - ! read number of tracers: - read(90,*,iostat = ierr) nq - write(*,*) "nq",nq - nqtot = nq ! set value of nqtot (in infotrac module) as nq - if (ierr /= 0) then - write(*,*) "pem1d: error reading number of tracers" - write(*,*) " (first line of traceur.def) " - stop - endif - if (nq < 1) then - write(*,*) "pem1d: error number of tracers" - write(*,*) "is nq=",nq," but must be >=1!" - stop - endif - endif - nq = nqtot + allocate(longitude(1),latitude(1),cell_area(1)) + call init_testphys1d(.true.,ngrid,nlayer,610.,nq,q_tmp,time_0,ps(1),ucov,vcov,teta,startfiles_1D,therestart1D,therestartfi, & + ndt,ptif,pks,dtphys,zqsat,qsurf_tmp,dq,dqdyn,day0,day,tsurf_tmp,gru,grv,w,q2_tmp,play,plev,tsoil_tmp, & + albedo_tmp,emis_tmp,latitude,longitude,cell_area,atm_wat_profile,atm_wat_tau) + ps(2) = ps(1) allocate(q(ip1jmp1,llm,nqtot)) - allocate(ap(nlayer + 1)) - allocate(bp(nlayer + 1)) - call init_pem1d(llm,nqtot,ucov,vcov,teta,q,ps,time_0,ap,bp) - pi = 2.*asin(1.) - g = 3.72 + q(ip1jmp1,:,:) = q_tmp(:,:) + deallocate(q_tmp,qsurf_tmp) nsplit_phys = 1 #endif @@ -356,8 +346,4 @@ ! Here we simply read them in the startfi_evol.nc file status =nf90_open(FILE_NAME, NF90_NOWRITE, ncid) - -allocate(longitude(ngrid)) -allocate(latitude(ngrid)) -allocate(cell_area(ngrid)) status = nf90_inq_varid(ncid,"longitude",lonvarid) @@ -414,5 +400,5 @@ qsurf_read_generic,cloudfrac,totcloudfrac,hice,rnat,pctsrf_sic, & tslab, tsea_ice,sea_ice) - call surfini(ngrid,nq,qsurf_read_generic,albedo_read_generic,albedo_bareground,albedo_snow_SPECTV,albedo_co2_ice_SPECTV) + call surfini(ngrid,nq,qsurf_read_generic,albedo_read_generic,albedo_bareground,albedo_snow_SPECTV,albedo_co2_ice_SPECTV) nslope = 1 @@ -452,5 +438,5 @@ endif if (noms(nnq) == "co2") igcm_co2 = nnq -enddo +enddo r = 8.314511*1000./mugaz @@ -473,7 +459,7 @@ allocate(flag_h2oflow_mesh(ngrid)) -flag_co2flow(:,:) = 0 +flag_co2flow(:,:) = 0 flag_co2flow_mesh(:) = 0 -flag_h2oflow(:,:) = 0 +flag_h2oflow(:,:) = 0 flag_h2oflow_mesh(:) = 0 @@ -614,5 +600,5 @@ global_ave_press_old = 0. do i = 1,ngrid - global_ave_press_old = global_ave_press_old + cell_area(i)*ps_start_GCM(i)/Total_surface + global_ave_press_old = global_ave_press_old + cell_area(i)*ps_start_GCM(i)/Total_surface enddo @@ -656,5 +642,5 @@ write(*,*) "Tot mass of H2O in the regolith=", totmassh2o_adsorbded endif ! adsorption -deallocate(tsurf_ave_yr1) +deallocate(tsurf_ave_yr1) !------------------------ @@ -678,5 +664,5 @@ !------------------------ ! II Run -! II_a Update pressure, ice and tracers +! II_a Update pressure, ice and tracers !------------------------ year_iter = 0 @@ -688,9 +674,9 @@ do i = 1,ngrid do islope = 1,nslope - global_ave_press_new = global_ave_press_new-g*cell_area(i)*tendencies_co2_ice(i,islope)*subslope_dist(i,islope)/cos(pi*def_slope_mean(islope)/180.)/Total_surface - enddo - enddo - call WRITEDIAGFI(ngrid,'ps_ave','Global average pressure','Pa',0,global_ave_press_new) - + global_ave_press_new = global_ave_press_new - g*cell_area(i)*tendencies_co2_ice(i,islope)*subslope_dist(i,islope)/cos(pi*def_slope_mean(islope)/180.)/Total_surface + enddo + enddo + call writediagfi(ngrid,'ps_ave','Global average pressure','Pa',0,global_ave_press_new) + if (adsorption_pem) then do i = 1,ngrid @@ -739,5 +725,5 @@ do ig = 1,ngrid - do t = 1, timelen + do t = 1,timelen q_co2_PEM_phys(ig,t) = q_co2_PEM_phys(ig,t)*(zplev_old_timeseries(ig,l,t) - zplev_old_timeseries(ig,l + 1,t))/ & (zplev_new_timeseries(ig,l,t) - zplev_new_timeseries(ig,l + 1,t)) & @@ -767,10 +753,10 @@ call evol_co2_ice_s(qsurf(:,igcm_co2,:),tendencies_co2_ice,iim,jjm_value,ngrid,cell_area,nslope) - do islope=1, nslope + do islope = 1,nslope write(str2(1:2),'(i2.2)') islope - call WRITEDIAGFI(ngrid,'h2o_ice_s_slope'//str2,'H2O ice','kg.m-2',2,qsurf(:,igcm_h2o_ice,islope)) - call WRITEDIAGFI(ngrid,'tendencies_h2o_ice_slope'//str2,'H2O ice tend','kg.m-2.year-1',2,tendencies_h2o_ice(:,islope)) - call WRITEDIAGFI(ngrid,'c2ice_slope'//str2,'CO2 ice','kg.m-2',2,qsurf(:,igcm_co2,islope)) - call WRITEDIAGFI(ngrid,'tendencies_co2_ice_slope'//str2,'CO2 ice tend','kg.m-2.year-1',2,tendencies_co2_ice(:,islope)) + call writediagfi(ngrid,'h2o_ice_s_slope'//str2,'H2O ice','kg.m-2',2,qsurf(:,igcm_h2o_ice,islope)) + call writediagfi(ngrid,'tendencies_h2o_ice_slope'//str2,'H2O ice tend','kg.m-2.year-1',2,tendencies_h2o_ice(:,islope)) + call writediagfi(ngrid,'c2ice_slope'//str2,'CO2 ice','kg.m-2',2,qsurf(:,igcm_co2,islope)) + call writediagfi(ngrid,'tendencies_co2_ice_slope'//str2,'CO2 ice tend','kg.m-2.year-1',2,tendencies_co2_ice(:,islope)) enddo @@ -783,5 +769,5 @@ if (co2glaciersflow) call co2glaciers_evol(timelen,ngrid,nslope,iflat,subslope_dist,def_slope_mean,vmr_co2_pem_phys,ps_timeseries, & global_ave_press_GCM,global_ave_press_new,qsurf(:,igcm_co2,:),flag_co2flow,flag_co2flow_mesh) - + write(*,*) "H2O glacier flows" @@ -790,7 +776,7 @@ do islope = 1,nslope write(str2(1:2),'(i2.2)') islope - call WRITEDIAGFI(ngrid,'co2ice_slope'//str2,'CO2 ice','kg.m-2',2,qsurf(:,igcm_co2,islope)) - call WRITEDIAGFI(ngrid,'tendencies_co2_ice_slope'//str2,'CO2 ice tend','kg.m-2.year-1',2,tendencies_co2_ice(:,islope)) - call WRITEDIAGFI(ngrid,'Flow_co2ice_slope'//str2,'CO2 ice flow','Boolean',2,flag_co2flow(:,islope)) + call writediagfi(ngrid,'co2ice_slope'//str2,'CO2 ice','kg.m-2',2,qsurf(:,igcm_co2,islope)) + call writediagfi(ngrid,'tendencies_co2_ice_slope'//str2,'CO2 ice tend','kg.m-2.year-1',2,tendencies_co2_ice(:,islope)) + call writediagfi(ngrid,'Flow_co2ice_slope'//str2,'CO2 ice flow','Boolean',2,flag_co2flow(:,islope)) enddo @@ -835,5 +821,5 @@ else albedo(ig,1,islope) = albedodat(ig) - albedo(ig,2,islope) = albedodat(ig) + albedo(ig,2,islope) = albedodat(ig) emis(ig,islope) = emissiv endif @@ -862,5 +848,5 @@ do t = 1,timelen - tsurf_GCM_timeseries(:,:,t) = tsurf_GCM_timeseries(:,:,t) + (tsurf_ave(:,:) - Tsurfave_before_saved(:,:)) + tsurf_GCM_timeseries(:,:,t) = tsurf_GCM_timeseries(:,:,t) + (tsurf_ave(:,:) - Tsurfave_before_saved(:,:)) enddo ! for the start @@ -873,5 +859,5 @@ do islope = 1,nslope write(str2(1:2),'(i2.2)') islope - call WRITEDIAGFI(ngrid,'tsurf_slope'//str2,'tsurf','K',2,tsurf(:,islope)) + call writediagfi(ngrid,'tsurf_slope'//str2,'tsurf','K',2,tsurf(:,islope)) enddo @@ -896,5 +882,5 @@ do isoil = 1,nsoilmx_PEM watersoil_density_PEM_timeseries(ig,isoil,islope,t) = exp(beta_clap_h2o/Tsoil_locslope(ig,isoil) + alpha_clap_h2o)/Tsoil_locslope(ig,isoil)*mmol(igcm_h2o_vap)/(mugaz*r) - if (isnan(Tsoil_locslope(ig,isoil))) call abort_pem("PEM - Update Tsoil","NaN detected in Tsoil ",1) + if (isnan(Tsoil_locslope(ig,isoil))) call abort_pem("PEM - Update Tsoil","NaN detected in Tsoil ",1) enddo enddo @@ -1001,5 +987,5 @@ endif - global_ave_press_old=global_ave_press_new + global_ave_press_old = global_ave_press_new enddo ! big time iteration loop of the pem @@ -1037,5 +1023,5 @@ water_sum = water_sum + qsurf(ig,igcm_h2o_ice,islope)*subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.) enddo - if (.not. watercaptag(ig)) then ! let's check if we have an 'infinite' source of water that has been forming. + if (.not. watercaptag(ig)) then ! let's check if we have an 'infinite' source of water that has been forming. if (water_sum > threshold_water_frost2perenial) then ! the overall mesh can be considered as an infite source of water. No need to change the albedo: done in II.d.1 watercaptag(ig) = .true. @@ -1071,15 +1057,10 @@ if (soil_pem) then call interpolate_TIPEM_TIGCM(ngrid,nslope,nsoilmx_PEM,nsoilmx,TI_PEM,inertiesoil) - tsoil(:,:,:) = tsoil_phys_PEM_timeseries(:,1:nsoilmx,:,timelen) + tsoil(:,:,:) = tsoil_phys_PEM_timeseries(:,1:nsoilmx,:,timelen) endif ! III_a.4 Pressure (for start) -do i = 1,ip1jmp1 - ps(i) = ps(i)*global_ave_press_new/global_ave_press_GCM -enddo - -do i = 1,ngrid - ps_start_GCM(i) = ps_start_GCM(i)*global_ave_press_new/global_ave_press_GCM -enddo +ps(:) = ps(:)*global_ave_press_new/global_ave_press_GCM +ps_start_GCM(:) = ps_start_GCM(:)*global_ave_press_new/global_ave_press_GCM ! III_a.5 Tracer (for start) @@ -1087,7 +1068,5 @@ do l = 1,nlayer + 1 - do ig = 1,ngrid - zplev_new(ig,l) = ap(l) + bp(l)*ps_start_GCM(ig) - enddo + zplev_new(:,l) = ap(l) + bp(l)*ps_start_GCM(:) enddo @@ -1104,5 +1083,5 @@ do ig = 1,ngrid q(ig,l,nnq) = q(ig,l,nnq)*(zplev_gcm(ig,l) - zplev_gcm(ig,l + 1))/(zplev_new(ig,l) - zplev_new(ig,l + 1)) & - + ((zplev_new(ig,l) - zplev_new(ig,l + 1)) - (zplev_gcm(ig,l) - zplev_gcm(ig,l + 1)))/ & + + ((zplev_new(ig,l) - zplev_new(ig,l + 1)) - (zplev_gcm(ig,l) - zplev_gcm(ig,l + 1)))/ & (zplev_new(ig,l) - zplev_new(ig,l + 1)) enddo @@ -1117,7 +1096,7 @@ do l = 1,llm - 1 if (q(ig,l,nnq) > 1 .and. (noms(nnq) /= "dust_number") .and. (noms(nnq) /= "ccn_number") .and. (noms(nnq) /= "stormdust_number") .and. (noms(nnq) /= "topdust_number")) then - extra_mass=(q(ig,l,nnq)-1)*(zplev_new(ig,l)-zplev_new(ig,l+1)) - q(ig,l,nnq)=1. - q(ig,l+1,nnq)=q(ig,l+1,nnq)+extra_mass*(zplev_new(ig,l+1)-zplev_new(ig,l+2)) + extra_mass = (q(ig,l,nnq) - 1)*(zplev_new(ig,l) - zplev_new(ig,l + 1)) + q(ig,l,nnq) = 1. + q(ig,l + 1,nnq) = q(ig,l + 1,nnq) + extra_mass*(zplev_new(ig,l + 1) - zplev_new(ig,l + 2)) write(*,*) 'extra ',noms(nnq),extra_mass, noms(nnq) /= "dust_number",noms(nnq) /= "ccn_number" endif @@ -1146,7 +1125,6 @@ write(*,*) "restart_evol.nc has been written" #else - do nnq = 1, nqtot - call writeprofile(nlayer,q(1,:,nnq),noms(nnq),nnq,qsurf) - enddo + call writerestart1D('restart1D_evol.txt',ps(1),tsurf,nlayer,teta,ucov,vcov,nq,noms,qsurf,q) + write(*,*) "restart1D_evol.txt has been written" #endif @@ -1198,3 +1176,2 @@ END PROGRAM pem -