C====================================================================== PROGRAM newstart c======================================================================= c c c Auteur: Christophe Hourdin/Francois Forget/Yann Wanherdrick c ------ c Adapted to Pluto: Tanguy Bertrand c Derniere modif : 06/2024 c c Objet: Create or modify the initial state for the LMD Mars GCM c ----- (fichiers NetCDF start et startfi) c c======================================================================= use mod_phys_lmdz_para, only: is_parallel, is_sequential, & is_mpi_root, is_omp_root, & is_master use infotrac, only: infotrac_init, nqtot, tname USE tracer_h, ONLY: igcm_n2,igcm_ch4_ice,igcm_co_ice,igcm_haze, & igcm_prec_haze,igcm_co_gas,igcm_ch4_gas,noms USE comsoil_h, ONLY: nsoilmx, layer, mlayer, inertiedat USE surfdat_h, ONLY: phisfi, albedodat, & zmea, zstd, zsig, zgam, zthe use datafile_mod, only: datadir, surfdir use ioipsl_getin_p_mod, only: getin_p use control_mod, only: day_step, iphysiq, anneeref, planet_type use phyredem, only: physdem0, physdem1 use iostart, only: open_startphy use filtreg_mod, only: inifilr USE mod_const_mpi, ONLY: COMM_LMDZ USE comvert_mod, ONLY: ap,bp,aps,bps,pa,preff USE comconst_mod, ONLY: lllm,daysec,dtvr,dtphys,cpp,kappa, . rad,omeg,g,r,pi USE serre_mod, ONLY: alphax USE temps_mod, ONLY: day_ini USE ener_mod, ONLY: etot0,ptot0,ztot0,stot0,ang0 use tabfi_mod, only: tabfi use dimphy, only: init_dimphy use iniphysiq_mod, only: iniphysiq use phys_state_var_mod, only: phys_state_var_init use phyetat0_mod, only: phyetat0 use exner_hyb_m, only: exner_hyb use geometry_mod, only: longitude, ! longitudes (rad) & latitude, ! latitudes (rad) & cell_area ! physics grid area (m2) implicit none include "dimensions.h" integer, parameter :: ngridmx = (2+(jjm-1)*iim - 1/jjm) include "paramet.h" include "comgeom2.h" include "comdissnew.h" include "netcdf.inc" c======================================================================= c Declarations c======================================================================= c Variables dimension du fichier "start_archive" c------------------------------------ CHARACTER relief*3 c Variables pour les lectures NetCDF des fichiers "start_archive" c-------------------------------------------------- INTEGER nid_dyn, nid_fi,nid,nvarid INTEGER length parameter (length = 100) INTEGER tab0 INTEGER NB_ETATMAX parameter (NB_ETATMAX = 100) REAL date REAL p_rad,p_omeg,p_g,p_cpp,p_mugaz,p_daysec c Variable histoire c------------------ REAL vcov(iip1,jjm,llm),ucov(iip1,jjp1,llm) ! vents covariants REAL phis(iip1,jjp1) REAL,ALLOCATABLE :: q(:,:,:,:) ! champs advectes c autre variables dynamique nouvelle grille c------------------------------------------ REAL pks(iip1,jjp1) REAL w(iip1,jjp1,llm+1) REAL pbaru(ip1jmp1,llm),pbarv(ip1jm,llm) REAL phi(iip1,jjp1,llm) integer klatdat,klongdat PARAMETER (klatdat=180,klongdat=360) c Physique sur grille scalaire c---------------------------- real zmeaS(iip1,jjp1),zstdS(iip1,jjp1) real zsigS(iip1,jjp1),zgamS(iip1,jjp1),ztheS(iip1,jjp1) c variable physique c------------------ REAL tsurf(ngridmx) ! surface temperature REAL,ALLOCATABLE :: tsoil(:,:) ! soil temperature REAL emis(ngridmx) ! surface emissivity real emisread ! added by RW REAL,ALLOCATABLE :: qsurf(:,:) REAL,ALLOCATABLE :: qsurf_tmp(:,:) REAL q2(ngridmx,llm+1) real alb(iip1,jjp1),albfi(ngridmx) ! albedos real,ALLOCATABLE :: ith(:,:,:),ithfi(:,:) ! thermal inertia (3D) real surfith(iip1,jjp1),surfithfi(ngridmx) ! surface thermal inertia (2D) REAL :: latfi(ngridmx),lonfi(ngridmx),airefi(ngridmx) INTEGER i,j,l,isoil,ig,idum,k real mugaz ! molar mass of the atmosphere integer ierr c Variables on the new grid along scalar points c------------------------------------------------------ REAL t(iip1,jjp1,llm) REAL tset(iip1,jjp1,llm) real phisold_newgrid(iip1,jjp1) real phisold(iip1,jjp1) REAL :: teta(iip1, jjp1, llm) REAL :: pk(iip1,jjp1,llm) REAL :: pkf(iip1,jjp1,llm) REAL :: ps(iip1, jjp1) REAL :: masse(iip1,jjp1,llm) REAL :: xpn,xps,xppn(iim),xpps(iim) REAL :: p3d(iip1, jjp1, llm+1) c Variable de l'ancienne grille c------------------------------ real time real tab_cntrl(100) real tab_cntrl_bis(100) c variables diverses c------------------- real choix_1,pp,choice character*80 fichnom character*250 filestring integer Lmodif,iq character modif*20 real z_reel(iip1,jjp1) real tsud,albsud,alb_bb,ith_bb,Tiso,Tabove real ptoto,pcap,patm,airetot,ptotn,patmn,psea character*1 yes logical :: flagtset=.false. , flagps0=.false. real val, val1, val2, val3, val4, dist, ref ! to store temporary variables real val5, val6, val7, val8, val9, val10,val11, val12 ! to store temporary variables real :: iceith=2000 ! thermal inertia of subterranean ice integer :: iref,jref,iref1,jref1,iref2,jref2 integer :: igref,igref1,igref2,igref3 INTEGER :: itau character(len=20) :: txt ! to store some text character(len=50) :: surfacefile ! "surface.nc" (or equivalent file) character(len=150) :: longline integer :: count real :: profile(llm+1) ! to store an atmospheric profile + surface value ! added by BC for equilibrium temperature startup real teque ! added by RW for nuketharsis real fact1 real fact2 c Special Pluto Map from Lellouch & Grundy c------------------------------------------ integer,parameter :: im_plu=360 ! from topo integer,parameter :: jm_plu=180 real latv_plu(jm_plu+1),lonu_plu(im_plu+1) real map_pluto_dat(im_plu,jm_plu+1) real N2_pluto_dat(im_plu,jm_plu+1) real CH4_pluto_dat(im_plu,jm_plu+1) real CO_pluto_dat(im_plu,jm_plu+1) real alb_dat(im_plu,jm_plu+1) real N2_pluto_rein(im_plu+1,jm_plu+1) real CH4_pluto_rein(im_plu+1,jm_plu+1) real CO_pluto_rein(im_plu+1,jm_plu+1) real alb_rein(im_plu+1,jm_plu+1) real N2_pluto_gcm(iip1,jjp1) real CH4_pluto_gcm(iip1,jjp1) real CO_pluto_gcm(iip1,jjp1) real alb_gcm(iip1,jjp1) c Special Topo Map mountain real lat0, lon0 integer ig0 real dist_pluto real radm,height, phisprev, temp real preffnew,panew c Special copy of terrain real actualmin,angle integer array_ind(ngridmx) real array_dist(ngridmx) real array_angle(ngridmx) c sortie visu pour les champs dynamiques CALL conf_gcm( 99, .TRUE. ) cpp = 0. preff = 2. pa = 0.2 ! to ensure disaster rather than minor error if we don`t ! make deliberate choice of these values elsewhere. planet_type="generic" ! initialize "serial/parallel" related stuff ! (required because we call tabfi() below, before calling iniphysiq) is_sequential=.true. is_parallel=.false. is_mpi_root=.true. is_omp_root=.true. is_master=.true. ! Load tracer number and names: call infotrac_init ! allocate arrays allocate(q(iip1,jjp1,llm,nqtot)) allocate(qsurf(ngridmx,nqtot)) allocate(qsurf_tmp(ngridmx,nqtot)) ! get value of nsoilmx and allocate corresponding arrays ! a default value of nsoilmx is set in comsoil_h call getin_p("nsoilmx",nsoilmx) allocate(tsoil(ngridmx,nsoilmx)) allocate(ith(iip1,jjp1,nsoilmx),ithfi(ngridmx,nsoilmx)) c======================================================================= c Choice of the start file(s) to use c======================================================================= write(*,*) 'From which kind of files do you want to create new', . 'start and startfi files' write(*,*) ' 0 - from a file start_archive' write(*,*) ' 1 - from files start and startfi' c----------------------------------------------------------------------- c Open file(s) to modify (start or start_archive) c----------------------------------------------------------------------- DO read(*,*,iostat=ierr) choix_1 if ((choix_1 /= 0).OR.(choix_1 /=1)) EXIT ENDDO c Open start_archive c ~~~~~~~~~~~~~~~~~~~~~~~~~~ if (choix_1.eq.0) then write(*,*) 'Creating start files from:' write(*,*) './start_archive.nc' write(*,*) fichnom = 'start_archive.nc' ierr = NF_OPEN (fichnom, NF_NOWRITE,nid) IF (ierr.NE.NF_NOERR) THEN write(6,*)' Problem opening file:',fichnom write(6,*)' ierr = ', ierr CALL ABORT ENDIF tab0 = 50 Lmodif = 1 c OR open start and startfi files c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ else write(*,*) 'Creating start files from:' write(*,*) './start.nc and ./startfi.nc' write(*,*) fichnom = 'start.nc' ierr = NF_OPEN (fichnom, NF_NOWRITE,nid_dyn) IF (ierr.NE.NF_NOERR) THEN write(6,*)' Problem opening file:',fichnom write(6,*)' ierr = ', ierr CALL ABORT ENDIF fichnom = 'startfi.nc' ierr = NF_OPEN (fichnom, NF_NOWRITE,nid_fi) IF (ierr.NE.NF_NOERR) THEN write(6,*)' Problem opening file:',fichnom write(6,*)' ierr = ', ierr CALL ABORT ENDIF tab0 = 0 Lmodif = 0 endif c======================================================================= c INITIALISATIONS DIVERSES c======================================================================= ! Initialize global tracer indexes (stored in tracer.h) ! ... this has to be done before phyetat0 ! and requires that "datadir" be correctly initialized call getin_p("datadir",datadir) call initracer(ngridmx,nqtot) ! Initialize dimphy module (klon,klev,..) call init_dimphy(ngridmx,llm) ! Allocate saved arrays (as in firstcall of physiq) call phys_state_var_init(nqtot) c----------------------------------------------------------------------- c Lecture du tableau des parametres du run (pour la dynamique) c----------------------------------------------------------------------- if (choix_1.eq.0) then write(*,*) 'reading tab_cntrl START_ARCHIVE' c ierr = NF_INQ_VARID (nid, "controle", nvarid) #ifdef NC_DOUBLE ierr = NF_GET_VAR_DOUBLE(nid, nvarid, tab_cntrl) #else ierr = NF_GET_VAR_REAL(nid, nvarid, tab_cntrl) #endif c else if (choix_1.eq.1) then write(*,*) 'reading tab_cntrl START' c ierr = NF_INQ_VARID (nid_dyn, "controle", nvarid) #ifdef NC_DOUBLE ierr = NF_GET_VAR_DOUBLE(nid_dyn, nvarid, tab_cntrl) #else ierr = NF_GET_VAR_REAL(nid_dyn, nvarid, tab_cntrl) #endif c write(*,*) 'reading tab_cntrl STARTFI' c ierr = NF_INQ_VARID (nid_fi, "controle", nvarid) #ifdef NC_DOUBLE ierr = NF_GET_VAR_DOUBLE(nid_fi, nvarid, tab_cntrl_bis) #else ierr = NF_GET_VAR_REAL(nid_fi, nvarid, tab_cntrl_bis) #endif c do i=1,50 tab_cntrl(i+50)=tab_cntrl_bis(i) enddo write(*,*) 'printing tab_cntrl', tab_cntrl do i=1,100 write(*,*) i,tab_cntrl(i) enddo write(*,*) 'Reading file START' fichnom = 'start.nc' CALL dynetat0(fichnom,vcov,ucov,teta,q,masse, . ps,phis,time) CALL iniconst CALL inigeom ! Initialize the physics CALL iniphysiq(iim,jjm,llm, & (jjm-1)*iim+2,comm_lmdz, & daysec,day_ini,dtphys, & rlatu,rlatv,rlonu,rlonv, & aire,cu,cv,rad,g,r,cpp, & 1) ! Lmodif set to 0 to disable modifications possibility in phyeta0 Lmodif=0 write(*,*) 'Reading file STARTFI' fichnom = 'startfi.nc' CALL phyetat0(.true.,ngridmx,llm,fichnom,tab0,Lmodif,nsoilmx, . nqtot,day_ini,time, . tsurf,tsoil,emis,q2,qsurf) !) ! temporary modif by RDW ! copy albedo and soil thermal inertia on (local) physics grid do i=1,ngridmx albfi(i) = albedodat(i) do j=1,nsoilmx ithfi(i,j) = inertiedat(i,j) enddo ! build a surfithfi(:) using 1st layer of ithfi(:), which might ! be needed later on if reinitializing soil thermal inertia surfithfi(i)=ithfi(i,1) enddo ! also copy albedo and soil thermal inertia on (local) dynamics grid ! so that options below can manipulate either (but must then ensure ! to correctly recast things on physics grid) call gr_fi_dyn(1,ngridmx,iip1,jjp1,albfi,alb) call gr_fi_dyn(nsoilmx,ngridmx,iip1,jjp1,ithfi,ith) call gr_fi_dyn(1,ngridmx,iip1,jjp1,surfithfi,surfith) endif c----------------------------------------------------------------------- c Initialisation des constantes dynamique c----------------------------------------------------------------------- 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,ztot0,stot0,ang0" write(*,*) kappa,etot0,ptot0,ztot0,stot0,ang0 ! for vertical coordinate preff=tab_cntrl(18) ! reference surface pressure pa=tab_cntrl(17) ! reference pressure at which coord is purely pressure !NB: in start_archive files tab_cntrl(17)=tab_cntrl(18)=0 write(*,*) "Newstart: preff=",preff," pa=",pa yes=' ' do while ((yes.ne.'y').and.(yes.ne.'n')) write(*,*) "Change the values of preff and pa? (y/n)" read(*,fmt='(a)') yes end do if (yes.eq.'y') then write(*,*)"New value of reference surface pressure preff?" write(*,*)" (for Mars, typically preff=610)" read(*,*) preff write(*,*)"New value of reference pressure pa for purely" write(*,*)"pressure levels (for hybrid coordinates)?" write(*,*)" (for Mars, typically pa=20)" read(*,*) pa endif c----------------------------------------------------------------------- c Lecture du tab_cntrl et initialisation des constantes physiques c - pour start: Lmodif = 0 => pas de modifications possibles c (modif dans le tabfi de readfi + loin) c - pour start_archive: Lmodif = 1 => modifications possibles c----------------------------------------------------------------------- if (choix_1.eq.0) then ! tabfi requires that input file be first opened by open_startphy(fichnom) fichnom = 'start_archive.nc' call open_startphy(fichnom) call tabfi (ngridmx,nid,Lmodif,tab0,day_ini,lllm,p_rad, . p_omeg,p_g,p_cpp,p_mugaz,p_daysec,time) else if (choix_1.eq.1) then fichnom = 'startfi.nc' call open_startphy(fichnom) Lmodif=1 ! Lmodif set to 1 to allow modifications in phyeta0 call tabfi (ngridmx,nid_fi,Lmodif,tab0,day_ini,lllm,p_rad, . p_omeg,p_g,p_cpp,p_mugaz,p_daysec,time) endif if (p_omeg .eq. -9999.) then p_omeg = 8.*atan(1.)/p_daysec print*,"new value of omega",p_omeg endif c Initialisation coordonnees /aires c ------------------------------- ! Note: rlatu(:) and rlonv(:) are commons defined in "comgeom.h" ! rlatu() and rlonv() are given in radians latfi(1)=rlatu(1) lonfi(1)=0. DO j=2,jjm DO i=1,iim latfi((j-2)*iim+1+i)=rlatu(j) lonfi((j-2)*iim+1+i)=rlonv(i) ENDDO ENDDO latfi(ngridmx)=rlatu(jjp1) lonfi(ngridmx)=0. CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,aire,airefi) rad = p_rad omeg = p_omeg g = p_g cpp = p_cpp mugaz = p_mugaz daysec = p_daysec kappa = 8.314*1000./(p_mugaz*p_cpp) ! added by RDW c======================================================================= c INITIALISATIONS DIVERSES c======================================================================= ! Load parameters from run.def file CALL defrun_new( 99, .TRUE. ) ! Initialize dynamics geometry and co. (which, when using start.nc may ! have changed because of modifications (tabi, preff, pa) above) CALL iniconst CALL inigeom idum=-1 idum=0 ! Initialize the physics for start_archive only IF (choix_1.eq.0) THEN CALL iniphysiq(iim,jjm,llm, & (jjm-1)*iim+2,comm_lmdz, & daysec,day_ini,dtphys, & rlatu,rlatv,rlonu,rlonv, & aire,cu,cv,rad,g,r,cpp, & 1) ENDIF c======================================================================= c lecture topographie, albedo, inertie thermique, relief sous-maille c======================================================================= if (choix_1.eq.0) then ! for start_archive files, ! where an external "surface.nc" file is needed c do while((relief(1:3).ne.'mol').AND.(relief(1:3).ne.'pla')) c write(*,*) c write(*,*) 'choix du relief (mola,pla)' c write(*,*) '(Topographie MGS MOLA, plat)' c read(*,fmt='(a3)') relief relief="mola" c enddo ! First get the correct value of datadir, if not already done: ! default 'datadir' is set in "datafile_mod" call getin_p("datadir",datadir) write(*,*) 'Available surface data files are:' filestring='ls '//trim(datadir)//'/'// & trim(surfdir)//' | grep .nc' call system(filestring) ! but in ye old days these files were in datadir, so scan it as well ! for the sake of retro-compatibility filestring='ls '//trim(datadir)//'/'//' | grep .nc' call system(filestring) write(*,*) '' write(*,*) 'Please choose the relevant file', & ' (e.g. "surface_mars.nc")' write(*,*) ' or "none" to not use any (and set planetary' write(*,*) ' albedo and surface thermal inertia)' read(*,fmt='(a50)') surfacefile if (surfacefile.ne."none") then CALL datareadnc(relief,surfacefile,phis,alb,surfith, & zmeaS,zstdS,zsigS,zgamS,ztheS) else ! specific case when not using a "surface.nc" file phis(:,:)=0 zmeaS(:,:)=0 zstdS(:,:)=0 zsigS(:,:)=0 zgamS(:,:)=0 ztheS(:,:)=0 write(*,*) "Enter value of albedo of the bare ground:" read(*,*) alb(1,1) alb(:,:)=alb(1,1) write(*,*) "Enter value of thermal inertia of soil:" read(*,*) surfith(1,1) surfith(:,:)=surfith(1,1) endif ! of if (surfacefile.ne."none") CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,phis,phisfi) CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,surfith,surfithfi) CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,alb,albfi) CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,zmeaS,zmea) CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,zstdS,zstd) CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,zsigS,zsig) CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,zgamS,zgam) CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,ztheS,zthe) endif ! of if (choix_1.eq.0) c======================================================================= c Lecture des fichiers (start ou start_archive) c======================================================================= if (choix_1.eq.0) then write(*,*) 'Reading file START_ARCHIVE' CALL lect_start_archive(ngridmx,llm, & date,tsurf,tsoil,emis,q2, & t,ucov,vcov,ps,teta,phisold_newgrid,q,qsurf, & surfith,nid) write(*,*) "OK, read start_archive file" ! copy soil thermal inertia ithfi(:,:)=inertiedat(:,:) ierr= NF_CLOSE(nid) else if (choix_1.eq.1) then !do nothing, start and startfi have already been read else CALL exit(1) endif dtvr = daysec/FLOAT(day_step) dtphys = dtvr * FLOAT(iphysiq) c======================================================================= c c======================================================================= do ! infinite loop on list of changes write(*,*) write(*,*) write(*,*) 'List of possible changes :' write(*,*) '~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~' write(*,*) write(*,*) 'flat : no topography ("aquaplanet")' write(*,*) 'set_ps_to_preff : used if changing preff with topo' write(*,*) 'qname : change tracer name' write(*,*) 't=profile : read temperature profile in profile.in' write(*,*) 'q=0 : ALL tracer =zero' write(*,*) 'q=x : give a specific uniform value to one tracer' write(*,*) 'qs=x : give a uniform value to a surface tracer' write(*,*) 'q=profile : specify a profile for a tracer' write(*,*) 'qnogcm: initial tracer for GCM from nogcm (CH4,CO)' write(*,*) 'isotherm : Isothermal Temperatures' write(*,*) 'tempprof : specific Temperature profile' write(*,*) 'initsoil : initialize soil Temperatures' write(*,*) 'ptot : change total pressure' write(*,*) 'therm_ini_s: set soil TI to reference surf. values' write(*,*) 'inert3d: give uniform value of thermal inertia' write(*,*) 'subsoilice_n: put deep underground ice in N. hemis' write(*,*) 'subsoilice_s: put deep underground ice in S. hemis' write(*,*) 'reservoir: increase/decrease reservoir of ice' write(*,*) 'reservoir_SP: increase/decrease reservoir in SP' write(*,*) 'plutomap: initialize surface like pluto from ...' write(*,*) 'subsoil_n2: diu-sea thermal inertia for N2 only' write(*,*) 'subsoil_ch4: diu-sea thermal inertia for CH4 only' write(*,*) 'subsoil_all: diu-sea thermal inertia for all terr' write(*,*) 'subsoil_alb: diu-sea thermal inertia from albedo' write(*,*) 'diurnal_TI: diurnal thermal inertia for all terr' write(*,*) 'initsurf: initial surface temperature (global)' write(*,*) 'modtsurf: initial surface temperature (global)' write(*,*) 'copylat: copy tsurf and tsoil latitude' write(*,*) 'copylatlon: copy tsurf and tsoil lat/lon' write(*,*) 'copylon: copy tsurf tsoil lat, n2surf and phisfi' write(*,*) 'tsurfice: temperature depending on N2 ice' write(*,*) 'icarus: option to change surface/soil temperature' write(*,*) 'icarus_ch4: special option to add ch4' write(*,*) 'delfrostch4: delete ch4 frost' write(*,*) 'modch4: remove/modify amount ch4 frost' write(*,*) 'modch4n2: modify amount ch4 frost according N2' write(*,*) 'modco: remove/modify amount co frost' write(*,*) 'modn2: remove/modify amount n2 ice' write(*,*) 'modcoch4: remove/modify co ch4 where no n2 ' write(*,*) 'checktsoil: change tsoil where no n2 ' write(*,*) 'non2noco: if no n2, no co ' write(*,*) 'modn2_topo: modify n2 ice topo according to topo' write(*,*) 'modwheren2: modify n2 where already n2' write(*,*) 'modn2HD: modify n2 for HD runs' write(*,*) 'modn2HD_SP: modify n2 for HD runs in SP' write(*,*) 'globch4co: add/remove global amount ch4co frost' write(*,*) 'source_ch4: add source ch4' write(*,*) 'nomountain: remove pluto mountains (!)' write(*,*) 'relief: add pluto crater or mountain' write(*,*) 'seticeNH: set ice for initial start with NH topo' write(*,*) 'topo_sp: change topography of SP' write(*,*) 'fill_sp: fill sp with N2 ice and adjust phisfi' write(*,*) 'fill_sp_inv: fill inverted sp and adjust' write(*,*) 'subsoil_SP: diu-sea thermal inertia for SP ' write(*,*) 'bladed: add ch4 on bladed terrains' write(*,*) 'cpbladed: add extension bladed terrains' write(*,*) 'slope: add slope over all long. (for triton)' write(*,*) 'digsp: specific to dig SP' write(*,*) 'bugn2: correct bug of warm n2 due to HighRes' write(*,*) 'bugch4: correct bug of warm ch4 due to HighRes' write(*,*) 'flatnosp : topo flat except SP (topo controled)' write(*,*) 'flatregion: topo flat for specific region' write(*,*) 'changetopo: change topo' write(*,*) 'asymtopo: N-S asym topo in tanh' write(*,*) 'corrtopo: correction topo bug' write(*,*) 'adjustphi: adjust phisfi according to N2 mass' write(*,*) 'correctphi: adjust phisfi' write(*,*) 'correctps: test to correct ps' write(*,*) 'toponoise: no flat topo' write(*,*) 'asymtriton: asymetry in longitude for triton' write(*,*) 'copytsoil: copy 2D tsoil field from nc file' write(*,*) 'albedomap: read in an albedomap albedo.nc' write(*,*) 'mod_distrib_ice: modify ice distrib. from albedo' write(*,*) write(*,*) 'Please note that emis and albedo are set in physiq' write(*,*) write(*,*) write(*,*) 'Change to perform ?' write(*,*) ' (enter keyword or return to end)' write(*,*) read(*,fmt='(a20)') modif if (modif(1:1) .eq. ' ')then exit ! exit loop on changes endif write(*,*) write(*,*) trim(modif) , ' : ' c 'flat : no topography ("aquaplanet")' c ------------------------------------- if (trim(modif) .eq. 'flat') then c set topo to zero z_reel(1:iip1,1:jjp1)=0 phis(1:iip1,1:jjp1)=g*z_reel(1:iip1,1:jjp1) CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,phis,phisfi) write(*,*) 'topography set to zero.' write(*,*) 'WARNING : the subgrid topography parameters', & ' were not set to zero ! => set calllott to F' c Choice of surface pressure yes=' ' do while ((yes.ne.'y').and.(yes.ne.'n')) write(*,*) 'Do you wish to choose homogeneous surface', & 'pressure (y) or let newstart interpolate ', & ' the previous field (n)?' read(*,fmt='(a)') yes end do if (yes.eq.'y') then flagps0=.true. write(*,*) 'New value for ps (Pa) ?' 201 read(*,*,iostat=ierr) patm if(ierr.ne.0) goto 201 write(*,*) patm if (patm.eq.-9999.) then patm = preff write(*,*) "we set patm = preff", preff endif write(*,*) write(*,*) ' new ps everywhere (Pa) = ', patm write(*,*) do j=1,jjp1 do i=1,iip1 ps(i,j)=patm enddo enddo end if c 'set_ps_to_preff' : used if changing preff with topo c ---------------------------------------------------- else if (trim(modif) .eq. 'set_ps_to_preff') then do j=1,jjp1 do i=1,iip1 ps(i,j)=preff enddo enddo c qname : change tracer name c -------------------------- 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,nqtot write(*,'(i3,a3,a20)')iq,' : ',trim(tname(iq)) enddo write(*,'(a35,i3)') & '(enter tracer number; between 1 and ',nqtot write(*,*)' or any other value to quit this option)' read(*,*) iq if ((iq.ge.1).and.(iq.le.nqtot)) 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') c q=0 : set tracers to zero c ------------------------- else if (trim(modif).eq.'q=0') then c mise a 0 des q (traceurs) write(*,*) 'Tracers set to 0 (1.E-30 in fact)' DO iq =1, nqtot DO l=1,llm DO j=1,jjp1 DO i=1,iip1 q(i,j,l,iq)=1.e-30 ENDDO ENDDO ENDDO ENDDO c set surface tracers to zero DO iq =1, nqtot DO ig=1,ngridmx qsurf(ig,iq)=0. ENDDO ENDDO c q=x : initialise tracer manually c -------------------------------- else if (trim(modif).eq.'q=x') then write(*,*) 'Which tracer do you want to modify ?' do iq=1,nqtot write(*,*)iq,' : ',trim(tname(iq)) enddo write(*,*) '(choose between 1 and ',nqtot,')' read(*,*) iq write(*,*)'mixing ratio of tracer ',trim(tname(iq)), & ' ? (kg/kg)' read(*,*) val DO l=1,llm DO j=1,jjp1 DO i=1,iip1 q(i,j,l,iq)=val ENDDO ENDDO ENDDO write(*,*) 'SURFACE value of tracer ',trim(tname(iq)), & ' ? (kg/m2)' read(*,*) val DO ig=1,ngridmx qsurf(ig,iq)=val ENDDO c qs=x : initialise surface tracer manually c -------------------------------- else if (trim(modif).eq.'qs=x') then write(*,*) 'Which tracer do you want to modify ?' do iq=1,nqtot write(*,*)iq,' : ',trim(tname(iq)) enddo write(*,*) '(choose between 1 and ',nqtot,')' read(*,*) iq write(*,*) 'SURFACE value of tracer ',trim(tname(iq)), & ' ? (kg/m2)' read(*,*) val DO ig=1,ngridmx qsurf(ig,iq)=val ENDDO c t=profile : initialize temperature with a given profile c ------------------------------------------------------- else if (trim(modif) .eq. 't=profile') then write(*,*) 'Temperature profile from ASCII file' write(*,*) "'profile.in' e.g. 1D output" 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)" txt="profile.in" open(33,file=trim(txt),status='old',form='formatted', & iostat=ierr) if (ierr.eq.0) then ! OK, found file 'profile_...', load the profile do l=1,llm+1 read(33,*,iostat=ierr) profile(l) write(*,*) profile(l) if (ierr.ne.0) then ! something went wrong exit ! quit loop endif enddo if (ierr.eq.0) then tsurf(1:ngridmx)=profile(1) tsoil(1:ngridmx,1:nsoilmx)=profile(1) do l=1,llm Tset(1:iip1,1:jjp1,l)=profile(l+1) flagtset=.true. enddo ucov(1:iip1,1:jjp1,1:llm)=0. vcov(1:iip1,1:jjm,1:llm)=0. q2(1:ngridmx,1:llm+1)=0. else write(*,*)'problem reading file ',trim(txt),' !' write(*,*)'No modifications to temperature' endif else write(*,*)'Could not find file ',trim(txt),' !' write(*,*)'No modifications to temperature' endif c q=profile : initialize tracer with a given profile c -------------------------------------------------- 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)" write(*,*) 'Which tracer do you want to set?' do iq=1,nqtot write(*,*)iq,' : ',trim(tname(iq)) enddo write(*,*) '(choose between 1 and ',nqtot,')' read(*,*) iq if ((iq.lt.1).or.(iq.gt.nqtot)) then ! wrong value for iq, go back to menu write(*,*) "wrong input value:",iq cycle endif ! 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,llm+1 read(41,*,iostat=ierr) profile(l) if (ierr.ne.0) then ! something went wrong exit ! quit loop endif enddo if (ierr.eq.0) then ! initialize tracer values qsurf(:,iq)=profile(1) do l=1,llm q(:,:,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 c qnogcm : initialise tracer from nogcm (CH4, CO) c -------------------------------- else if (modif(1:len_trim(modif)).eq.'qnogcm') then write(*,*) 'Which tracer do you want to modify ?' write(*,*) 'Should be ch4_gas and co_gas' do iq=1,nqtot write(*,*)iq,' : ',trim(noms(iq)),' : ',q(1,1,1,iq) enddo write(*,*) '(choose between 1 and ',nqtot,')' read(*,*) iq DO l=1,llm DO j=1,jjp1 DO i=1,iip1 q(i,j,l,iq)=q(i,j,1,iq) ENDDO ENDDO ENDDO c isothermal temperatures c ------------------------------------------------ else if (modif(1:len_trim(modif)) .eq. 'isotherm') then write(*,*)'Isothermal temperature of the atmosphere' write(*,*) 'Value of THIS temperature ? :' 203 read(*,*,iostat=ierr) Tset(1,1,1) if(ierr.ne.0) goto 203 flagtset=.true. DO l=1,llm DO j=1,jjp1 DO i=1,iip1 Tset(i,j,l)=Tset(1,1,1) ENDDO ENDDO ENDDO write(*,*) 'atmospheric temp =', Tset(2,2,2) c specific temperature profile c ------------------------------------------------ else if (modif(1:len_trim(modif)) .eq. 'tempprof') then write(*,*) 'phi=' write(*,*) phi(1,1,:) write(*,*) 'temperature profile in the atmosphere' write(*,*) 'Value of THIS temperature ? :' 204 read(*,*,iostat=ierr) Tset(1,1,1) if(ierr.ne.0) goto 204 flagtset=.true. DO l=1,llm DO j=1,jjp1 DO i=1,iip1 Tset(i,j,l)=Tset(1,1,1) ENDDO ENDDO ENDDO write(*,*) 'atmospheric temp =', Tset(2,2,2) c initsoil: subsurface temperature c --------------------------------- else if (modif(1:len_trim(modif)) .eq. 'initsoil') then write(*,*)'Isothermal temperature of the subsurface' write(*,*) 'Value of this temperature ? :' 303 read(*,*,iostat=ierr) Tiso if(ierr.ne.0) goto 303 do l=1,nsoilmx do ig=1, ngridmx tsoil(ig,l) = Tiso end do end do c icarus : changing tsoil tsurf on latitudinal bands c ----------------------------------------------------- else if (modif(1:len_trim(modif)) .eq. 'icarus') then write(*,*) 'At which lat lon do you want to extract the & reference tsurf/tsoil profile ?' 407 read(*,*,iostat=ierr) val,val2 if(ierr.ne.0) goto 407 write(*,*) 'You want lat =',val write(*,*) 'You want lon =',val2 dist=0 ref=1000 do j=1,ngridmx-1 dist=sqrt((lonfi(j)*180./pi-val2)**2+ & (latfi(j)*180./pi-val)**2) IF (dist.lt.ref) then ref=dist jref=j ENDIF ENDDO write(*,*)'Will use nearest grid point which is:', & latfi(jref)*180./pi,lonfi(jref)*180./pi ! Extraction of the profile write(*,*) 'Profile is =',tsoil(jref,:) write(*,*) 'tsurf is =',tsurf(jref) write(*,*) 'Choice lat for latitudinal band with this profile' 408 read(*,*,iostat=ierr) val3 if(ierr.ne.0) goto 408 write(*,*) 'Choice delta lat for this latitudinal band' 409 read(*,*,iostat=ierr) val4 if(ierr.ne.0) goto 409 write(*,*) 'Choice delta temp (K) for this latitudinal band' 410 read(*,*,iostat=ierr) val5 if(ierr.ne.0) goto 410 write(*,*) 'How much N2 ice should I put on this band (kg/m2)' 415 read(*,*,iostat=ierr) choice if(ierr.ne.0) goto 415 write(*,*) 'Choice lat2' 411 read(*,*,iostat=ierr) val6 if(ierr.ne.0) goto 411 write(*,*) 'Choice delta lat for this latitudinal band' 412 read(*,*,iostat=ierr) val7 if(ierr.ne.0) goto 412 write(*,*) 'Choice delta temp (K) for this latitudinal band' 413 read(*,*,iostat=ierr) val8 if(ierr.ne.0) goto 413 DO ig=1,ngridmx IF ( (latfi(ig)*180./pi.ge.(val3-val4/2.)) .and. & (latfi(ig)*180./pi.le.(val3+val4/2.)) .and. & (qsurf(ig,igcm_n2).lt.0.001) ) then tsurf(ig)=tsurf(jref)+val5 DO l=1,nsoilmx tsoil(ig,l)=tsoil(jref,l)+val5 ENDDO qsurf(ig,igcm_n2)=choice ENDIF IF ( (latfi(ig)*180./pi.ge.(val6-val7/2.)) .and. & (latfi(ig)*180./pi.le.(val6+val7/2.)) .and. & (qsurf(ig,igcm_n2).lt.0.001) ) then tsurf(ig)=tsurf(jref)+val8 DO l=1,nsoilmx tsoil(ig,l)=tsoil(jref,l)+val8 ENDDO ENDIF ENDDO c icarus_ch4 : changing tsoil tsurf and adding ch4 c ----------------------------------------------------- else if (modif(1:len_trim(modif)) .eq. 'icarus_ch4') then write(*,*) 'At which lat lon do you want to extract the & reference tsurf/tsoil profile ?' 416 read(*,*,iostat=ierr) val,val2 if(ierr.ne.0) goto 416 write(*,*) 'You want lat =',val write(*,*) 'You want lon =',val2 dist=0 ref=1000 do j=1,ngridmx-1 dist=sqrt((lonfi(j)*180./pi-val2)**2+ & (latfi(j)*180./pi-val)**2) IF (dist.lt.ref) then ref=dist jref=j ENDIF ENDDO write(*,*)'Will use nearest grid point which is:', & latfi(jref)*180./pi,lonfi(jref)*180./pi ! Extraction of the profile write(*,*) 'Profile is =',tsoil(jref,:) write(*,*) 'tsurf is =',tsurf(jref) write(*,*) 'Choice band : lat1 and lat2 ?' 417 read(*,*,iostat=ierr) val3,val4 if(ierr.ne.0) goto 417 write(*,*) 'Choice temp coefficient for this latitudinal band' 418 read(*,*,iostat=ierr) val5 if(ierr.ne.0) goto 418 write(*,*) 'How much CH4 ice do I put on this band (kg/m2)' 419 read(*,*,iostat=ierr) choice if(ierr.ne.0) goto 419 DO ig=1,ngridmx IF ( (latfi(ig)*180./pi.ge.val3) .and. & (latfi(ig)*180./pi.le.val4) .and. & (qsurf(ig,igcm_ch4_ice).lt.0.001) ) then tsurf(ig)=tsurf(jref)*val5 DO l=1,nsoilmx tsoil(ig,l)=tsoil(jref,l)*val5 ENDDO qsurf(ig,igcm_ch4_ice)=choice ENDIF ENDDO c globch4co : adding/remove global ch4 co c ---------------------------------------------------------- else if (modif(1:len_trim(modif)) .eq. 'globch4co') then write(*,*) 'Adding or removing ch4 co ' write(*,*) 'Choice amount add/less ch4 ice (0 = rm all)' 431 read(*,*,iostat=ierr) val3 if(ierr.ne.0) goto 431 write(*,*) 'Choice amount add/less co ice (0 = rm all)' 432 read(*,*,iostat=ierr) val4 if(ierr.ne.0) goto 432 IF (val3.eq.0.) then DO ig=1,ngridmx qsurf(ig,igcm_ch4_ice)=0. ENDDO ELSE DO ig=1,ngridmx qsurf(ig,igcm_ch4_ice)=qsurf(ig,igcm_ch4_ice)+val3 ENDDO ENDIF IF (val4.eq.0.) then DO ig=1,ngridmx qsurf(ig,igcm_co_ice)=0. ENDDO ELSE DO ig=1,ngridmx qsurf(ig,igcm_co_ice)=qsurf(ig,igcm_co_ice)+val4 ENDDO ENDIF c bladed : adding/remove ch4 on bladed terrains c ---------------------------------------------------------- else if (modif(1:len_trim(modif)).eq.'bladed') then write(*,*) 'Adding or removing ch4 on the bladed terrains' write(*,*) 'Choice band : lat1 and lat2 ?' 450 read(*,*,iostat=ierr) val,val2 if(ierr.ne.0) goto 450 write(*,*) 'Choice band : lon1 and lon2 ?' 451 read(*,*,iostat=ierr) val3,val4 if(ierr.ne.0) goto 451 write(*,*) 'Choice phisfi minimum ?' 454 read(*,*,iostat=ierr) choice if(ierr.ne.0) goto 454 write(*,*) 'Choice multiplicative factor' 452 read(*,*,iostat=ierr) val5 if(ierr.ne.0) goto 452 write(*,*) 'Choice additional ch4' 453 read(*,*,iostat=ierr) val6 if(ierr.ne.0) goto 453 write(*,*) 'Choice index for tsurf tsoil' 449 read(*,*,iostat=ierr) iref if(ierr.ne.0) goto 449 ! shift DO ig=1,ngridmx IF ( (latfi(ig)*180./pi.ge.val) .and. & (latfi(ig)*180./pi.le.val2) .and. & (lonfi(ig)*180./pi.ge.val3) .and. & (lonfi(ig)*180./pi.le.val4) .and. & (phisfi(ig).gt.choice) ) then qsurf(ig,igcm_ch4_ice)=qsurf(ig,igcm_ch4_ice)*val5 qsurf(ig,igcm_ch4_ice)=qsurf(ig,igcm_ch4_ice)+val6 tsurf(ig)=tsurf(iref) DO l=1,nsoilmx tsoil(ig,l)=tsoil(iref,l) ENDDO ENDIF ENDDO c cpbladed : add extension bladed terrains c ---------------------------------------------------------- else if (modif(1:len_trim(modif)).eq.'cpbladed') then write(*,*) 'Choice lat,lon, centre of band to be copied ?' 560 read(*,*,iostat=ierr) val,val2 if(ierr.ne.0) goto 560 write(*,*) 'Choice distance (km) from there defining the area' 561 read(*,*,iostat=ierr) val3 if(ierr.ne.0) goto 561 write(*,*) 'Nb of copies ?' 562 read(*,*,iostat=ierr) l if(ierr.ne.0) goto 562 ! Get index correponding to central points ! 1/ Reference point igref=1 actualmin=1000. do ig=1,ngridmx val6=dist_pluto(latfi(ig),lonfi(ig), & val*pi/180.,val2*pi/180.) if (val6.lt.actualmin) then actualmin=val6 igref=ig endif enddo do k=1,l write(*,*) 'Choice lat,lon where terrain is copied' 563 read(*,*,iostat=ierr) val4,val5 if(ierr.ne.0) goto 563 if (val5.gt.180.) then val5=val5-360. endif ! 2/ Target point igref2=1 actualmin=1000. do ig=1,ngridmx val6=dist_pluto(latfi(ig),lonfi(ig), & val4*pi/180.,val5*pi/180.) if (val6.lt.actualmin) then actualmin=val6 igref2=ig endif enddo ! for each point within A1, get distance and angle ! save in a table i=1 array_ind(:)=0 array_dist(:)=1000. array_angle(:)=0. do ig=1,ngridmx val6=dist_pluto(latfi(ig),lonfi(ig), & val*pi/180.,val2*pi/180.) if (val6.lt.val3) then array_ind(i)=ig array_dist(i)=val6 array_angle(i)=atan2(val-latfi(ig)*180./pi, & val2-lonfi(ig)*180./pi) i=i+1 endif enddo ! for each point within A2, get distance and angle ! and look where fit with previous table, and set value do ig=1,ngridmx val6=dist_pluto(latfi(ig),lonfi(ig), & val4*pi/180.,val5*pi/180.) if (val6.lt.val3) then ! dist = val6 ! get angle: angle=atan2(val4-latfi(ig)*180./pi, & val5-lonfi(ig)*180./pi) ! Loop where min actualmin=1000. do j=1,i if ( (array_angle(j).lt.angle+0.52).and. & (array_angle(j).gt.angle-0.52).and. & (array_dist(j)-val6.lt.actualmin) ) then actualmin=array_dist(j)-val6 igref3=j endif enddo phisfi(ig)=phisfi(array_ind(igref3)) qsurf(ig,igcm_ch4_ice)= & qsurf(array_ind(igref3),igcm_ch4_ice) tsurf(ig)=tsurf(array_ind(igref3)) endif enddo enddo CALL gr_fi_dyn(1,ngridmx,iip1,jjp1,phisfi,phis) c delfrostch4/ delete ch4 frost on a latitudinal band c ---------------------------------------------------------- else if (modif(1:len_trim(modif)) .eq. 'delfrostch4') then write(*,*) 'removing ch4 latitudinally' write(*,*) 'Choice band : lat1 and lat2 ?' read(*,*,iostat=ierr) val,val2 write(*,*) 'Choice band : lon1 and lon2 ?' 522 read(*,*,iostat=ierr) val4,val5 if(ierr.ne.0) goto 522 write(*,*) 'Choice amount max below whcih ch4 is removed' 521 read(*,*,iostat=ierr) val3 if(ierr.ne.0) goto 521 DO ig=1,ngridmx IF ( (latfi(ig)*180./pi.ge.val) .and. & (latfi(ig)*180./pi.le.val2) .and. & (lonfi(ig)*180./pi.ge.val4) .and. & (lonfi(ig)*180./pi.lt.val5) .and. & (qsurf(ig,igcm_ch4_ice).lt.val3) ) then qsurf(ig,igcm_ch4_ice)=0. ENDIF ENDDO c modch4 : adding/remove ch4 frost on a latitudinal band c ---------------------------------------------------------- else if (modif(1:len_trim(modif)) .eq. 'modch4') then write(*,*) 'Adding or removing ch4 latitudinally' write(*,*) 'Choice band : lat1 and lat2 ?' read(*,*,iostat=ierr) val,val2 write(*,*) 'Choice band : lon1 and lon2 ?' 422 read(*,*,iostat=ierr) val4,val5 if(ierr.ne.0) goto 422 write(*,*) 'Choice multiplicative factor' 421 read(*,*,iostat=ierr) val3 if(ierr.ne.0) goto 421 write(*,*) 'Choice additional ch4' 423 read(*,*,iostat=ierr) val6 if(ierr.ne.0) goto 423 DO ig=1,ngridmx IF ( (latfi(ig)*180./pi.ge.val) .and. & (latfi(ig)*180./pi.le.val2) .and. & (lonfi(ig)*180./pi.ge.val4) .and. & (lonfi(ig)*180./pi.lt.val5)) then ! & (qsurf(ig,igcm_n2).gt.50)) then ! & (qsurf(ig,igcm_ch4_ice).lt.10) ) then qsurf(ig,igcm_ch4_ice)=qsurf(ig,igcm_ch4_ice)*val3 qsurf(ig,igcm_ch4_ice)=qsurf(ig,igcm_ch4_ice)+val6 ENDIF ENDDO c modch4n2 : adding/remove ch4 frost on a latitudinal band c ---------------------------------------------------------- else if (modif(1:len_trim(modif)) .eq. 'modch4n2') then write(*,*) 'Adding or removing ch4 latitudinally' write(*,*) 'Choice band : lat1 and lat2 ?' read(*,*,iostat=ierr) val,val2 write(*,*) 'Choice band : lon1 and lon2 ?' read(*,*,iostat=ierr) val4,val5 write(*,*) 'Choice range n2 for modif' read(*,*,iostat=ierr) val7,val8 write(*,*) 'Choice multiplicative factor ch4' read(*,*,iostat=ierr) val3 write(*,*) 'Choice additional ch4' read(*,*,iostat=ierr) val6 DO ig=1,ngridmx IF ( (latfi(ig)*180./pi.ge.val) .and. & (latfi(ig)*180./pi.le.val2) .and. & (lonfi(ig)*180./pi.ge.val4) .and. & (lonfi(ig)*180./pi.lt.val5) .and. & (qsurf(ig,igcm_n2).gt.val7) .and. & (qsurf(ig,igcm_n2).lt.val8) ) then qsurf(ig,igcm_ch4_ice)=qsurf(ig,igcm_ch4_ice)*val3 qsurf(ig,igcm_ch4_ice)=qsurf(ig,igcm_ch4_ice)+val6 ENDIF ENDDO c non2noco : if no n2 no co c ---------------------------------------------------------- else if (modif(1:len_trim(modif)) .eq. 'non2noco') then DO ig=1,ngridmx IF ( (qsurf(ig,igcm_n2).lt.0.001)) then qsurf(ig,igcm_co_ice)=0. ENDIF ENDDO c modco : adding/remove co frost on a latitudinal band c ---------------------------------------------------------- else if (modif(1:len_trim(modif)) .eq. 'modco') then write(*,*) 'Adding or removing co where N2 is present ' write(*,*) 'Choice limit mini n2 pour mettre co?' 460 read(*,*,iostat=ierr) val7 if(ierr.ne.0) goto 460 write(*,*) 'Choice band : lat1 and lat2 ?' 461 read(*,*,iostat=ierr) val,val2 if(ierr.ne.0) goto 461 write(*,*) 'Choice band : lon1 and lon2 ?' 462 read(*,*,iostat=ierr) val4,val5 if(ierr.ne.0) goto 462 write(*,*) 'Choice multiplicative factor' 463 read(*,*,iostat=ierr) val3 if(ierr.ne.0) goto 463 write(*,*) 'Choice additional co' 464 read(*,*,iostat=ierr) val6 if(ierr.ne.0) goto 464 DO ig=1,ngridmx IF ( (latfi(ig)*180./pi.ge.val) .and. & (latfi(ig)*180./pi.le.val2) .and. & (lonfi(ig)*180./pi.ge.val4) .and. & (lonfi(ig)*180./pi.lt.val5) .and. & (qsurf(ig,igcm_n2).lt.val7)) then qsurf(ig,igcm_co_ice)=qsurf(ig,igcm_co_ice)*val3 qsurf(ig,igcm_co_ice)=qsurf(ig,igcm_co_ice)+val6 ENDIF ENDDO c modn2 : adding/remove n2 ice c ---------------------------------------------------------- else if (modif(1:len_trim(modif)) .eq. 'modn2') then write(*,*) 'Adding or removing n2 ' write(*,*) 'Choice band : lat1 and lat2 ?' 425 read(*,*,iostat=ierr) val,val2 if(ierr.ne.0) goto 425 write(*,*) 'Choice band : lon1 and lon2 ?' 426 read(*,*,iostat=ierr) val4,val5 if(ierr.ne.0) goto 426 write(*,*) 'Choice multiplicative factor' 427 read(*,*,iostat=ierr) val3 if(ierr.ne.0) goto 427 write(*,*) 'Choice additional n2' 428 read(*,*,iostat=ierr) val6 if(ierr.ne.0) goto 428 DO ig=1,ngridmx IF ( (latfi(ig)*180./pi.ge.val) .and. & (latfi(ig)*180./pi.le.val2) .and. & (lonfi(ig)*180./pi.ge.val4) .and. & (lonfi(ig)*180./pi.lt.val5) ) then c & (qsurf(ig,igcm_n2).lt.50)) then qsurf(ig,igcm_n2)=qsurf(ig,igcm_n2)*val3 qsurf(ig,igcm_n2)=qsurf(ig,igcm_n2)+val6 ENDIF ! IF ((phisfi(ig).gt.-1000.)) then ! qsurf(ig,igcm_n2)=0. ! ENDIF ENDDO c modcoch4 : adding/remove co and ch4 frost where co without n2 c ---------------------------------------------------------- else if (modif(1:len_trim(modif)) .eq. 'modcoch4') then write(*,*) 'Adding or removing co where N2 is not present ' write(*,*) 'Choice band : lat1 and lat2 ?' 491 read(*,*,iostat=ierr) val,val2 if(ierr.ne.0) goto 491 write(*,*) 'Choice band : lon1 and lon2 ?' 492 read(*,*,iostat=ierr) val4,val5 if(ierr.ne.0) goto 492 write(*,*) 'Choice multiplicative factor' 493 read(*,*,iostat=ierr) val3 if(ierr.ne.0) goto 493 write(*,*) 'Choice additional co ch4' 494 read(*,*,iostat=ierr) val6 if(ierr.ne.0) goto 494 DO ig=1,ngridmx IF ( (latfi(ig)*180./pi.ge.val) .and. & (latfi(ig)*180./pi.le.val2) .and. & (lonfi(ig)*180./pi.ge.val4) .and. & (lonfi(ig)*180./pi.le.val5) .and. & (qsurf(ig,igcm_n2).lt.10.)) then qsurf(ig,igcm_co_ice)=qsurf(ig,igcm_co_ice)*val3 qsurf(ig,igcm_co_ice)=qsurf(ig,igcm_co_ice)+val6 qsurf(ig,igcm_ch4_ice)=qsurf(ig,igcm_ch4_ice)*val3 qsurf(ig,igcm_ch4_ice)=qsurf(ig,igcm_ch4_ice)+val6 ENDIF ENDDO c modn2HD : adding/remove n2 ice for HD runs c ---------------------------------------------------------- else if (modif(1:len_trim(modif)) .eq. 'modn2HD') then write(*,*) 'Adding or removing n2 ' write(*,*) 'Choice band : lat1 and lat2 ?' 480 read(*,*,iostat=ierr) val,val1 if(ierr.ne.0) goto 480 write(*,*) 'Choice band : lon1 and lon2 ?' 481 read(*,*,iostat=ierr) val2,val3 if(ierr.ne.0) goto 481 write(*,*) 'Choice amount N2 min max where to modify' 482 read(*,*,iostat=ierr) val4,val5 if(ierr.ne.0) goto 482 write(*,*) 'Choice phisfi min max where change n2' 483 read(*,*,iostat=ierr) val6,val11 if(ierr.ne.0) goto 483 write(*,*) 'Choice multiplicative factor' 484 read(*,*,iostat=ierr) val7 if(ierr.ne.0) goto 484 write(*,*) 'Choice additional n2' 485 read(*,*,iostat=ierr) val8 if(ierr.ne.0) goto 485 ! write(*,*) 'Choice ind lon equivalent for change tsurf tsoil' ! 486 read(*,*,iostat=ierr) val9 ! if(ierr.ne.0) goto 486 DO ig=1,ngridmx IF ( (latfi(ig)*180./pi.ge.val) .and. & (latfi(ig)*180./pi.le.val1) .and. & (lonfi(ig)*180./pi.ge.val2) .and. & (lonfi(ig)*180./pi.lt.val3) .and. & (qsurf(ig,igcm_n2).ge.val4) .and. & (qsurf(ig,igcm_n2).lt.val5) .and. & (phisfi(ig).ge.val6) .and. & (phisfi(ig).le.val11) ) then qsurf(ig,igcm_n2)=qsurf(ig,igcm_n2)*val7 qsurf(ig,igcm_n2)=qsurf(ig,igcm_n2)+val8 !qsurf(ig,igcm_ch4_ice)=0. qsurf(ig,igcm_co_ice)=0. ! index of ig !if (val9.eq.0.) then ! iref=2546 !else ! val10=modulo((ig-1),96) ! choice=ig+int(96/2)-val10 ! !choice=int(1+96*(val9-1)+val10) !endif !tsurf(ig)=tsurf(iref) !DO l=1,nsoilmx ! tsoil(ig,l)=tsoil(iref,l) !ENDDO tsurf(ig)=tsurf(ig)+4 DO l=1,nsoilmx tsoil(ig,l)=tsoil(ig,l)+4 ENDDO ENDIF ! IF ((phisfi(ig).gt.-1000.)) then ! qsurf(ig,igcm_n2)=0. ! ENDIF ENDDO c modn2HD_SP : adding/remove n2 ice for HD runs c ---------------------------------------------------------- else if (modif(1:len_trim(modif)) .eq. 'modn2HD_SP') then write(*,*) 'Adding or removing n2 ' write(*,*) 'Choice band : lat1 and lat2 ?' 501 read(*,*,iostat=ierr) val,val1 if(ierr.ne.0) goto 501 write(*,*) 'Choice band : lon1 and lon2 ?' 502 read(*,*,iostat=ierr) val2,val3 if(ierr.ne.0) goto 502 write(*,*) 'Choice amount N2 min max where to modify' 503 read(*,*,iostat=ierr) val4,val5 if(ierr.ne.0) goto 503 write(*,*) 'Choice phisfi min max where change n2' 504 read(*,*,iostat=ierr) val6,val11 if(ierr.ne.0) goto 504 write(*,*) 'Choice multiplicative factor' 505 read(*,*,iostat=ierr) val7 if(ierr.ne.0) goto 505 write(*,*) 'Choice additional n2' 506 read(*,*,iostat=ierr) val8 if(ierr.ne.0) goto 506 write(*,*) 'Choice ind lon equivalent for change tsurf tsoil' 507 read(*,*,iostat=ierr) iref if(ierr.ne.0) goto 507 DO ig=1,ngridmx IF ( (latfi(ig)*180./pi.ge.val) .and. & (latfi(ig)*180./pi.le.val1) .and. & (lonfi(ig)*180./pi.ge.val2) .and. & (lonfi(ig)*180./pi.lt.val3) .and. & (qsurf(ig,igcm_n2).ge.val4) .and. & (qsurf(ig,igcm_n2).lt.val5) .and. & (phisfi(ig).ge.val6) .and. & (phisfi(ig).le.val11) ) then qsurf(ig,igcm_n2)=qsurf(ig,igcm_n2)*val7 qsurf(ig,igcm_n2)=qsurf(ig,igcm_n2)+val8 qsurf(ig,igcm_ch4_ice)=0. qsurf(ig,igcm_co_ice)=0. ! index of ig !if (val9.eq.0.) then ! iref=2546 !else !val10=modulo((ig-1),96) !choice=ig+int(96/2)-val10 !choice=int(1+96*(val9-1)+val10) if (iref.ge.1) then tsurf(ig)=tsurf(iref) DO l=1,nsoilmx tsoil(ig,l)=tsoil(iref,l) ENDDO else if (iref.eq.0) then choice=int(ig/96.)*96+84 print*, 'choice=',choice tsurf(ig)=tsurf(int(choice)) DO l=1,nsoilmx tsoil(ig,l)=tsoil(int(choice),l) ENDDO endif ENDIF ENDDO c modn2_topo : adding/remove n2 ice c ---------------------------------------------------------- else if (modif(1:len_trim(modif)) .eq. 'modn2_topo') then write(*,*) 'Adding or removing n2 ' write(*,*) 'Choice band : lat1 and lat2 ?' 441 read(*,*,iostat=ierr) val,val2 if(ierr.ne.0) goto 441 write(*,*) 'Choice band : lon1 and lon2 ?' 442 read(*,*,iostat=ierr) val4,val5 if(ierr.ne.0) goto 442 write(*,*) 'Choice topo 1 et 2 (phi) where change is active' 443 read(*,*,iostat=ierr) val7,val8 if(ierr.ne.0) goto 443 write(*,*) 'Choice multiplicative factor' 444 read(*,*,iostat=ierr) val3 if(ierr.ne.0) goto 444 write(*,*) 'Choice additional n2' 445 read(*,*,iostat=ierr) val6 if(ierr.ne.0) goto 445 DO ig=1,ngridmx IF ( (latfi(ig)*180./pi.ge.val) .and. & (latfi(ig)*180./pi.le.val2) .and. & (lonfi(ig)*180./pi.ge.val4) .and. & (lonfi(ig)*180./pi.lt.val5) .and. & (phisfi(ig).le.val8) .and. & (phisfi(ig).ge.val7) ) then qsurf(ig,igcm_n2)=qsurf(ig,igcm_n2)*val3 qsurf(ig,igcm_n2)=qsurf(ig,igcm_n2)+val6 ENDIF ENDDO c modwheren2 : adding/remove n2 ice where already n2 c ---------------------------------------------------------- else if (modif(1:len_trim(modif)) .eq. 'modwheren2') then write(*,*) 'Choice multiplicative factor' 471 read(*,*,iostat=ierr) val3 if(ierr.ne.0) goto 471 write(*,*) 'Choice additional n2' 472 read(*,*,iostat=ierr) val6 if(ierr.ne.0) goto 472 DO ig=1,ngridmx IF ( qsurf(ig,igcm_n2).gt.0.001 ) then qsurf(ig,igcm_n2)=qsurf(ig,igcm_n2)*val3 qsurf(ig,igcm_n2)=qsurf(ig,igcm_n2)+val6 ENDIF ENDDO c checktsoil : changing tsoil if no n2 c ---------------------------------------------------------- else if (modif(1:len_trim(modif)) .eq. 'checktsoil') then write(*,*) 'selecting area where tsoil to be check ' write(*,*) 'Choice band : lat1 and lat2 ?' 511 read(*,*,iostat=ierr) val,val1 if(ierr.ne.0) goto 511 write(*,*) 'Choice band : lon1 and lon2 ?' 512 read(*,*,iostat=ierr) val2,val3 if(ierr.ne.0) goto 512 write(*,*) 'Choice temp min max in which change is made' 513 read(*,*,iostat=ierr) val4,val5 if(ierr.ne.0) goto 513 write(*,*) 'Choice phisfi min max where change n2' 514 read(*,*,iostat=ierr) val6,val11 if(ierr.ne.0) goto 514 DO ig=1,ngridmx IF ( (latfi(ig)*180./pi.ge.val) .and. & (latfi(ig)*180./pi.le.val1) .and. & (lonfi(ig)*180./pi.ge.val2) .and. & (lonfi(ig)*180./pi.le.val3) .and. & (((tsurf(ig).ge.val4) .and. & (tsurf(ig).le.val5)) .or. & ((tsoil(ig,nsoilmx).ge.val4) .and. & (tsoil(ig,nsoilmx).le.val5))) .and. & (phisfi(ig).ge.val6) .and. & (phisfi(ig).le.val11) .and. & (qsurf(ig,igcm_n2).le.0.001) ) then ! DO i=1,ngridmx ! IF ( (latfi(i).eq.latfi(ig)) .and. ! & (lonfi(i).eq.0.) ) then ! tsurf(ig)=50. !qsurf(ig,igcm_ch4_ice)=qsurf(i,igcm_ch4_ice) ! DO l=1,nsoilmx tsoil(ig,l)=50. !tsoil(i,l) ENDDO !ENDIF !ENDDO ENDIF ENDDO c asymtriton : changing ice, tsurf and tsoil on triton c ---------------------------------------------------------- else if (modif(1:len_trim(modif)) .eq. 'asymtriton') then write(*,*) 'selecting area where tsoil to be check ' write(*,*) 'Choice center latitude of sinuisoid distrib?' 531 read(*,*,iostat=ierr) val1 if(ierr.ne.0) goto 531 write(*,*) 'Choice maginutde in latitude?' 532 read(*,*,iostat=ierr) val2 if(ierr.ne.0) goto 532 write(*,*) 'Choice center longitude?' 533 read(*,*,iostat=ierr) val3 if(ierr.ne.0) goto 533 ! write(*,*) 'iip1,jjp1=',iip1,jjp1,ngridmx DO ig=1,ngridmx ! Latitude of the sinusoid: val11=val1+val2*cos(lonfi(ig)*1.57079633*2/pi-val3*pi/180.) ! If above line and ice: remove ice IF ( (latfi(ig)*180./pi.ge.val11) .and. & (latfi(ig)*180./pi.le.val1+val2) .and. & (qsurf(ig,igcm_n2).gt.0.) ) then ! Look for same longitude point where no ice val5=1. jref=ig ! 1 ! ... iip1 ... x (jjp1-2) 32 x 23 ! 1 do while (val5.ge.1..and.jref.gt.iip1) ! northward point jref=jref-iip1+1 ! ice at the point val5=qsurf(jref,igcm_n2) ! write(*,*) jref, ! & latfi(jref)*180./pi,lonfi(jref)*180/pi,val5 enddo if (val5.ge.1) write(*,*) 'NO POINT FOUND WITH NO ICE' ! Copy point in the selected area tsurf(ig)=tsurf(jref) qsurf(ig,igcm_n2)=qsurf(jref,igcm_n2) qsurf(ig,igcm_ch4_ice)=qsurf(jref,igcm_ch4_ice) qsurf(ig,igcm_co_ice)=qsurf(jref,igcm_co_ice) DO l=1,nsoilmx tsoil(ig,l)=tsoil(jref,l) ENDDO ENDIF ! If below line and no ice: add ice IF ( (latfi(ig)*180./pi.le.val11) .and. & (latfi(ig)*180./pi.le.val1+val2) .and. & (qsurf(ig,igcm_n2).eq.0.) ) then ! Look for same longitude point where ice val5=1. jref=ig do while (val5.le.1.and.jref.lt.ngridmx-iip1) ! southward point jref=jref+iip1-1 ! ice at the point val5=qsurf(jref,igcm_n2) write(*,*) jref, & latfi(jref)*180./pi,lonfi(jref)*180/pi,val5 enddo if (val5.le.1) write(*,*) 'NO POINT FOUND WITH ICE' ! Copy point in the selected area tsurf(ig)=tsurf(jref) qsurf(ig,igcm_n2)=qsurf(jref,igcm_n2) qsurf(ig,igcm_ch4_ice)=qsurf(jref,igcm_ch4_ice) qsurf(ig,igcm_co_ice)=qsurf(jref,igcm_co_ice) DO l=1,nsoilmx tsoil(ig,l)=tsoil(jref,l) ENDDO ENDIF ENDDO c source_ch4 : adding source ch4 c ---------------------------------------------------------- else if (modif(1:len_trim(modif)) .eq. 'source_ch4') then write(*,*) 'Adding ch4 ice latitudinally ' write(*,*) 'Choice : lat1 and lat2 ?' 433 read(*,*,iostat=ierr) val,val2 if(ierr.ne.0) goto 433 write(*,*) 'Choice : lon1 and lon2 ?' 434 read(*,*,iostat=ierr) val3,val4 if(ierr.ne.0) goto 434 write(*,*) 'Choice amount (kg/m2)' 435 read(*,*,iostat=ierr) val5 if(ierr.ne.0) goto 435 DO ig=1,ngridmx IF ( (latfi(ig)*180./pi.ge.val) .and. & (latfi(ig)*180./pi.le.val2) .and. & (lonfi(ig)*180./pi.ge.val3) .and. & (lonfi(ig)*180./pi.lt.val4) ) then qsurf(ig,igcm_ch4_ice)=qsurf(ig,igcm_ch4_ice)+val5 ENDIF ENDDO c source_co : adding source co c ---------------------------------------------------------- else if (modif(1:len_trim(modif)) .eq. 'source_co') then write(*,*) 'Adding co ice latitudinally ' write(*,*) 'Choice : lat1 and lat2 ?' 436 read(*,*,iostat=ierr) val,val2 if(ierr.ne.0) goto 436 write(*,*) 'Choice : lon1 and lon2 ?' 437 read(*,*,iostat=ierr) val3,val4 if(ierr.ne.0) goto 437 write(*,*) 'Choice amount (kg/m2)' 438 read(*,*,iostat=ierr) val5 if(ierr.ne.0) goto 438 DO ig=1,ngridmx IF ( (latfi(ig)*180./pi.ge.val) .and. & (latfi(ig)*180./pi.le.val2) .and. & (lonfi(ig)*180./pi.ge.val3) .and. & (lonfi(ig)*180./pi.lt.val4) ) then qsurf(ig,igcm_co_ice)=qsurf(ig,igcm_co_ice)+val5 ENDIF ENDDO ! therm_ini_s: (re)-set soil thermal inertia to reference surface values ! ---------------------------------------------------------------------- else if (modif(1:len_trim(modif)).eq.'therm_ini_s') then ! write(*,*)"surfithfi(1):",surfithfi(1) do isoil=1,nsoilmx inertiedat(1:ngridmx,isoil)=surfithfi(1:ngridmx) enddo write(*,*)'OK: Soil thermal inertia has been reset to referenc &e surface values' write(*,*)"inertiedat(1,1):",inertiedat(1,1) ithfi(:,:)=inertiedat(:,:) ! recast ithfi() onto ith() call gr_fi_dyn(nsoilmx,ngridmx,iip1,jjp1,ithfi,ith) ! inert3d: set soil thermal inertia to specific uniform value ! ---------------------------------------------------------------------- else if (modif(1:len_trim(modif)).eq.'inert3d') then write(*,*) 'Actual value of surf Thermal inertia at ig=1: ', & inertiedat(1,1), ' SI' write(*,*) 'Value of Thermal inertia (SI) ? ' read(*,*) val do isoil=1,nsoilmx do ig=1,ngridmx inertiedat(ig,isoil)=val enddo enddo write(*,*)'OK: Soil TI set to ',val,' SI everywhere' ithfi(:,:)=inertiedat(:,:) ! recast ithfi() onto ith() call gr_fi_dyn(nsoilmx,ngridmx,iip1,jjp1,ithfi,ith) ! subsoilice_n: Put deep ice layer in northern hemisphere soil ! ------------------------------------------------------------ else if (modif(1:len_trim(modif)).eq.'subsoilice_n') then write(*,*)'From which latitude (in deg.), up to the north pole, &should we put subterranean ice?' ierr=1 do while (ierr.ne.0) read(*,*,iostat=ierr) val if (ierr.eq.0) then ! got a value ! do a sanity check if((val.lt.0.).or.(val.gt.90)) then write(*,*)'Latitude should be between 0 and 90 deg. !!!' ierr=1 else ! find corresponding jref (nearest latitude) ! note: rlatu(:) contains decreasing values of latitude ! starting from PI/2 to -PI/2 do j=1,jjp1 if ((rlatu(j)*180./pi.ge.val).and. & (rlatu(j+1)*180./pi.le.val)) then ! find which grid point is nearest to val: if (abs(rlatu(j)*180./pi-val).le. & abs((rlatu(j+1)*180./pi-val))) then jref=j else jref=j+1 endif write(*,*)'Will use nearest grid latitude which is:', & rlatu(jref)*180./pi endif enddo ! of do j=1,jjp1 endif ! of if((val.lt.0.).or.(val.gt.90)) endif !of if (ierr.eq.0) enddo ! of do while ! Build layers() (as in soil_settings.F) val2=sqrt(mlayer(0)*mlayer(1)) val3=mlayer(1)/mlayer(0) do isoil=1,nsoilmx layer(isoil)=val2*(val3**(isoil-1)) enddo write(*,*)'At which depth (in m.) does the ice layer begin?' write(*,*)'(currently, the deepest soil layer extends down to:' & ,layer(nsoilmx),')' ierr=1 do while (ierr.ne.0) read(*,*,iostat=ierr) val2 ! write(*,*)'val2:',val2,'ierr=',ierr if (ierr.eq.0) then ! got a value, but do a sanity check if(val2.gt.layer(nsoilmx)) then write(*,*)'Depth should be less than ',layer(nsoilmx) ierr=1 endif if(val2.lt.layer(1)) then write(*,*)'Depth should be more than ',layer(1) ierr=1 endif endif enddo ! of do while ! find the reference index iref the depth corresponds to ! if (val2.lt.layer(1)) then ! iref=1 ! else do isoil=1,nsoilmx-1 if((val2.gt.layer(isoil)).and.(val2.lt.layer(isoil+1))) & then iref=isoil exit endif enddo ! thermal inertia of the ice: ierr=1 do while (ierr.ne.0) write(*,*)'What is the value of subterranean ice thermal inert &ia? (e.g.: 2000)' read(*,*,iostat=ierr)iceith enddo ! of do while ! recast ithfi() onto ith() call gr_fi_dyn(nsoilmx,ngridmx,iip1,jjp1,ithfi,ith) do j=1,jref ! write(*,*)'j:',j,'rlatu(j)*180./pi:',rlatu(j)*180./pi do i=1,iip1 ! loop on longitudes ! Build "equivalent" thermal inertia for the mixed layer ith(i,j,iref+1)=sqrt((layer(iref+1)-layer(iref))/ & (((val2-layer(iref))/(ith(i,j,iref)**2))+ & ((layer(iref+1)-val2)/(iceith)**2))) ! Set thermal inertia of lower layers do isoil=iref+2,nsoilmx ith(i,j,isoil)=iceith ! ice enddo enddo ! of do i=1,iip1 enddo ! of do j=1,jjp1 CALL gr_dyn_fi(nsoilmx,iip1,jjp1,ngridmx,ith,ithfi) ! subsoilice_s: Put deep ice layer in southern hemisphere soil ! ------------------------------------------------------------ else if (modif(1:len_trim(modif)).eq.'subsoilice_s') then write(*,*)'From which latitude (in deg.), down to the south pol &e, should we put subterranean ice?' ierr=1 do while (ierr.ne.0) read(*,*,iostat=ierr) val if (ierr.eq.0) then ! got a value ! do a sanity check if((val.gt.0.).or.(val.lt.-90)) then write(*,*)'Latitude should be between 0 and -90 deg. !!!' ierr=1 else ! find corresponding jref (nearest latitude) ! note: rlatu(:) contains decreasing values of latitude ! starting from PI/2 to -PI/2 do j=1,jjp1 if ((rlatu(j)*180./pi.ge.val).and. & (rlatu(j+1)*180./pi.le.val)) then ! find which grid point is nearest to val: if (abs(rlatu(j)*180./pi-val).le. & abs((rlatu(j+1)*180./pi-val))) then jref=j else jref=j+1 endif write(*,*)'Will use nearest grid latitude which is:', & rlatu(jref)*180./pi endif enddo ! of do j=1,jjp1 endif ! of if((val.lt.0.).or.(val.gt.90)) endif !of if (ierr.eq.0) enddo ! of do while ! Build layers() (as in soil_settings.F) val2=sqrt(mlayer(0)*mlayer(1)) val3=mlayer(1)/mlayer(0) do isoil=1,nsoilmx layer(isoil)=val2*(val3**(isoil-1)) enddo write(*,*)'At which depth (in m.) does the ice layer begin?' write(*,*)'(currently, the deepest soil layer extends down to:' & ,layer(nsoilmx),')' ierr=1 do while (ierr.ne.0) read(*,*,iostat=ierr) val2 ! write(*,*)'val2:',val2,'ierr=',ierr if (ierr.eq.0) then ! got a value, but do a sanity check if(val2.gt.layer(nsoilmx)) then write(*,*)'Depth should be less than ',layer(nsoilmx) ierr=1 endif if(val2.lt.layer(1)) then write(*,*)'Depth should be more than ',layer(1) ierr=1 endif endif enddo ! of do while ! find the reference index iref the depth corresponds to do isoil=1,nsoilmx-1 if((val2.gt.layer(isoil)).and.(val2.lt.layer(isoil+1))) & then iref=isoil exit endif enddo ! thermal inertia of the ice: ierr=1 do while (ierr.ne.0) write(*,*)'What is the value of subterranean ice thermal inert &ia? (e.g.: 2000)' read(*,*,iostat=ierr)iceith enddo ! of do while ! recast ithfi() onto ith() call gr_fi_dyn(nsoilmx,ngridmx,iip1,jjp1,ithfi,ith) do j=jref,jjp1 ! write(*,*)'j:',j,'rlatu(j)*180./pi:',rlatu(j)*180./pi do i=1,iip1 ! loop on longitudes ! Build "equivalent" thermal inertia for the mixed layer ith(i,j,iref+1)=sqrt((layer(iref+1)-layer(iref))/ & (((val2-layer(iref))/(ith(i,j,iref)**2))+ & ((layer(iref+1)-val2)/(iceith)**2))) ! Set thermal inertia of lower layers do isoil=iref+2,nsoilmx ith(i,j,isoil)=iceith ! ice enddo enddo ! of do i=1,iip1 enddo ! of do j=jref,jjp1 CALL gr_dyn_fi(nsoilmx,iip1,jjp1,ngridmx,ith,ithfi) c ---------------------------------------------------------- c 'reservoir_SP' : increase or decrase ices reservoir in SP by factor c ---------------------------------------------------------- else if (modif(1:len_trim(modif)).eq.'reservoir_SP') then write(*,*) "Increase/Decrease N2 reservoir by factor :" read(*,*) val7 write(*,*) "Increase/Decrease CH4 reservoir by factor :" read(*,*) val8 write(*,*) "Increase/Decrease CO reservoir by factor :" read(*,*) val9 ! Definition SP: val3=-33. !lat1 val4=50. !lat2 val5=140. ! lon1 val6=-155. ! lon2 choice=-50. ! min phisfi write(*,*) 'def SP :' write(*,*) 'lat :',val3,val4 write(*,*) 'lon :',val5,'180 / -180',val6 write(*,*) 'choice phisfi min :',choice DO ig=1,ngridmx IF ( (latfi(ig)*180./pi.ge.val3) .and. & (latfi(ig)*180./pi.le.val4) .and. & ( ((lonfi(ig)*180./pi.ge.-180.) .and. & (lonfi(ig)*180./pi.le.val6)) .or. & ((lonfi(ig)*180./pi.ge.val5) .and. & (lonfi(ig)*180./pi.le.180.))) ) then IF ((phisfi(ig).le.choice)) then !.and. c & (qsurf(ig,igcm_n2).ge.50)) then qsurf(ig,igcm_n2)=qsurf(ig,igcm_n2)*val7 qsurf(ig,igcm_ch4_ice)=qsurf(ig,igcm_ch4_ice)*val8 qsurf(ig,igcm_co_ice)=qsurf(ig,igcm_co_ice)*val9 ENDIF ENDIF ENDDO c ---------------------------------------------------------- c 'reservoir' : increase or decrase ices reservoir by factor c ---------------------------------------------------------- else if (modif(1:len_trim(modif)).eq.'reservoir') then write(*,*) "Increase/Decrease N2 reservoir by factor :" read(*,*) val write(*,*) "Increase/Decrease CH4 reservoir by factor :" read(*,*) val2 write(*,*) "Increase/Decrease CO reservoir by factor :" read(*,*) val3 DO ig=1,ngridmx qsurf(ig,igcm_n2)=qsurf(ig,igcm_n2)*val qsurf(ig,igcm_ch4_ice)=qsurf(ig,igcm_ch4_ice)*val2 qsurf(ig,igcm_co_ice)=qsurf(ig,igcm_co_ice)*val3 ENDDO c -------------------------------------------------------- c 'plutomap' : initialize pluto ices distribution c -------------------------------------------------------- else if (modif(1:len_trim(modif)).eq.'plutomap') then write(*,*) 'pluto_map.dat has to be in your simulation file !!' open(27,file='pluto_map.dat') N2_pluto_dat(:,:) =0. CH4_pluto_dat(:,:) =0. CO_pluto_dat(:,:) =0. ! Dimension file pluto_map IF (jm_plu.ne.180) then write(*,*) 'Newstart : you should set jm_plu to 180' stop ENDIF ! Read values do j=1,jm_plu+1 read(27,*) (map_pluto_dat(i,j),i=1,im_plu) do i=1,im_plu !!!!!! BTD_v2 if (map_pluto_dat(i,j).eq.3) then N2_pluto_dat(i,j)=100000. elseif (map_pluto_dat(i,j).eq.4) then CH4_pluto_dat(i,j)=100000. elseif (map_pluto_dat(i,j).eq.0) then alb_dat(i,j)=0.3 elseif (map_pluto_dat(i,j).eq.6) then alb_dat(i,j)=0.6 elseif (map_pluto_dat(i,j).eq.7) then alb_dat(i,j)=0.1 endif !!!!!! BTD !if (map_pluto_dat(i,j).eq.3) then ! CH4_pluto_dat(i,j)=100000. !endif end do end do close (27) ! Interpolate !! latitude and longitude in REindexed pluto map : latv_plu(1)=90. *pi/180. do j=2,jm_plu latv_plu(j)= (90-(j-1 -0.5)*180./(jm_plu-1))*pi/180. end do latv_plu(jm_plu+1)= -90. *pi/180. do i=1,im_plu+1 lonu_plu(i)=(-180+ (i-1)*360./(im_plu) ) *pi/180. enddo !Reindexation to shift the longitude grid like a LMD GCM grid... do j=1,jm_plu+1 N2_pluto_rein(1,j)=(N2_pluto_dat(1,j)+ & N2_pluto_dat(im_plu,j))/2 CH4_pluto_rein(1,j)=(CH4_pluto_dat(1,j)+ & CH4_pluto_dat(im_plu,j))/2 CO_pluto_rein(1,j)=(CO_pluto_dat(1,j)+ & CO_pluto_dat(im_plu,j))/2 alb_rein(1,j)=(alb_dat(1,j)+ & alb_dat(im_plu,j))/2 do i=2,im_plu N2_pluto_rein(i,j)=(N2_pluto_dat(i-1,j)+ & N2_pluto_dat(i,j))/2 CH4_pluto_rein(i,j)=(CH4_pluto_dat(i-1,j)+ & CH4_pluto_dat(i,j))/2 CO_pluto_rein(i,j)=(CO_pluto_dat(i-1,j)+ & CO_pluto_dat(i,j))/2 alb_rein(i,j)=(alb_dat(i-1,j)+ & alb_dat(i,j))/2 end do N2_pluto_rein(im_plu+1,j) = N2_pluto_rein(1,j) CH4_pluto_rein(im_plu+1,j) = CH4_pluto_rein(1,j) CO_pluto_rein(im_plu+1,j) = CO_pluto_rein(1,j) alb_rein(im_plu+1,j) = alb_rein(1,j) end do call interp_horiz(N2_pluto_rein,N2_pluto_gcm, & im_plu,jm_plu,iim,jjm,1,lonu_plu,latv_plu,rlonu,rlatv) call interp_horiz(CH4_pluto_rein,CH4_pluto_gcm, & im_plu,jm_plu,iim,jjm,1,lonu_plu,latv_plu,rlonu,rlatv) call interp_horiz(CO_pluto_rein,CO_pluto_gcm, & im_plu,jm_plu,iim,jjm,1,lonu_plu,latv_plu,rlonu,rlatv) call interp_horiz(alb_rein,alb_gcm, & im_plu,jm_plu,iim,jjm,1,lonu_plu,latv_plu,rlonu,rlatv) ! Switch dyn grid to fi grid qsurf_tmp(:,:) =0. CALL gr_dyn_fi(1,iip1,jjp1,ngridmx, & N2_pluto_gcm,qsurf_tmp(1,igcm_n2)) if (igcm_ch4_ice.ne.0) then CALL gr_dyn_fi(1,iip1,jjp1,ngridmx, & CH4_pluto_gcm,qsurf_tmp(1,igcm_ch4_ice)) endif if (igcm_co_ice.ne.0) then CALL gr_dyn_fi(1,iip1,jjp1,ngridmx, & CO_pluto_gcm,qsurf_tmp(1,igcm_co_ice)) endif alb=alb_gcm CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,alb_gcm,albfi) !print*, 'N2dat=',N2_pluto_gcm !print*, 'COdat=',CO_pluto_gcm !print*, 'CH4dat=',CH4_pluto_gcm print*, 'N2dat=',qsurf_tmp(:,igcm_n2) print*, 'COdat=',qsurf_tmp(:,igcm_co_ice) print*, 'CH4dat=',qsurf_tmp(:,igcm_ch4_ice) ! Fill DO ig=1,ngridmx qsurf(ig,igcm_n2)=qsurf(ig,igcm_n2)+qsurf_tmp(ig,igcm_n2) if (igcm_ch4_ice.ne.0) then qsurf(ig,igcm_ch4_ice)=qsurf(ig,igcm_ch4_ice)+ & qsurf_tmp(ig,igcm_ch4_ice) endif if (igcm_co_ice.ne.0) then qsurf(ig,igcm_co_ice)=qsurf(ig,igcm_co_ice)+ & qsurf_tmp(ig,igcm_co_ice) endif ENDDO endif enddo ! of do ! infinite loop on liste of changes 999 continue c======================================================================= c Correct pressure on the new grid (menu 0) c======================================================================= if ((choix_1.eq.0).and.(.not.flagps0)) then r = 1000.*8.31/mugaz do j=1,jjp1 do i=1,iip1 ps(i,j) = ps(i,j) * . exp((phisold_newgrid(i,j)-phis(i,j)) / . (t(i,j,1) * r)) end do end do c periodicite de ps en longitude do j=1,jjp1 ps(1,j) = ps(iip1,j) end do end if c======================================================================= c Initialisation de la physique / ecriture de newstartfi : c======================================================================= CALL inifilr CALL pression(ip1jmp1, ap, bp, ps, p3d) c----------------------------------------------------------------------- c Initialisation pks: c----------------------------------------------------------------------- CALL exner_hyb(ip1jmp1, ps, p3d, pks, pk, pkf) ! Calcul de la temperature potentielle teta if (flagtset) then DO l=1,llm DO j=1,jjp1 DO i=1,iim teta(i,j,l) = Tset(i,j,l) * cpp/pk(i,j,l) ENDDO teta (iip1,j,l)= teta (1,j,l) ENDDO ENDDO else if (choix_1.eq.0) then DO l=1,llm DO j=1,jjp1 DO i=1,iim teta(i,j,l) = t(i,j,l) * cpp/pk(i,j,l) ENDDO teta (iip1,j,l)= teta (1,j,l) ENDDO ENDDO endif C Calcul intermediaire c if (choix_1.eq.0) then CALL massdair( p3d, masse ) c print *,' ALPHAX ',alphax DO l = 1, llm DO i = 1, iim xppn(i) = aire( i, 1 ) * masse( i , 1 , l ) xpps(i) = aire( i,jjp1 ) * masse( i , jjp1 , l ) ENDDO xpn = SUM(xppn)/apoln xps = SUM(xpps)/apols DO i = 1, iip1 masse( i , 1 , l ) = xpn masse( i , jjp1 , l ) = xps ENDDO ENDDO endif phis(iip1,:) = phis(1,:) itau=0 if (choix_1.eq.0) then day_ini=int(date) endif c CALL geopot ( ip1jmp1, teta , pk , pks, phis , phi ) CALL caldyn0( itau,ucov,vcov,teta,ps,masse,pk,phis , * phi,w, pbaru,pbarv,day_ini+time ) CALL dynredem0("restart.nc",day_ini,phis) CALL dynredem1("restart.nc",0.0,vcov,ucov,teta,q,masse,ps) C C Ecriture etat initial physique C call physdem0("restartfi.nc",longitude,latitude,nsoilmx,ngridmx, & llm, & nqtot,dtphys,real(day_ini),0.0, & cell_area,albfi,ithfi,zmea,zstd,zsig,zgam,zthe) call physdem1("restartfi.nc",nsoilmx,ngridmx,llm,nqtot, & dtphys,real(day_ini), & tsurf,tsoil,emis,q2,qsurf) ! & cloudfrac,totalfrac,hice, ! & rnat,pctsrf_sic,tslab,tsea_ice,sea_ice) c======================================================================= c Formats c======================================================================= 1 FORMAT(//10x,'la valeur de im =',i4,2x,'lue sur le fichier de dema *rrage est differente de la valeur parametree iim =',i4//) 2 FORMAT(//10x,'la valeur de jm =',i4,2x,'lue sur le fichier de dema *rrage est differente de la valeur parametree jjm =',i4//) 3 FORMAT(//10x,'la valeur de lllm =',i4,2x,'lue sur le fichier demar *rage est differente de la valeur parametree llm =',i4//) write(*,*) "newstart: All is well that ends well." end !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! function dist_pluto(lat1,lon1,lat2,lon2) implicit none real dist_pluto real lat1,lon1,lat2,lon2 real dlon, dlat real a,c real rad rad=1190. ! km dlon = lon2 - lon1 dlat = lat2 - lat1 a = (sin(dlat/2))**2 + cos(lat1) * cos(lat2) * (sin(dlon/2))**2 c = 2 * atan2( sqrt(a), sqrt(1-a) ) dist_pluto = rad * c return end