Index: /trunk/LMDZ.COMMON/libf/evolution/criterion_pem_stop_mod.F90 =================================================================== --- /trunk/LMDZ.COMMON/libf/evolution/criterion_pem_stop_mod.F90 (revision 3049) +++ /trunk/LMDZ.COMMON/libf/evolution/criterion_pem_stop_mod.F90 (revision 3050) @@ -63,9 +63,9 @@ present_surf.GT.ini_surf*(1+water_ice_criterion)) then STOPPING=.TRUE. - print *, "Reason of stopping : The surface of water ice sublimating reach the threshold:" - print *, "Current surface of water ice sublimating=", present_surf - print *, "Initial surface of water ice sublimating=", ini_surf - print *, "Percentage of change accepted=", water_ice_criterion*100 - print *, "present_surf=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: unk/LMDZ.COMMON/libf/evolution/init_phys_1d_mod.F90 =================================================================== --- /trunk/LMDZ.COMMON/libf/evolution/init_phys_1d_mod.F90 (revision 3049) +++ (revision ) @@ -1,587 +1,0 @@ -module init_phys_1d_mod - -implicit none - -contains - -#ifdef CPP_1D - -subroutine init_phys_1d(llm_in,nqtot_in,v,u,temp,q,psurf,time,ap_out,bp_out) - -!-------------- init_phys_1d ------------- -! -! 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 not found' + write(*,*) 'get_field_rgen: Field <'//field_name//'> not found' CALL abort_physic("get_field_rgen","Failed to get field",1) ENDIF @@ -329,5 +329,5 @@ IF (ierr/=NF90_NOERR) THEN ! La variable exist dans le fichier mais la lecture a echouee. - PRINT*, 'get_field_rgen: Failed reading <'//field_name//'>' + write(*,*) 'get_field_rgen: Failed reading <'//field_name//'>' ! IF (field_name=='CLWCON' .OR. field_name=='RNEBCON' .OR. field_name=='RATQS') THEN @@ -336,8 +336,8 @@ ! ierr=NF90_GET_VAR(nid_start,varid,field_glo(1:klon_glo)) ! IF (ierr/=NF90_NOERR) THEN -! PRINT*, 'phyetat0: Lecture echouee aussi en 2D pour <'//field_name//'>' +! write(*,*) 'phyetat0: Lecture echouee aussi en 2D pour <'//field_name//'>' ! CALL abort ! ELSE -! PRINT*, 'phyetat0: La variable <'//field_name//'> lu sur surface seulement'!, selon ancien format, le reste mis a zero' +! write(*,*) 'phyetat0: La variable <'//field_name//'> lu sur surface seulement'!, selon ancien format, le reste mis a zero' ! END IF ! ELSE @@ -436,5 +436,5 @@ ierr=NF90_GET_VAR(nid_start,varid,var) IF (ierr/=NF90_NOERR) THEN - PRINT*, 'phyetat0: Failed loading <'//trim(var_name)//'>' + write(*,*) 'phyetat0: Failed loading <'//trim(var_name)//'>' CALL abort_physic("get_var_rgen","Failed to read variable",1) ENDIF @@ -456,5 +456,5 @@ ELSE IF (.NOT. tmp_found) THEN - PRINT*, 'phyetat0: Variable <'//trim(var_name)//'> not found' + write(*,*) 'phyetat0: Variable <'//trim(var_name)//'> not found' CALL abort_physic("get_var_rgen","Failed to read variable",1) ENDIF @@ -900,5 +900,5 @@ ELSE - PRINT *, "Error phyredem(put_field_rgen) : wrong dimension for ",trim(field_name) + write(*,*) "Error phyredem(put_field_rgen) : wrong dimension for ",trim(field_name) write(*,*) " field_size =",field_size CALL abort_physic("put_field_rgen","wrong field dimension",1) @@ -1026,5 +1026,5 @@ idim1d=idim9 ELSE - PRINT *, "put_var_rgen error : wrong dimension" + write(*,*) "put_var_rgen error : wrong dimension" write(*,*) " var_size =",var_size CALL abort_physic("get_var_rgen","Wrong variable dimension",1) Index: /trunk/LMDZ.COMMON/libf/evolution/nb_time_step_GCM.F90 =================================================================== --- /trunk/LMDZ.COMMON/libf/evolution/nb_time_step_GCM.F90 (revision 3049) +++ /trunk/LMDZ.COMMON/libf/evolution/nb_time_step_GCM.F90 (revision 3050) @@ -66,5 +66,5 @@ CALL err(NF90_CLOSE(fID),"close",fichnom) - print *, "The number of timestep of the PCM run data=", timelen + write(*,*) "The number of timestep of the PCM run data=", timelen CONTAINS Index: /trunk/LMDZ.COMMON/libf/evolution/orbit_param_criterion_mod.F90 =================================================================== --- /trunk/LMDZ.COMMON/libf/evolution/orbit_param_criterion_mod.F90 (revision 3049) +++ /trunk/LMDZ.COMMON/libf/evolution/orbit_param_criterion_mod.F90 (revision 3050) @@ -79,5 +79,5 @@ close(73) if (last_ilask == 0 .or. last_ilask == size(yearlask,1) + 1) then - write(*,*) 'The current year could not be find in the file "obl_ecc_lsp.asc".' + write(*,*) 'The current year could not be found in the file "obl_ecc_lsp.asc".' stop else @@ -196,7 +196,7 @@ enddo -if (.not. found_obl) write(*,*) 'The maximum reachable year for obl. could not be find in the file "obl_ecc_lsp.asc".' -if (.not. found_ecc) write(*,*) 'The maximum reachable year for ecc. could not be find in the file "obl_ecc_lsp.asc".' -if (.not. found_lsp) write(*,*) 'The maximum reachable year for Lsp could not be find in the file "obl_ecc_lsp.asc".' +if (.not. found_obl) write(*,*) 'The maximum reachable year for obl. could not be found in the file "obl_ecc_lsp.asc".' +if (.not. found_ecc) write(*,*) 'The maximum reachable year for ecc. could not be found in the file "obl_ecc_lsp.asc".' +if (.not. found_lsp) write(*,*) 'The maximum reachable year for Lsp could not be found in the file "obl_ecc_lsp.asc".' write(*,*) 'Max. number of iterations for the obl. parameter =', max_obl_iter Index: /trunk/LMDZ.COMMON/libf/evolution/pem.F90 =================================================================== --- /trunk/LMDZ.COMMON/libf/evolution/pem.F90 (revision 3049) +++ /trunk/LMDZ.COMMON/libf/evolution/pem.F90 (revision 3050) @@ -96,5 +96,5 @@ use physics_distribution_mod, only: init_physics_distribution use mod_grid_phy_lmdz, only: regular_lonlat - use init_phys_1d_mod, only: init_phys_1d + use init_pem1d_mod, only: init_pem1d #endif @@ -316,7 +316,7 @@ allocate(ap(nlayer + 1)) allocate(bp(nlayer + 1)) - call init_phys_1d(llm,nqtot,vcov,ucov,teta,q,ps,time_0,ap,bp) - pi = 2.*asin(1.) - g = 3.72 + call init_pem1d(llm,nqtot,ucov,vcov,teta,q,ps,time_0,ap,bp) + pi = 2.*asin(1.) + g = 3.72 nsplit_phys = 1 #endif @@ -663,7 +663,7 @@ !------------------------ #ifndef CPP_STD - call iniorbit(aphelie,periheli,year_day,peri_day,obliquit) + call iniorbit(aphelie,periheli,year_day,peri_day,obliquit) #else - call iniorbit(apoastr,periastr,year_day,peri_day,obliquit) + call iniorbit(apoastr,periastr,year_day,peri_day,obliquit) #endif @@ -691,5 +691,4 @@ enddo enddo - write(*,*) 'Global average pressure old time step',global_ave_press_old call WRITEDIAGFI(ngrid,'ps_ave','Global average pressure','Pa',0,global_ave_press_new) @@ -697,8 +696,8 @@ do i = 1,ngrid global_ave_press_new = global_ave_press_new - g*cell_area(i)*delta_co2_adsorbded(i)/Total_surface - enddo - write(*,*) 'Global average pressure old time step',global_ave_press_old - write(*,*) 'Global average pressure new time step',global_ave_press_new - endif + enddo + endif + write(*,*) 'Global average pressure old time step',global_ave_press_old + write(*,*) 'Global average pressure new time step',global_ave_press_new ! II.a.2. Old pressure levels for the timeseries, this value is deleted when unused and recreated each time (big memory consuption) @@ -1180,6 +1179,4 @@ call pemdem0("restartfi_PEM.nc",longitude,latitude,cell_area,nsoilmx_PEM,ngrid, & float(day_ini),0.,nslope,def_slope,subslope_dist) - - call pemdem1("restartfi_PEM.nc",i_myear,nsoilmx_PEM,ngrid,nslope,tsoil_PEM, & TI_PEM, porefillingice_depth,porefillingice_thickness, & Index: /trunk/LMDZ.COMMON/libf/evolution/read_data_GCM.F90 =================================================================== --- /trunk/LMDZ.COMMON/libf/evolution/read_data_GCM.F90 (revision 3049) +++ /trunk/LMDZ.COMMON/libf/evolution/read_data_GCM.F90 (revision 3050) @@ -80,5 +80,5 @@ B=1/m_noco2 - print *, "Opening ", fichnom, "..." + write(*,*) "Opening ", fichnom, "..." ! Open initial state NetCDF file @@ -86,24 +86,24 @@ CALL err(NF90_OPEN(var,NF90_NOWRITE,fID),"open",var) - print *, "Downloading data for vmr co2..." + write(*,*) "Downloading data for vmr co2..." CALL get_var3("co2_layer1" ,q_co2_dyn) - print *, "Downloading data for vmr co2 done" - print *, "Downloading data for vmr h20..." + write(*,*) "Downloading data for vmr co2 done" + write(*,*) "Downloading data for vmr h20..." CALL get_var3("h2o_layer1" ,q_h2o_dyn) - print *, "Downloading data for vmr h2o done" - print *, "Downloading data for surface pressure ..." + write(*,*) "Downloading data for vmr h2o done" + write(*,*) "Downloading data for surface pressure ..." CALL get_var3("ps" ,ps_GCM) - print *, "Downloading data for surface pressure done" - print *, "nslope=", nslope + write(*,*) "Downloading data for surface pressure done" + write(*,*) "nslope=", nslope if(nslope.gt.1) then - print *, "Downloading data for co2ice_slope ..." + write(*,*) "Downloading data for co2ice_slope ..." DO islope=1,nslope @@ -112,6 +112,6 @@ ENDDO - print *, "Downloading data for co2ice_slope done" - print *, "Downloading data for h2o_ice_s_slope ..." + write(*,*) "Downloading data for co2ice_slope done" + write(*,*) "Downloading data for h2o_ice_s_slope ..." DO islope=1,nslope @@ -120,9 +120,9 @@ ENDDO - print *, "Downloading data for h2o_ice_s_slope done" + write(*,*) "Downloading data for h2o_ice_s_slope done" #ifndef CPP_STD - print *, "Downloading data for watercap_slope ..." + write(*,*) "Downloading data for watercap_slope ..." DO islope=1,nslope write(num,fmt='(i2.2)') islope @@ -130,9 +130,9 @@ ! watercap_slope(:,:,:,:)= 0. ENDDO - print *, "Downloading data for watercap_slope done" + write(*,*) "Downloading data for watercap_slope done" #endif - print *, "Downloading data for tsurf_slope ..." + write(*,*) "Downloading data for tsurf_slope ..." DO islope=1,nslope @@ -141,5 +141,5 @@ ENDDO - print *, "Downloading data for tsurf_slope done" + write(*,*) "Downloading data for tsurf_slope done" #ifndef CPP_STD @@ -147,5 +147,5 @@ if(soil_pem) then - print *, "Downloading data for soiltemp_slope ..." + write(*,*) "Downloading data for soiltemp_slope ..." DO islope=1,nslope @@ -154,7 +154,7 @@ ENDDO - print *, "Downloading data for soiltemp_slope done" - - print *, "Downloading data for watersoil_density ..." + write(*,*) "Downloading data for soiltemp_slope done" + + write(*,*) "Downloading data for watersoil_density ..." DO islope=1,nslope @@ -163,7 +163,7 @@ ENDDO - print *, "Downloading data for watersoil_density done" - - print *, "Downloading data for watersurf_density ..." + write(*,*) "Downloading data for watersoil_density done" + + write(*,*) "Downloading data for watersurf_density ..." DO islope=1,nslope @@ -172,5 +172,5 @@ ENDDO - print *, "Downloading data for watersurf_density done" + write(*,*) "Downloading data for watersurf_density done" endif !soil_pem @@ -196,12 +196,12 @@ ! Compute the minimum over the year for each point - print *, "Computing the min of h2o_ice_slope" + write(*,*) "Computing the min of h2o_ice_slope" min_h2o_ice_dyn(:,:,:)=minval(h2o_ice_s_dyn+watercap_slope,4) - print *, "Computing the min of co2_ice_slope" + write(*,*) "Computing the min of co2_ice_slope" min_co2_ice_dyn(:,:,:)=minval(co2_ice_slope_dyn,4) !Compute averages - print *, "Computing average of tsurf" + write(*,*) "Computing average of tsurf" tsurf_ave_dyn(:,:,:)=SUM(tsurf_gcm_dyn(:,:,:,:),4)/timelen @@ -215,7 +215,7 @@ if(soil_pem) then - print *, "Computing average of tsoil" + write(*,*) "Computing average of tsoil" tsoil_ave_dyn(:,:,:,:)=SUM(tsoil_gcm_dyn(:,:,:,:,:),5)/timelen - print *, "Computing average of watersurf_density" + write(*,*) "Computing average of watersurf_density" watersurf_density_ave(:,:) = SUM(watersurf_density(:,:,:),3)/timelen endif Index: /trunk/LMDZ.COMMON/libf/evolution/reshape_XIOS_output.F90 =================================================================== --- /trunk/LMDZ.COMMON/libf/evolution/reshape_XIOS_output.F90 (revision 3049) +++ /trunk/LMDZ.COMMON/libf/evolution/reshape_XIOS_output.F90 (revision 3050) @@ -101,5 +101,5 @@ do i=1,nvars status = nf90_inquire_variable(ncid1, varids(i), name=namevar, xtype=xtype_var, ndims = numdims,natts = numatts) - print *, "namevar00= ", namevar + write(*,*) "namevar00= ", namevar if(status /= nf90_noerr) call handle_err(status) allocate(dimid_var(numdims)) Index: /trunk/LMDZ.MARS/changelog.txt =================================================================== --- /trunk/LMDZ.MARS/changelog.txt (revision 3049) +++ /trunk/LMDZ.MARS/changelog.txt (revision 3050) @@ -4205,2 +4205,5 @@ Improvements of scripts to launch the chained simulations of GCM and PEM runs. Correction of a case where maximum admissible change of orbital parameters could not be found, in particular for Lsp because of modulo + some improvements. + +== 26/09/2023 == JBC +Minor changes concerning the form of the code in the PEM. Index: /trunk/LMDZ.MARS/deftank/pem/launch_pem.sh =================================================================== --- /trunk/LMDZ.MARS/deftank/pem/launch_pem.sh (revision 3049) +++ /trunk/LMDZ.MARS/deftank/pem/launch_pem.sh (revision 3050) @@ -51,12 +51,4 @@ if [ ! -f "run_GCM.def" ]; then echo "Error: file \"run_GCM.def\" does not exist in $dir!" - exit 1 -fi -if [ ! -f "diagfi_PEM.def" ]; then - echo "Error: file \"diagfi_PEM.def\" does not exist in $dir!" - exit 1 -fi -if [ ! -f "diagfi_GCM.def" ]; then - echo "Error: file \"diagfi_GCM.def\" does not exist in $dir!" exit 1 fi @@ -118,5 +110,5 @@ cp run_GCM.def run.def rm diagfi.def - if [ ! -f "diagfi_GCM.def" ]; then + if [ -f "diagfi_GCM.def" ]; then cp diagfi_GCM.def diagfi.def fi @@ -146,4 +138,8 @@ cp restart.nc starts/restart${iGCM}.nc mv restart.nc start.nc + fi + if [ -f "restart1D.txt" ]; then + cp restart1D.txt starts/restart1D${iGCM}.txt + mv restart1D.txt start1D.txt fi if ls profile_out_* 1> /dev/null 2>&1; then @@ -165,12 +161,18 @@ done echo "Done!" - #--- Running PEM + #--- Running PEM echo "Run PEM $iPEM..." cp run_PEM.def run.def rm diagfi.def - if [ ! -f "diagfi_PEM.def" ]; then + if [ -f "diagfi_PEM.def" ]; then cp diagfi_PEM.def diagfi.def fi mv startfi.nc startfi_evol.nc + if [ -f "start.nc" ]; then + cp start.nc start_evol.nc + fi + #if [ -f "start1D.txt" ]; then + # cp start1D.txt start1D.txt + #fi ./$exePEM > out_runPEM${iPEM} 2>&1 if [ ! -f "restartfi_evol.nc" ]; then # Check if run ended abnormally