MODULE newstart_mod !========================================= ! Author: J. Mauxion (JM in the following) ! Adapted from newstart program !========================================= IMPLICIT NONE PUBLIC :: newstart CONTAINS SUBROUTINE newstart(ngrid,nlayer,nq,pplev,pplay,pphi,& pu_old,pv_old,pt_old,pq_old,& pu_new,pv_new,pt_new,pq_new,ps_new) ! from dynamico USE omp_para, ONLY: is_omp_first_level, ll_begin, & ll_end, ll_endp1 USE dimensions, ONLY: ii_begin, ii_end, jj_begin, jj_end, iim, jjm, & z_rup, z_up, z_lup, z_rdown, z_down, z_ldown, & t_right, t_rup, t_lup, t_ldown, & u_right, u_rup, u_lup, u_left, u_ldown, u_rdown, & swap_dimensions USE abort_mod, ONLY: abort_acc, dynamico_abort ! from LMDZ USE mod_phys_lmdz_para, ONLY: is_parallel, is_sequential,& is_mpi_root,is_omp_root,is_master USE initracer_mod, ONLY: initracer USE phyredem, only: physdem0, physdem1 USE iostart, only: open_startphy USE infotrac, ONLY: infotrac_init, tname USE inichim_newstart_mod, ONLY: inichim_newstart USE tabfi_mod, ONLY: tabfi USE comslope_mod, ONLY: nslope USE ener_mod, ONLY: etot0,ptot0,ztot0,stot0,ang0 USE comvert_mod, ONLY: ap,bp!,pa,preff USE comconst_mod, ONLY: lllm,daysec,dtphys,dtvr,cpp,kappa,rad,omeg,g,r,pi USE comsoil_h, ONLY: nsoilmx, nqsoil,inertiedat,inertiesoil USE surfdat_h, ONLY: z0_default, z0, albedodat USE geometry_mod, ONLY: longitude,latitude,cell_area USE dust_param_mod, ONLY: tauscaling USE tracer_mod, ONLY: noms, mmol,& igcm_dust_number, igcm_dust_mass,& igcm_ccn_number, igcm_ccn_mass,& igcm_h2o_vap, igcm_h2o_ice, igcm_co2,& igcm_hdo_vap, igcm_hdo_ice,& igcm_n2, igcm_ar, igcm_o2, igcm_co,& igcm_o, igcm_h2 USE netcdf USE phyetat0_mod, ONLY: phyetat0 ! inputs: array structure and old dynamical variables ! --------------------------------------------------------------- INTEGER,INTENT(in) :: ngrid ! number of atmospheric columns INTEGER,INTENT(in) :: nlayer ! number of atmospheric layers INTEGER,INTENT(in) :: nq ! number of tracers REAL,INTENT(IN) :: pplev(ngrid,nlayer+1) ! inter-layer pressure (Pa) REAL,INTENT(IN) :: pplay(ngrid,nlayer) ! mid-layer pressure (Pa) REAL,INTENT(IN) :: pphi(ngrid,nlayer) ! geopotential at mid-layer (m2s-2) REAL,INTENT(in) :: pu_old(ngrid,nlayer) ! old zonal wind component (m/s) REAL,INTENT(in) :: pv_old(ngrid,nlayer) ! old meridional wind component (m/s) REAL,INTENT(in) :: pt_old(ngrid,nlayer) ! old temperature (K) REAL,INTENT(in) :: pq_old(ngrid,nlayer,nq) ! old tracers (.../kg_of_air) ! outputs: new dynamical variables ! ---------------------------------------------------------------- REAL,INTENT(out) :: pu_new(ngrid,nlayer) ! new zonal wind (m/s) REAL,INTENT(out) :: pv_new(ngrid,nlayer) ! new meridional wind (m/s) REAL,INTENT(out) :: pt_new(ngrid,nlayer) ! new temperature (K) REAL,INTENT(out) :: pq_new(ngrid,nlayer,nq) ! new tracers (../kg) REAL,INTENT(out) :: ps_new(ngrid) ! new surface pressure (Pa) ! Local variables: ! ---------------- INTEGER :: iq,l,j,i,ij,Lmodif,ierr,islope, tab0 ! some indices and status variables CHARACTER(len=20) :: modif='Default' ! to store which modif to apply CHARACTER(len=3) :: yes ! to check whether input command is correct !!! Physics related variables (mostly used to read with phyetat0) !!! Old for reading, new for storing. !!! Need to ensure coherence between dyn and phy REAL :: tsurf_phy_old(ngrid,nslope) ! surface temperature REAL :: tsoil_phy_old(ngrid,nsoilmx,nslope) ! soil temperature REAL :: albedo_phy_old(ngrid,2,nslope) ! surface albedo REAL :: emis_phy_old(ngrid,nslope) ! surface emissivity REAL :: q2_phy_old(ngrid,nlayer+1) REAL :: qsurf_phy_old(ngrid,nq,nslope) ! tracers on surface REAL :: qsoil_phy_old(ngrid,nsoilmx,nqsoil,nslope) ! tracers in the subsurface REAL :: tauscaling_phy_old(ngrid) ! dust conversion factor REAL :: totcloudfrac_phy_old(ngrid) ! total cloud fraction REAL :: wstar_phy_old(ngrid) ! Max vertical velocity in thermals (m/s) REAL :: watercap_phy_old(ngrid,nslope) ! h2o_ice cover REAL :: perennial_co2ice_phy_old(ngrid,nslope) ! perennial co2 ice(kg/m^2) REAL :: def_slope_phy_old(nslope+1) ! boundaries for bining of the slopes REAL :: def_slope_mean_phy_old(nslope) ! Mean between boundaries REAL :: subslope_dist_phy_old(ngrid,nslope) ! undermesh statistics REAL :: ith_phy_old(ngrid,nsoilmx) ! thermal inertia of the soil REAL :: albsurf_phy_old(ngrid) ! Reduced version of albedo REAL :: ithsurf_phy_old(ngrid) ! Reduce surface version of ith REAL :: tsurf_phy_new(ngrid,nslope) ! surface temperature REAL :: tsoil_phy_new(ngrid,nsoilmx,nslope) ! soil temperature REAL :: albedo_phy_new(ngrid,2,nslope) ! surface albedo REAL :: emis_phy_new(ngrid,nslope) ! surface emissivity REAL :: q2_phy_new(ngrid,nlayer+1) REAL :: qsurf_phy_new(ngrid,nq,nslope) ! tracers on surface REAL :: qsoil_phy_new(ngrid,nsoilmx,nqsoil,nslope) ! tracers in the subsurface REAL :: tauscaling_phy_new(ngrid) ! dust conversion factor REAL :: totcloudfrac_phy_new(ngrid) ! total cloud fraction REAL :: wstar_phy_new(ngrid) ! Max vertical velocity in thermals (m/s) REAL :: watercap_phy_new(ngrid,nslope) ! h2o_ice cover REAL :: perennial_co2ice_phy_new(ngrid,nslope) ! perennial co2 ice(kg/m^2) REAL :: def_slope_phy_new(nslope+1) ! boundaries for bining of the slopes REAL :: def_slope_mean_phy_new(nslope) ! Mean between boundaries REAL :: subslope_dist_phy_new(ngrid,nslope) ! undermesh statistics REAL :: ith_phy_new(ngrid,nsoilmx) ! thermal inertia of the soil REAL :: albsurf_phy_new(ngrid) ! Reduced surface version of albedo REAL :: ithsurf_phy_new(ngrid) ! Reduce surface version of ith REAL :: phis_new(ngrid) ! new geopotential at ground level (both physics and dyn ?) REAL :: ps_old(ngrid) ! Local variable for old ps REAL :: p_new(ngrid,nlayer+1) ! Local variable for new interlayer pressure ! Controle variables (spread over both start files) REAL :: tab_cntrl(100) REAL :: tab_cntrl_bis(100) !!! Variables related to newstart modif !!! ! Variables used to change composition REAL :: latref,lonref,dlatmin,dlonmin ! reference position for composition REAL :: Svmr,Smmr,Smmr_old,Smmr_new,n,Sn REAL :: Mair_old,Mair_new,vmr_old,vmr_new REAL :: coefvmr(nq) ! Correction coefficient when changing composition REAL :: maxq ! Value of most abundant tracer after composition change INTEGER :: iqmax ! Index of most abundant tracer after composition change INTEGER :: iref,jref,ijref ! column index where to perform change ! Others REAL DoverH !D/H ratio REAL :: tsud,albsud,alb_bb,ith_bb,Tiso,val ! to update polar caps INTEGER :: flagthermo, flagh2o ! to ini chemistry LOGICAL :: flagiso ! to set isothermal setup REAL :: ptoto,pcap,patmo,airetot,ptotn,patmn ! Variables to handle pressure changes REAL :: profile(nlayer+1) ! to store an atmospheric profile + surface value REAL :: taufactor !! scalar dust conversion factor REAL :: tauscadyn(ngrid) ! dust conversion factor on the dynamics grid ! Miscellaneous INTEGER :: nid_dyn,nid_fi,nvarid ! variables to handle NETCDF reading CHARACTER(len=50) :: txt ! to store some text REAL :: hour_ini INTEGER :: day_ini REAL :: p_rad,p_omeg,p_g,p_mugaz,p_daysec ! some parameters to read from tabfi REAL :: mugaz ! molar mass of the atmosphere (g/mol) INTEGER :: nlayer_tab ! local copy for tabfi's INTENT(OUT)/(INOUT) nlayer argument ! Initialize new dynamical variables to old ones pu_new(:,:)=pu_old(:,:) pv_new(:,:)=pv_old(:,:) pt_new(:,:)=pt_old(:,:) pq_new(:,:,:)=pq_old(:,:,:) ps_new(:)=pplev(:,1) ps_old(:)=pplev(:,1) ! Load tracer number and names: CALL infotrac_init !----------------------------------------------------------------------- ! Open file(s) to modify (start by default as reso is fixed) !----------------------------------------------------------------------- WRITE(*,*) 'Creating start files from:' WRITE(*,*) './start_icosa.nc and ./startfi.nc' WRITE(*,*) ierr = NF90_OPEN ('start_icosa.nc',NF90_NOWRITE,nid_dyn) IF (ierr.NE.NF90_NOERR) THEN CALL dynamico_abort("newstart: opening start_icosa.nc raised an error.") ENDIF ierr = NF90_OPEN ('startfi.nc', NF90_NOWRITE,nid_fi) IF (ierr.NE.NF90_NOERR) THEN CALL dynamico_abort("newstart: opening startfi.nc raised an error.") ENDIF ! First reading tabfi, no modif. tab0 = 0 Lmodif = 0 !----------------------------------------------------------------------- ! Lecture du tableau des parametres du run (pour la dynamique) !----------------------------------------------------------------------- WRITE(*,*) 'Reading tab_cntrl START_ICOSA' ierr = NF90_INQ_VARID (nid_dyn, "controle", nvarid) ierr = NF90_GET_VAR(nid_dyn, nvarid, tab_cntrl) WRITE(*,*) 'Reading tab_cntrl STARTFI' ierr = NF90_INQ_VARID (nid_fi, "controle", nvarid) ierr = NF90_GET_VAR(nid_fi, nvarid, tab_cntrl_bis) DO i=1,50 tab_cntrl(i+50)=tab_cntrl_bis(i) ENDDO ierr = NF90_CLOSE(nid_dyn) WRITE(*,*) 'Printing tab_cntrl:' DO i=1,100 WRITE(*,*) i,tab_cntrl(i) ENDDO !----------------------------------------------------------------------- ! Initialisation des constantes dynamique !----------------------------------------------------------------------- hour_ini = tab_cntrl(4) kappa = tab_cntrl(9) etot0 = tab_cntrl(12) ptot0 = tab_cntrl(13) ztot0 = tab_cntrl(14) stot0 = tab_cntrl(15) ang0 = tab_cntrl(16) WRITE(*,*) "Newstart: kappa,etot0,ptot0,ztot0albfi" !----------------------------------------------------------------------- ! Lecture du tab_cntrl et initialisation des constantes physiques !----------------------------------------------------------------------- nlayer_tab = nlayer ! Initialize "serial/parallel" related stuff: ! (required because we call tabfi() below, before calling iniphysiq) is_sequential=.false. is_parallel=.true. is_mpi_root=.true. ! correct because we compile parallel but use one core is_omp_root=.true. ! correct because we compile parallel but use one core is_master=.true. ! correct because we compile parallel but use one core CALL open_startphy('startfi.nc') CALL tabfi(nid_fi,Lmodif,tab0,day_ini,nlayer_tab,p_rad,& p_omeg,p_g,p_mugaz,p_daysec,hour_ini) ! JM: I use a local version of tabfi to allow using it despite ! parallelization (at compilation level only, usage is meant ! to be done on one core only). rad = p_rad omeg = p_omeg g = p_g mugaz = p_mugaz daysec = p_daysec !======================================================================= ! Lecture des fichiers (start ou start_archive) !======================================================================= Lmodif=1 ! Lmodif set to 1 to allow modifications in phyetat0 WRITE(*,*) 'Loading file STARTFI' CALL phyetat0 ('startfi.nc',tab0,Lmodif,nsoilmx,ngrid,nlayer,nq,& nqsoil,day_ini,hour_ini,tsurf_phy_old,tsoil_phy_old,albedo_phy_old,emis_phy_old,& q2_phy_old,qsurf_phy_old,qsoil_phy_old,tauscaling_phy_old,totcloudfrac_phy_old,& wstar_phy_old,watercap_phy_old,perennial_co2ice_phy_old,& def_slope_phy_old,def_slope_mean_phy_old,subslope_dist_phy_old) ! Initialize new from old (default is no change) tsurf_phy_new(:,:)=tsurf_phy_old(:,:) tsoil_phy_new(:,:,:)=tsoil_phy_old(:,:,:) albedo_phy_new(:,:,:)=albedo_phy_old(:,:,:) emis_phy_new(:,:)=emis_phy_old(:,:) q2_phy_new(:,:)=q2_phy_old(:,:) qsurf_phy_new(:,:,:)=qsurf_phy_old(:,:,:) qsoil_phy_new(:,:,:,:)=qsoil_phy_old(:,:,:,:) tauscaling_phy_new(:)=tauscaling_phy_old(:) totcloudfrac_phy_new(:)=totcloudfrac_phy_old(:) wstar_phy_new(:)=wstar_phy_old(:) watercap_phy_new(:,:)=watercap_phy_old(:,:) perennial_co2ice_phy_new(:,:)=perennial_co2ice_phy_old(:,:) def_slope_phy_new(:)=def_slope_phy_old(:) def_slope_mean_phy_new(:)=def_slope_mean_phy_old(:) subslope_dist_phy_new(:,:)=subslope_dist_phy_old(:,:) ! Copy albedo (surface) and soil thermal inertia albsurf_phy_old(:)=albedodat(:) albsurf_phy_new(:)=albsurf_phy_old(:) ith_phy_old(:,:)=inertiedat(:,:) ith_phy_new(:,:)=ith_phy_old(:,:) ! Store a surface thermal inertia value ithsurf_phy_old(:)=ith_phy_old(:,1) ithsurf_phy_new(:)=ithsurf_phy_old(:) !======================================================================= ! Initialize tracer names noms(:) and indexes (igcm_co2, igcm_h2o_vap, ...) !======================================================================= CALL initracer(ngrid,nq) !======================================================================= ! MAIN LOOP: looping as long as a change !======================================================================= DO WHILE (modif.ne."") ! JM: possible changes that where in the lonlat version and not ! yet adapted in the present dynamico version are commented out WRITE(*,*) WRITE(*,*) WRITE(*,*) 'List of possible changes :' WRITE(*,*) '~~~~~~~~~~~~~~~~~~~~~~~~~~' WRITE(*,*) !WRITE(*,*) 'flat : no topography ("aquaplanet")' WRITE(*,*) 'bilball : uniform albedo and thermal inertia' WRITE(*,*) 'z0 : set a uniform surface roughness length' WRITE(*,*) 'coldspole : cold subsurface and high albedo at S.Pole' WRITE(*,*) 'qname : change tracer name' WRITE(*,*) 'q=0 : ALL tracer =zero' WRITE(*,*) 'q=factor : change tracer value by a multiplicative factor' WRITE(*,*) 'q=x : give a specific uniform value to one tracer' WRITE(*,*) 'q=profile : specify a profile for a tracer' WRITE(*,*) 'freedust : rescale dust to a true value' WRITE(*,*) 'ini_q : tracers initialization for chemistry and water vapour' WRITE(*,*) 'ini_q-h2o : tracers initialization for chemistry only' WRITE(*,*) 'composition : change atm main composition: CO2,N2,Ar,O2,CO' WRITE(*,*) 'inihdo : initialize HDO' WRITE(*,*) 'ini_h2osurf : reinitialize surface water ice ' WRITE(*,*) 'noglacier : Remove tropical H2O ice if |lat|<45' WRITE(*,*) 'watercapn : H20 ice on permanent N polar cap ' WRITE(*,*) 'watercaps : H20 ice on permanent S polar cap ' WRITE(*,*) 'wetstart : start with a wet atmosphere' WRITE(*,*) 'isotherm : Isothermal Temperatures, wind set to zero' WRITE(*,*) 'co2ice=0 : remove CO2 polar cap i.e. qsurf(co2)=0 ' !WRITE(*,*) 'ptot : change total pressure' !WRITE(*,*) 'therm_ini_s : set soil thermal inertia to reference surface values' !WRITE(*,*) 'subsoilice_n : put deep underground ice layer in northern hemisphere' !WRITE(*,*) 'subsoilice_s : put deep underground ice layer in southern hemisphere' !WRITE(*,*) 'mons_ice : put underground ice layer according to MONS derived data' !WRITE(*,*) 'nslope : Change the number of subgrid scale slope' WRITE(*,*) WRITE(*,*) 'Change to perform ?' WRITE(*,*) ' (enter keyword or return to end)' WRITE(*,*) READ(*,fmt='(a20)') modif WRITE(*,*) WRITE(*,*) trim(modif) , ' : ' ! bilball : albedo, inertie thermique unIForme ! -------------------------------------------- IF (trim(modif).eq.'bilball') THEN WRITE(*,*) 'Constant albedo and thermal inertia:' WRITE(*,*) 'New value for albedo (ex: 0.25) ?' READ(*,*,iostat=ierr) alb_bb DO WHILE (ierr.ne.0) READ(*,*,iostat=ierr) alb_bb ENDDO WRITE(*,*) 'Uniform albedo (new value):',alb_bb WRITE(*,*) 'New value for thermal inertia (eg: 247) ?' READ(*,*,iostat=ierr) ith_bb DO WHILE (ierr.ne.0) READ(*,*,iostat=ierr) ith_bb ENDDO WRITE(*,*) 'Uniform thermal inertia (new value):',ith_bb albedo_phy_new(:,:,:)=alb_bb albsurf_phy_new(:)=alb_bb ith_phy_new(:,:)=ith_bb ithsurf_phy_new(:)=ith_bb ! Also reset surface roughness length to default value WRITE(*,*) 'Surface roughness length set to:',z0_default z0(:)=z0_default ! z0 : set surface roughness length to a constant value ! ----------------------------------------------------- ELSE IF (trim(modif).eq.'z0') THEN WRITE(*,*) 'Set a uniform surface roughness length:' WRITE(*,*) 'Value for z0_default (ex: ',z0_default,') ?' READ(*,*,iostat=ierr) z0_default DO WHILE (ierr.ne.0) READ(*,*,iostat=ierr) z0_default ENDDO z0(:)=z0_default ! coldspole : always cold south polar sub-soil cap ! ------------------------------------------------ ELSE IF (trim(modif).eq.'coldspole') THEN WRITE(*,*) 'Set always cold south polar sub-soil cap:' WRITE(*,*)'New value for the subsurface temperature'//& ' beneath the permanent southern polar cap ? (eg: 141 K)' READ(*,*,iostat=ierr) tsud DO WHILE (ierr.ne.0) READ(*,*,iostat=ierr) tsud ENDDO WRITE(*,*) 'New value of the subsurface temperature:',tsud DO l=2,nsoilmx tsoil_phy_new(:,l,:) = tsud ENDDO WRITE(*,*)'New value for the albedo '//& 'of the permanent southern polar cap ? (eg: 0.75)' READ(*,*,iostat=ierr) albsud DO WHILE (ierr.ne.0) READ(*,*,iostat=ierr) albsud ENDDO albedo_phy_new(:,:,:)=albsud albsurf_phy_new(:)=albsud WRITE(*,*) 'ig=',ngrid,'albedo perennial cap=',albsurf_phy_new(ngrid) ! ptot : modification of the total pressure=ice+current atmosphere ! ------------------------------------------------------------------- ELSE IF (trim(modif).eq.'ptot') THEN ! calcul de la pression totale glace + atm actuelle patmo=0. airetot=0. pcap=0. DO ij=1,ngrid patmo=patmo+ps_old(ij)*cell_area(ij) airetot=airetot+cell_area(ij) DO islope=1,nslope pcap = pcap+cell_area(ij)*qsurf_phy_old(ij,igcm_co2,islope)*g*& subslope_dist_phy_old(ij,islope)/cos(def_slope_mean_phy_old(islope)*pi/180.) ENDDO ENDDO ptoto = pcap + patmo WRITE(*,*) 'Current total pressure at surface (co2 ice + atm):',& ptoto/airetot WRITE(*,*) 'New value ?' READ(*,*) ptotn ptotn=ptotn*airetot patmn=ptotn-pcap WRITE(*,*) 'ptoto,patmo,ptotn,patmn:' WRITE(*,*) ptoto,patmo,ptotn,patmn WRITE(*,*)'Mult. factor for pressure (atm only):', patmn/patmo ps_new(:)=ps_old(:)*patmn/patmo ! correction pour la conservation des traceurs yes=' ' DO WHILE((yes.ne.'y').and.(yes.ne.'n')) WRITE(*,*) 'Do you wish to conserve tracer total mass (y)'//& ' or tracer mixing ratio (n) ?' READ(*,fmt='(a)') yes ENDDO IF (yes.eq.'y') THEN WRITE(*,*) 'OK: conservation of tracer total mass' pq_new(:,:,:)=pq_old(:,:,:)*patmo/patmn ELSE WRITE(*,*) 'OK: conservation of tracer mixing ratio' ENDIF ! qname : change tracer name ! -------------------------- ELSE IF (trim(modif).eq.'qname') THEN yes='y' DO WHILE (yes.eq.'y') WRITE(*,*) 'Which tracer name do you want to change ?' DO iq=1,nq WRITE(*,'(i3,a3,a20)') iq,':',trim(tname(iq)) ENDDO WRITE(*,'(a35,i3)') '(enter tracer number; between 1 and ',nq WRITE(*,*) ' or any other value to quit this option)' READ(*,*) iq IF ((iq.ge.1).and.(iq.le.nq)) THEN WRITE(*,*) 'Change tracer name ',trim(tname(iq)),' to ?' READ(*,*) txt tname(iq)=txt WRITE(*,*) 'Do you want to change another tracer name (y/n)?' READ(*,'(a)') yes ELSE ! inapropiate value of iq; quit this option yes='n' ENDIF ! of if ((iq.ge.1).and.(iq.le.nqtot)) ENDDO ! of do while (yes.ne.'y') ! q=0 : set tracers to zero ! ------------------------- ELSE IF (trim(modif).eq.'q=0') THEN ! mise a 0 des q (traceurs) WRITE(*,*) 'Tracers set to 0 (1.E-30 in fact)' pq_new(:,:,:)=0. ! set surface tracers to zero qsurf_phy_new(:,:,:)=0. ! q=factor : change value of tracer by a multiplicative factor ! ------------------------------------------------------------ ELSE IF (trim(modif).eq.'q=factor') THEN yes='y' DO WHILE (yes.eq.'y') WRITE(*,*) 'Which tracer do you want to modify ?' DO iq=1,nq WRITE(*,*) iq,':',trim(tname(iq)) ENDDO WRITE(*,*) '(choose between 1 and ',nq,')' READ(*,*) iq IF ((iq.ge.1).and.(iq.le.nq)) THEN WRITE(*,*) "Factor to multiply current mixing ratio by ?" READ(*,*) val pq_new(:,:,iq)=pq_old(:,:,iq)*val qsurf_phy_new(:,iq,:)=qsurf_phy_old(:,iq,:)*val WRITE(*,*) 'Do you want to change another tracer (y/n)?' READ(*,'(a)') yes ELSE ! inapropiate value of iq; quit this option yes='n' ENDIF ENDDO ! q=x : initialise tracer manually ! -------------------------------- ELSE IF (trim(modif).eq.'q=x') THEN yes='y' DO WHILE (yes.eq.'y') WRITE(*,*) 'Which tracer do you want to modify ?' DO iq=1,nq WRITE(*,*) iq,':',trim(tname(iq)) ENDDO WRITE(*,*) '(choose between 1 and ',nq,')' READ(*,*) iq IF ((iq.ge.1).and.(iq.le.nq)) THEN WRITE(*,*)'Mixing ratio of tracer ',trim(tname(iq))//& ' ? (kg/kg)' READ(*,*) val pq_new(:,:,iq)=val WRITE(*,*) 'SURFACE value of tracer ',trim(tname(iq))//& ' ? (kg/m2)' READ(*,*) val qsurf_phy_new(:,iq,:)=val WRITE(*,*) 'Do you want to change another tracer (y/n)?' READ(*,'(a)') yes ELSE ! inapropiate value of iq; quit this option yes='n' ENDIF ENDDO ! q=profile : initialize tracer with a given profile ! -------------------------------------------------- ELSE IF (trim(modif).eq.'q=profile') THEN WRITE(*,*) 'Tracer profile will be sought in ASCII file' WRITE(*,*) "'profile_tracer' where 'tracer' is tracer name" WRITE(*,*) "(one value per line in file; starting with" WRITE(*,*) "surface value, the 1st atmospheric layer" WRITE(*,*) "followed by 2nd, etc. up to top of atmosphere)" yes='y' DO WHILE (yes.eq.'y') WRITE(*,*) 'Which tracer do you want to set?' DO iq=1,nq WRITE(*,*)iq,':',trim(tname(iq)) ENDDO WRITE(*,*) '(choose between 1 and ',nq,')' READ(*,*) iq IF ((iq.ge.1).and.(iq.le.nq)) THEN ! look for input file 'profile_tracer' txt="profile_"//trim(tname(iq)) OPEN(41,file=trim(txt),status='old',form='formatted',iostat=ierr) IF (ierr.eq.0) THEN ! OK, found file 'profile_...', load the profile DO l=1,nlayer+1 READ(41,*,iostat=ierr) profile(l) IF (ierr.ne.0) THEN ! something went wrong CALL dynamico_abort('newstart: something went wrong while reading the profile.') ENDIF ENDDO IF (ierr.eq.0) THEN ! initialize tracer values qsurf_phy_new(:,iq,:)=profile(1) DO l=1,nlayer pq_new(:,l,iq)=profile(l+1) ENDDO WRITE(*,*)'OK, tracer ',trim(tname(iq)),& ' initialized ','using values from file ',trim(txt) ELSE WRITE(*,*)'Problem reading file ',trim(txt),' !' WRITE(*,*)'No modifications to tracer ',trim(tname(iq)) ENDIF ELSE WRITE(*,*)'Could not find file ',trim(txt),' !' WRITE(*,*)'No modifications to tracer ',trim(tname(iq)) ENDIF WRITE(*,*) 'Do you want to set another tracer profile (y/n)?' READ(*,'(a)') yes ELSE ! inapropiate value of iq; quit this option yes='n' ENDIF ENDDO ! freedust: convert dust from virtual to true values ! -------------------------------------------------- ELSE IF (trim(modif).eq.'freedust') THEN IF (minval(tauscaling_phy_old).lt.0) THEN WRITE(*,*) 'WARNING conversion factor negative' WRITE(*,*) 'This is probably because it was not present'//& 'in the file' WRITE(*,*) 'A constant conversion is used instead.' taufactor=1.e-3 ELSE taufactor=1. ENDIF IF (igcm_dust_number.ne.0) pq_new(:,:,igcm_dust_number)=pq_old(:,:,igcm_dust_number)*taufactor IF (igcm_dust_mass.ne.0) pq_new(:,:,igcm_dust_mass)=pq_old(:,:,igcm_dust_mass)*taufactor IF (igcm_ccn_number.ne.0) pq_new(:,:,igcm_ccn_number)=pq_old(:,:,igcm_ccn_number)*taufactor IF (igcm_ccn_mass.ne.0) pq_new(:,:,igcm_ccn_mass)=pq_old(:,:,igcm_ccn_mass)*taufactor tauscaling_phy_new(:) = 1. WRITE(*,*) 'Done rescaling to true value' ! ini_q : Initialize tracers for chemistry ! ----------------------------------------------- ELSE IF (trim(modif).eq.'ini_q') THEN flagh2o=1 flagthermo=0 yes=' ' ! For more than 32 layers, possible to initiate thermosphere only IF (nlayer.gt.32) THEN DO WHILE ((yes.ne.'y').and.(yes.ne.'n')) WRITE(*,*)'Initialisation for thermosphere only? (y/n)' READ(*,fmt='(a)') yes IF (yes.eq.'y') THEN flagthermo=1 ELSE flagthermo=0 ENDIF ENDDO ENDIF CALL inichim_newstart(ngrid,nq,pq_new,qsurf_phy_new,ps_old,flagh2o,flagthermo) WRITE(*,*) 'inichim_newstart: chemical species and '//& 'water vapour initialised' ! ini_q-h2o : as above except for the water vapour tracer ! ------------------------------------------------------ ELSE IF (trim(modif).eq.'ini_q-h2o') THEN flagh2o = 0 flagthermo = 0 yes=' ' ! for more than 32 layers, possible to initiate thermosphere only IF (nlayer.gt.32) THEN DO WHILE ((yes.ne.'y').and.(yes.ne.'n')) WRITE(*,*)'Initialisation for thermosphere only? (y/n)' READ(*,fmt='(a)') yes IF (yes.eq.'y') THEN flagthermo=1 ELSE flagthermo=0 ENDIF ENDDO ENDIF CALL inichim_newstart(ngrid,nq,pq_new,qsurf_phy_new,ps_old,& flagh2o, flagthermo) WRITE(*,*) 'inichim_newstart: chemical species initialised'//& '(except water vapour)' ! inihdo : initialize HDO with user D/H value ! -------------------------------------------------------- ELSE IF (trim(modif).eq.'inihdo') THEN ! check that there is indeed a water vapor tracer IF (igcm_h2o_vap.eq.0) THEN CALL dynamico_abort("newstart : no water vapour tracer! Can't use this option.") ENDIF WRITE(*,*)'Input D/H ratio (in SMOW)?' WRITE(*,*)'If value is <0 THEN HDO=H2O' READ(*,*, iostat=ierr) DoverH DO WHILE (ierr.ne.0) READ(*,*, iostat=ierr) DoverH ENDDO DoverH = DoverH * 2 * 155.76e-6 IF (DoverH.lt.0) THEN DoverH = 1. ENDIF ! D/H (SMOW) = 155.76e-6 so HDO/H2O is twice that qsurf_phy_new(:,igcm_h2o_ice,:)=max(0.,qsurf_phy_old(:,igcm_h2o_ice,:)) ! Update the hDO tracers pq_new(:,:,igcm_hdo_vap)=pq_old(:,:,igcm_h2o_vap)*DoverH pq_new(:,:,igcm_hdo_ice)=pq_old(:,:,igcm_h2o_ice)*DoverH qsurf_phy_new(:,igcm_hdo_ice,:)=qsurf_phy_old(:,igcm_h2o_ice,:)*DoverH ! composition : change main composition: CO2,N2,Ar,O2,CO (FF 03/2014) ! -------------------------------------------------------- ELSE IF (trim(modif).eq.'composition') THEN WRITE(*,*) "Lat (degN) lon (degE) of the reference site ?" WRITE(*,*) "e.g. MSL : -4.5 137. " READ(*,*,iostat=ierr) latref, lonref DO WHILE (ierr.ne.0) READ(*,*,iostat=ierr) latref, lonref ENDDO ! select GCM point close to reference site dlonmin=90. dlatmin=45. iref=1 jref=1 DO ij=1,ngrid IF (abs(longitude(ij)*180./pi-lonref).lt.dlonmin) THEN iref=ij dlonmin=abs(longitude(ij)*180./pi-lonref) ENDIF IF (abs(latitude(ij)*180./pi-latref).lt.dlatmin) THEN jref=ij dlonmin=abs(longitude(ij)*180./pi-lonref) dlatmin=abs(latitude(ij)*180./pi-latref) ENDIF ENDDO IF (iref.ne.jref) THEN WRITE(*,*) "INFO: closest lon and lat are different." WRITE(*,*) "iref=",iref,", jref=",jref WRITE(*,*) "lon=",longitude(iref)*180./pi,", lat=",& latitude(jref)*180./pi WRITE(*,*) "dlonmin=", dlonmin,", dlatmin=",dlatmin IF (dlonmin.le.dlatmin) THEN WRITE(*,*) "Keeping iref." ijref=iref ELSE WRITE(*,*) "Keeping jref." ijref=jref ENDIF ELSE ijref=iref ENDIF WRITE(*,*) "In GCM : lat=",latitude(ijref)*180./pi WRITE(*,*) "In GCM : lon=",longitude(ijref)*180./pi WRITE(*,*) ! compute air molar mass at reference site Smmr=0. Sn=0. WRITE(*,*) 'igcm_co2=',igcm_co2 WRITE(*,*) 'igcm_n2=',igcm_n2 WRITE(*,*) 'igcm_ar=',igcm_ar WRITE(*,*) 'igcm_o2=',igcm_o2 WRITE(*,*) 'igcm_co=',igcm_co WRITE(*,*) noms DO iq=1,nq IF ((iq.eq.igcm_co2).or.(iq.eq.igcm_n2).or.& (iq.eq.igcm_ar).or.(iq.eq.igcm_o2).or.& (iq.eq.igcm_co)) THEN Smmr=Smmr+pq_old(ijref,1,iq) Sn=Sn+pq_old(ijref,1,iq)/mmol(iq) ENDIF ENDDO ! Special case : implicit non-co2 gases ! JN 11/2019 IF ((igcm_n2.eq.0).or.(igcm_ar.eq.0)) THEN WRITE(*,*) "Warning: non-co2 gases are implicit." WRITE(*,*) "At reference site:" WRITE(*,*) "q=", pq_old(ijref,1,igcm_co2) WRITE(*,*) "Sum of mass mix. ratio (ie MMR(co2))=",Smmr Mair_old = 44.0*Smmr+33.87226017157708*(1-Smmr) ! 33.87226017157708 is the molar mass of non-co2 ! atmosphere measured by MSL at Ls ~184 ELSE ! Assume co2/n2/ar/o2/co are available Mair_old=(pq_old(ijref,1,igcm_co2)*mmol(igcm_co2)& +pq_old(ijref,1,igcm_n2)*mmol(igcm_n2)& +pq_old(ijref,1,igcm_ar)*mmol(igcm_ar)& +pq_old(ijref,1,igcm_o2)*mmol(igcm_o2)& +pq_old(ijref,1,igcm_co)*mmol(igcm_co))/Smmr ENDIF WRITE(*,*) "Air molar mass (g/mol) at reference site=",Mair_old ! Ask for new volume mixing ratio at reference site Svmr =0. Sn=0. DO iq=1,nq coefvmr(iq)=1. IF ((iq.eq.igcm_n2).or.(iq.eq.igcm_ar).or.& (iq.eq.igcm_o2).or.(iq.eq.igcm_co)) THEN vmr_old=pq_old(ijref,1,iq)*Mair_old/mmol(iq) WRITE(*,*) "Previous vmr("//trim(tname(iq))//")=", vmr_old IF (iq.eq.igcm_n2) THEN WRITE(*,*) "New vmr(n2)? (MSL: 2.8e-02 at Ls~180, Trainer et al. 2019)" ENDIF IF (iq.eq.igcm_ar) THEN WRITE(*,*) "New vmr(ar)? (MSL: 2.1e-02 at Ls~180)" ENDIF IF (iq.eq.igcm_o2) THEN WRITE(*,*) "New vmr(o2)? (MSL: 1.7e-03 at Ls~180)" ENDIF IF (iq.eq.igcm_co) THEN WRITE(*,*) "New vmr(co)? (ACS: 1e-03 at Ls~180)" ENDIF READ(*,*,iostat=ierr) vmr_new DO WHILE (ierr.ne.0) READ(*,*,iostat=ierr) vmr_new ENDDO WRITE(*,*) "New vmr("//trim(tname(iq))//")=",vmr_new WRITE(*,*) coefvmr(iq)=vmr_new/vmr_old Svmr=Svmr+vmr_new Sn=Sn+vmr_new*mmol(iq) ENDIF ENDDO ! of do iq=1,nqtot IF ((igcm_n2.eq.0).or.(igcm_ar.eq.0)) THEN WRITE(*,*) "Warning: non-co2 gases are implicit." vmr_old=pq_old(ijref,1,igcm_co2)*Mair_old/mmol(igcm_co2) WRITE(*,*) "Previous vmr(co2)=", vmr_old WRITE(*,*) "New vmr(co2) ? (MSL: 0.947 at Ls~180) Trainer et al. 2019)" READ(*,*,iostat=ierr) vmr_new DO WHILE (ierr.ne.0) READ(*,*,iostat=ierr) vmr_new ENDDO coefvmr(igcm_co2)=vmr_new/vmr_old Svmr=Svmr+vmr_new Sn=vmr_new*mmol(igcm_co2)+(1-vmr_new)*33.87226017157708 ! Molar mass of non-co2 atm from MSL ENDIF ! Estimation of new Air molar mass at reference site (assuming vmr_co2 = 1-Svmr) Mair_new=Sn+(1-Svmr)*mmol(igcm_co2) ! Estimation of new Air molar mass when non-co2 gases are implicit IF ((igcm_n2.eq.0).or.(igcm_ar.eq.0)) THEN Mair_new=vmr_new*mmol(igcm_co2)+(1-vmr_new)& *33.87226017157708 ! Molar mass of non-co2 atm from MSL WRITE(*,*) "We consider non-co2 gases vmr measured from Curiosity" ENDIF WRITE(*,*) "NEW Air molar mass (g/mol) at reference site=",Mair_new ! Compute mass mixing ratio changes DO iq=1,nq IF ((iq.eq.igcm_co2).or.(iq.eq.igcm_n2).or.(iq.eq.igcm_ar)& .or.(iq.eq.igcm_o2).or.(iq.eq.igcm_co)) THEN WRITE(*,*) "Everywhere mmr("//trim(tname(iq))//& ") is multiplied by ",coefvmr(iq)*Mair_old/Mair_new ENDIF ENDDO ! Recompute mass mixing ratios everywhere, and adjust mmr of most abundant species ! to keep sum of mmr constant. DO l=1,nlayer DO ij=1,ngrid Smmr_old = 0. Smmr_new = 0. DO iq=1,nq IF ((iq.eq.igcm_co2).or.(iq.eq.igcm_n2).or.(iq.eq.igcm_ar)& .or.(iq.eq.igcm_o2).or.(iq.eq.igcm_co)& .or.(iq.eq.igcm_o).or.(iq.eq.igcm_h2)) THEN Smmr_old=Smmr_old+pq_old(ij,l,iq) ! sum of old mmr pq_new(ij,l,iq)=pq_old(ij,l,iq)*coefvmr(iq)*Mair_old/Mair_new Smmr_new=Smmr_new+pq_new(ij,l,iq) ! sum of new mmr ENDIF ENDDO iqmax=0 maxq=0. DO iq=1,nq IF ((iq.eq.igcm_co2).or.(iq.eq.igcm_n2).or.(iq.eq.igcm_ar)& .or.(iq.eq.igcm_o2).or.(iq.eq.igcm_co)& .or.(iq.eq.igcm_o).or.(iq.eq. igcm_h2)) THEN IF (pq_new(ij,l,iq).gt.maxq) THEN maxq=pq_new(ij,l,iq) iqmax=iq ENDIF ENDIF ENDDO pq_new(ij,l,iqmax)=pq_new(ij,l,iqmax)+Smmr_old-Smmr_new ENDDO ENDDO WRITE(*,*) "The most abundant species is modified everywhere to keep "//& "sum of mmr constant" WRITE(*,*) 'At reference site vmr(CO2)=',pq_new(ijref,1,igcm_co2)*Mair_new/mmol(igcm_co2) WRITE(*,*) "Compared to MSL observation: vmr(CO2)= 0.947 at Ls=180" Sn = pq_new(ijref,1,igcm_co2)*Mair_new/mmol(igcm_co2)& + pq_new(ijref,1,igcm_n2)*Mair_new/mmol(igcm_n2)& + pq_new(ijref,1,igcm_ar)*Mair_new/mmol(igcm_ar)& + pq_new(ijref,1,igcm_o2)*Mair_new/mmol(igcm_o2)& + pq_new(ijref,1,igcm_co)*Mair_new/mmol(igcm_co) WRITE(*,*) 'Sum of volume mixing ratios=', Sn ! wetstart : wet atmosphere with a north to south gradient ! -------------------------------------------------------- ELSE IF (trim(modif).eq.'wetstart') THEN ! check that there is indeed a water vapor tracer IF (igcm_h2o_vap.eq.0) THEN CALL dynamico_abort("newstart: no water vapour tracer ! Can't use this option.") ENDIF DO ij=1,ngrid pq_new(ij,:,igcm_h2o_vap)=150.e-6*(latitude(ij)+pi/2.)/pi ENDDO WRITE(*,*) 'Water mass mixing ratio at north pole=',150.e-6 WRITE(*,*) '---------------------------south pole=',0. ! ini_h2osurf : reinitialize surface water ice ! -------------------------------------------------- ELSE IF (trim(modif).eq.'ini_h2osurf') THEN WRITE(*,*)'max surface ice left?(e.g. 0.2 kg/m2=200microns)' READ(*,*,iostat=ierr) val DO WHILE (ierr.ne.0) READ(*,*,iostat=ierr) val ENDDO WRITE(*,*)'Also set negative values of surf ice to 0' qsurf_phy_new(:,igcm_h2o_ice,:)=min(val,qsurf_phy_old(:,igcm_h2o_ice,:)) qsurf_phy_new(:,igcm_h2o_ice,:)=max(0.,qsurf_phy_old(:,igcm_h2o_ice,:)) ! noglacier : remove tropical water ice (to initialize high res sim) ! -------------------------------------------------- ELSE IF (trim(modif).eq.'noglacier') THEN WRITE(*,*) 'OK: remove surface ice for |lat|<45' DO ij=1,ngrid IF (abs(latitude(ij)*180./pi).lt.45.) THEN qsurf_phy_new(ij,igcm_h2o_ice,:)=0. ENDIF ENDDO ! watercapn : H20 ice on permanent northern cap ! -------------------------------------------------- ELSE IF (trim(modif).eq.'watercapn') THEN DO ij=1,ngrid IF (latitude(ij)*180./pi.gt.80.) THEN qsurf_phy_new(ij,igcm_h2o_ice,:)=1.e5 WRITE(*,*) '(ig,islope)=(',ij,',:)', 'H2O ice mass (kg/m2)= ',& qsurf_phy_new(ij,igcm_h2o_ice,:) WRITE(*,*)'==> Ice mesh North boundary (deg)= ',latitude(ij)*180./pi ENDIF ENDDO ! watercaps : H20 ice on permanent southern cap ! ------------------------------------------------- ELSE IF (trim(modif).eq.'watercaps') THEN DO ij=1,ngrid IF (latitude(ij)*180./pi.lt.-80.) THEN qsurf_phy_new(ij,igcm_h2o_ice,:)=1.e5 WRITE(*,*) '(ig,islope)=(',ij,',:)', 'H2O ice mass (kg/m2)= ',& qsurf_phy_new(ij,igcm_h2o_ice,:) WRITE(*,*)'==> Ice mesh South boundary (deg)= ',latitude(ij)*180./pi ENDIF ENDDO ! isotherm : isothermal temperatures and no winds ! ------------------------------------------------ ELSE IF (trim(modif).eq.'isotherm') THEN WRITE(*,*) 'Isothermal temperature of the atmosphere, surface and subsurface' WRITE(*,*) 'Value of this temperature ?' READ(*,*,iostat=ierr) Tiso DO WHILE (ierr.ne.0) READ(*,*,iostat=ierr) Tiso ENDDO tsurf_phy_new(:,:)=Tiso tsoil_phy_new(:,:,:)=Tiso pt_new(:,:)=Tiso flagiso=.true. pu_new(:,:)=0 pv_new(:,:)=0 q2_phy_new(:,:)=0 ! co2ice=0 : remove CO2 polar ice caps' ! ------------------------------------------------ ELSE IF (trim(modif).eq.'co2ice=0') THEN qsurf_phy_new(:,igcm_co2,:)=0 DO ij=1,ngrid emis_phy_new(ij,:)=emis_phy_old(ngrid/2,:) ! JM: raw copy from original but I am not sure why END DO ENDIF ! of if modif=... ENDDO ! of do while modif.ne."" !======================================================================= ! MAIN LOOP: END MAIN LOOP !======================================================================= !======================================================================= ! WRITING RESTARTFI !======================================================================= CALL physdem0("restartfi.nc",longitude,latitude,& nsoilmx,ngrid,nlayer,& nq,dtphys,REAL(day_ini),hour_ini,cell_area,& albsurf_phy_new,ith_phy_new,def_slope_phy_new,& subslope_dist_phy_new) CALL physdem1("restartfi.nc",nsoilmx,ngrid,nlayer,nq,nqsoil,& dtphys,hour_ini,& tsurf_phy_new,tsoil_phy_new,inertiesoil,albedo_phy_new,& emis_phy_new,q2_phy_new,qsurf_phy_new,qsoil_phy_new,& tauscaling_phy_new,totcloudfrac_phy_new,wstar_phy_new,& watercap_phy_new,perennial_co2ice_phy_new) WRITE(*,*) "newstart: All is well that ends well." END SUBROUTINE newstart END MODULE newstart_mod