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,nvarid_input,nvarid_inputs INTEGER nid_fi_input 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) REAL field_input(ngridmx) REAL,ALLOCATABLE :: field_inputs(:,:) 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 integer randtab(ngridmx) 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 tsurf_bb,tsurf_bb2 real ptoto,pcap,patm,airetot,ptotn,patmn,psea,geop real tempsoil(24),levsoil(24) 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(field_inputs(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 ptot : Modification of the total pressure: ice + current atmosphere c ------------------------------------------------------------------- else if (modif(1:len_trim(modif)).eq.'ptot') then c calcul de la pression totale glace + atm actuelle patm=0. airetot=0. pcap=0. DO j=1,jjp1 DO i=1,iim ig=1+(j-2)*iim +i if(j.eq.1) ig=1 if(j.eq.jjp1) ig=ngridmx patm = patm + ps(i,j)*aire(i,j) airetot= airetot + aire(i,j) ENDDO ENDDO ptoto = pcap + patm print*,'Current total pressure at surface ', & ptoto/airetot print*,'new value?' read(*,*) ptotn ptotn=ptotn*airetot patmn=ptotn-pcap print*,'ptoto,patm,ptotn,patmn' print*,ptoto,patm,ptotn,patmn print*,'Mult. factor for pressure (atm only)', patmn/patm do j=1,jjp1 do i=1,iip1 ps(i,j)=ps(i,j)*patmn/patm enddo enddo c 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 end do if (yes.eq.'y') then write(*,*) 'OK : conservation of tracer total mass' DO iq =1, nqtot DO l=1,llm DO j=1,jjp1 DO i=1,iip1 q(i,j,l,iq)=q(i,j,l,iq)*patm/patmn ENDDO ENDDO ENDDO ENDDO else write(*,*) 'OK : conservation of tracer mixing ratio' end if 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 c write(*,*) 'SURFACE value of tracer ',trim(tname(iq)), c & ' ? (kg/m2)' c read(*,*) val c DO ig=1,ngridmx c qsurf(ig,iq)=val c 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).gt.0.001) ) then ! DO i=1,ngridmx ! IF ( (latfi(i).eq.latfi(ig)) .and. ! & (lonfi(i).eq.0.) ) then ! tsurf(ig)=34.7 !qsurf(ig,igcm_ch4_ice)=qsurf(i,igcm_ch4_ice) ! DO l=1,nsoilmx tsoil(ig,l)=34.7 !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 subsoil_SP : choice of thermal inertia values for SP c ---------------------------------------------------------------- else if (modif(1:len_trim(modif)) .eq. 'subsoil_SP') then write(*,*) 'New value for subsoil thermal inertia in SP ?' 510 read(*,*,iostat=ierr) ith_bb if(ierr.ne.0) goto 510 write(*,*) 'thermal inertia (new value):',ith_bb write(*,*)'At which depth (in m.) does the ice layer begin?' write(*,*)'(currently, the deepest soil layer extends down to:' & ,layer(1),' - ',layer(nsoilmx),')' write(*,*)'write 0 for uniform value for all subsurf levels?' 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 if(val2.eq.0) then write(*,*)'Thermal inertia set for all subsurface layers' ierr=0 else write(*,*)'Depth should be more than ',layer(1) ierr=1 endif endif endif enddo ! of do while ! find the reference index iref the depth corresponds to if(val2.eq.0) then iref=1 write(*,*)'Level selected is first level: ',layer(iref),' m' else do isoil=1,nsoilmx-1 if((val2.gt.layer(isoil)).and.(val2.lt.layer(isoil+1))) & then iref=isoil write(*,*)'Level selected : ',layer(isoil),' m' exit endif enddo endif ! Definition SP: val3=-50. !lat1 val4=60. !lat2 val5=130. ! lon1 val6=-140. ! lon2 choice=-100. ! 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) .and. & (qsurf(ig,igcm_n2).ge.1000)) then DO j=iref,nsoilmx ithfi(ig,j)=ith_bb ENDDO ENDIF ENDIF ENDDO CALL gr_fi_dyn(nsoilmx,ngridmx,iip1,jjp1,ithfi,ith) c topo_sp : change topo of SP c ----------------------------------------------------- else if (modif(1:len_trim(modif)) .eq. 'topo_sp') then ! def SP: val2=-33. !lat1 val3=50. !lat2 val4=140. ! lon1 val5=-155. ! lon2 write(*,*) 'def SP :' write(*,*) 'lat :',val2,val3 write(*,*) 'lon :',val4,'180 / -180',val5 write(*,*) 'choice phisfi min (ex: -100):' 606 read(*,*,iostat=ierr) val6 if(ierr.ne.0) goto 606 write(*,*) 'choice coefficient for phisfi (ex: 2):' 605 read(*,*,iostat=ierr) choice if(ierr.ne.0) goto 605 write(*,*) 're scale the pressure :' r = 8.314511E+0 *1000.E+0/mugaz temp=40. write(*,*) 'r, sale height temperature =',r,temp do j=1,jjp1 do i=1,iip1 phisprev= phis(i,j) IF ( (rlatu(j)*180./pi.ge.val2) .and. & (rlatu(j)*180./pi.le.val3) .and. & ( ((rlonv(i)*180./pi.ge.-180.) .and. & (rlonv(i)*180./pi.le.val5)) .or. & ((rlonv(i)*180./pi.ge.val4) .and. & (rlonv(i)*180./pi.le.180.))) ) then IF (phis(i,j).le.val6) then phis(i,j)=phis(i,j)*choice ps(i,j) = ps(i,j) * . exp((phisprev-phis(i,j))/(temp * r)) ENDIF ENDIF end do end do c periodicity of surface ps in longitude do j=1,jjp1 ps(1,j) = ps(iip1,j) end do CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,phis,phisfi) c fill_sp inv: fill inverted SP with N2 ice and adjust phisfi (flat SP) c ----------------------------------------------------- else if (modif(1:len_trim(modif)) .eq. 'fill_sp_inv') then ! def SP: val2=-33. !lat1 val3=50. !lat2 val4=-40. ! lon1 val5=25. ! lon2 write(*,*) 'def SP :' write(*,*) 'lat :',val2,val3 write(*,*) 'lon :',val4,val5 write(*,*) 'choice phisfi ref SP (ex: -800):' 706 read(*,*,iostat=ierr) val6 if(ierr.ne.0) goto 706 write(*,*) 're scale the pressure :' r = 8.314511E+0 *1000.E+0/mugaz temp=40. write(*,*) 'r, sale height temperature =',r,temp,g qsurf=0. write(*,*) 'latfi=',latfi !CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,phis,phisfi) write(*,*) 'phis=',phis write(*,*) 'phisfi=',phisfi !CALL gr_fi_dyn(1,ngridmx,iip1,jjp1,phisfi,phis) phisold = phis write(*,*) 'phisold=',phisold do ig=1,ngridmx write(*,*) 'lat=',latfi(ig)*180./pi write(*,*) 'lon=',lonfi(ig)*180./pi write(*,*) 'phisfi=',phisfi(ig) IF ( (latfi(ig)*180./pi.ge.val2) .and. & (latfi(ig)*180./pi.le.val3) .and. & (lonfi(ig)*180./pi.ge.val4) .and. & (lonfi(ig)*180./pi.le.val5) ) then write(*,*) 'hello SP',phisfi(ig),ig IF (phisfi(ig).le.val6) then qsurf(ig,igcm_n2)=(val6-phisfi(ig))/g*1000. phisfi(ig)=val6 write(*,*) 'changes',phisfi(ig),qsurf(ig,igcm_n2) ENDIF ENDIF end do c update new phis and ps CALL gr_fi_dyn(1,ngridmx,iip1,jjp1,phisfi,phis) do j=1,jjp1 do i=1,iip1 ps(i,j) = ps(i,j) * . exp((phisold(i,j)-phis(i,j))/(temp * r)) enddo enddo c periodicity of surface ps in longitude do j=1,jjp1 ps(1,j) = ps(iip1,j) end do c adjust phisfi according to the amount of N2 ice c ----------------------------------------------------- else if (modif(1:len_trim(modif)) .eq. 'adjustphi') then phisold = phis do ig=1,ngridmx phisfi(ig)=phisfi(ig)+qsurf(ig,igcm_n2)*g/1000. end do c update new phis and ps CALL gr_fi_dyn(1,ngridmx,iip1,jjp1,phisfi,phis) r = 8.314511E+0 *1000.E+0/mugaz temp=37. do j=1,jjp1 do i=1,iip1 ps(i,j) = ps(i,j) * . exp((phisold(i,j)-phis(i,j))/(temp * r)) enddo enddo c periodicity of surface ps in longitude do j=1,jjp1 ps(1,j) = ps(iip1,j) end do c correct phisfi c ----------------------------------------------------- else if (modif(1:len_trim(modif)) .eq. 'correctphi') then write(*,*) 'choice latitudes:' 537 read(*,*,iostat=ierr) val1,val2 if(ierr.ne.0) goto 537 write(*,*) 'choice longitudes:' 538 read(*,*,iostat=ierr) val3,val4 if(ierr.ne.0) goto 538 write(*,*) 'choice phi min max' 539 read(*,*,iostat=ierr) val5,val6 if(ierr.ne.0) goto 539 write(*,*) 'choice factor phi' 535 read(*,*,iostat=ierr) val8 if(ierr.ne.0) goto 535 write(*,*) 'choice add phi' 536 read(*,*,iostat=ierr) val7 if(ierr.ne.0) goto 536 do ig=1,ngridmx IF ( (latfi(ig)*180./pi.ge.val1) .and. & (latfi(ig)*180./pi.le.val2) .and. & (lonfi(ig)*180./pi.ge.val3) .and. & (lonfi(ig)*180./pi.le.val4) ) then IF ((phisfi(ig).le.val6).and. & (phisfi(ig).ge.val5)) then phisfi(ig)=phisfi(ig)*val8 phisfi(ig)=phisfi(ig)+val7 ENDIF ENDIF enddo phisold = phis c update new phis and ps CALL gr_fi_dyn(1,ngridmx,iip1,jjp1,phisfi,phis) r = 8.314511E+0 *1000.E+0/mugaz temp=37. do j=1,jjp1 do i=1,iip1 ps(i,j) = ps(i,j) * . exp((phisold(i,j)-phis(i,j))/(temp * r)) enddo enddo c periodicity of surface ps in longitude do j=1,jjp1 do i=1,iip1 ps(i,j) = ps(i,j) * . exp((phisold(i,j)-phis(i,j))/(temp * r)) enddo enddo c periodicity of surface ps in longitude do j=1,jjp1 ps(1,j) = ps(iip1,j) end do CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,phis,phisfi) c correct phisfi c ----------------------------------------------------- else if (modif(1:len_trim(modif)) .eq. 'correctps') then phisold = phis c update new phis and ps CALL gr_fi_dyn(1,ngridmx,iip1,jjp1,phisfi,phis) r = 8.314511E+0 *1000.E+0/mugaz temp=37. do j=1,jjp1 do i=1,iip1 ps(i,j) = ps(i,j) * . exp((phisold(i,j)-phis(i,j))/(temp * r)) enddo enddo c periodicity of surface ps in longitude do j=1,jjp1 do i=1,iip1 ps(i,j) = ps(i,j) * . exp((phisold(i,j)-phis(i,j))/(temp * r)) enddo enddo c periodicity of surface ps in longitude do j=1,jjp1 ps(1,j) = ps(iip1,j) end do CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,phis,phisfi) c No flat topo c ----------------------------------------------------- else if (modif(1:len_trim(modif)) .eq. 'toponoise') then write(*,*) 'choice latitudes:' 587 read(*,*,iostat=ierr) val1,val2 if(ierr.ne.0) goto 587 write(*,*) 'choice longitudes:' 588 read(*,*,iostat=ierr) val3,val4 if(ierr.ne.0) goto 588 write(*,*) 'choice phi min max' 589 read(*,*,iostat=ierr) val5,val6 if(ierr.ne.0) goto 589 write(*,*) 'choice amplitude min max new phi' 585 read(*,*,iostat=ierr) val7,val8 if(ierr.ne.0) goto 585 c write(*,*) 'choice nb of random cases' c586 read(*,*,iostat=ierr) val7 c if(ierr.ne.0) goto 586 do ig=1,ngridmx IF ( (latfi(ig)*180./pi.ge.val1) .and. & (latfi(ig)*180./pi.le.val2) .and. & (lonfi(ig)*180./pi.ge.val3) .and. & (lonfi(ig)*180./pi.le.val4) ) then IF ((phisfi(ig).le.val6).and. & (phisfi(ig).ge.val5)) then CALL RANDOM_NUMBER(val9) phisfi(ig)=val7+(val8-val7)*val9 ENDIF ENDIF enddo phisold = phis c update new phis and ps CALL gr_fi_dyn(1,ngridmx,iip1,jjp1,phisfi,phis) r = 8.314511E+0 *1000.E+0/mugaz temp=37. do j=1,jjp1 do i=1,iip1 ps(i,j) = ps(i,j) * . exp((phisold(i,j)-phis(i,j))/(temp * r)) enddo enddo c periodicity of surface ps in longitude do j=1,jjp1 do i=1,iip1 ps(i,j) = ps(i,j) * . exp((phisold(i,j)-phis(i,j))/(temp * r)) enddo enddo c periodicity of surface ps in longitude do j=1,jjp1 ps(1,j) = ps(iip1,j) end do CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,phis,phisfi) c fill_sp : fill SP with N2 ice and adjust phisfi (flat SP) c ----------------------------------------------------- else if (modif(1:len_trim(modif)) .eq. 'fill_sp') then ! def SP: val2=-33. !lat1 val3=50. !lat2 val4=140. ! lon1 val5=-155. ! lon2 write(*,*) 'def SP :' write(*,*) 'lat :',val2,val3 write(*,*) 'lon :',val4,'180 / -180',val5 write(*,*) 'choice phisfi ref SP (ex: -800):' 607 read(*,*,iostat=ierr) val6 if(ierr.ne.0) goto 607 write(*,*) 're scale the pressure :' r = 8.314511E+0 *1000.E+0/mugaz temp=40. write(*,*) 'r, sale height temperature =',r,temp,g qsurf=0. write(*,*) 'latfi=',latfi !CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,phis,phisfi) write(*,*) 'phis=',phis write(*,*) 'phisfi=',phisfi !CALL gr_fi_dyn(1,ngridmx,iip1,jjp1,phisfi,phis) phisold = phis write(*,*) 'phisold=',phisold do ig=1,ngridmx write(*,*) 'lat=',latfi(ig)*180./pi write(*,*) 'lon=',lonfi(ig)*180./pi write(*,*) 'phisfi=',phisfi(ig) IF ( (latfi(ig)*180./pi.ge.val2) .and. & (latfi(ig)*180./pi.le.val3) .and. & ( ((lonfi(ig)*180./pi.ge.-180.) .and. & (lonfi(ig)*180./pi.le.val5)) .or. & ((lonfi(ig)*180./pi.ge.val4) .and. & (lonfi(ig)*180./pi.le.180.))) ) then write(*,*) 'hello SP',phisfi(ig),ig IF (phisfi(ig).le.val6) then qsurf(ig,igcm_n2)=(val6-phisfi(ig))/g*1000. phisfi(ig)=val6 write(*,*) 'changes',phisfi(ig),qsurf(ig,igcm_n2) ENDIF ENDIF end do c update new phis and ps CALL gr_fi_dyn(1,ngridmx,iip1,jjp1,phisfi,phis) do j=1,jjp1 do i=1,iip1 ps(i,j) = ps(i,j) * . exp((phisold(i,j)-phis(i,j))/(temp * r)) enddo enddo c periodicity of surface ps in longitude do j=1,jjp1 ps(1,j) = ps(iip1,j) end do c bugch4 correct bug warm ch4 due to change of resolution c ----------------------------------------------------- else if (modif(1:len_trim(modif)) .eq. 'bugch4') then write(*,*) 'Ok we are going to try to solve this bug' write(*,*) 'First extract a profile of tsoil and tsurf' write(*,*) 'that you want to set as reference :' write(*,*) 'tsoil_ref and tsurf_ref ' write(*,*) 'number of points to check ' 546 read(*,*,iostat=ierr) jref1 if(ierr.ne.0) goto 546 write(*,*) 'choice acceptable tsurf range for ch4: t1,t2:' 547 read(*,*,iostat=ierr) val1,val2 if(ierr.ne.0) goto 547 ! Check tsurf at all locations DO ig=1,ngridmx IF (qsurf(ig,igcm_ch4_ice).gt.0.001.and. & qsurf(ig,igcm_n2).eq.0.) then IF ((tsurf(ig).lt.val1) .or. & (tsurf(ig).ge.val2)) then write(*,*) 'Pb tsurf point ',ig do jref=1,jref1 IF ((ig-jref.ge.1).and.qsurf(ig,igcm_n2).eq.0..and. & (qsurf(int(ig-jref),igcm_ch4_ice).gt.0.001).and. & (tsurf(ig-jref).gt.val1 & .and.tsurf(ig-jref).lt.val2)) then tsurf(ig)=tsurf(int(ig-jref)) DO l=1,nsoilmx tsoil(ig,l)=tsoil(int(ig-jref),l) ENDDO goto 548 ELSEIF ((ig+jref.le.ngridmx).and. & qsurf(ig,igcm_n2).eq.0..and. & (qsurf(int(ig+jref),igcm_ch4_ice).gt.0.001).and. & (tsurf(ig+jref).gt.val1 & .and.tsurf(ig+jref).lt.val2)) then tsurf(ig)=tsurf(int(ig+jref)) DO l=1,nsoilmx tsoil(ig,l)=tsoil(int(ig+jref),l) ENDDO goto 548 ENDIF enddo 548 write(*,*) 'found point with ch4' ENDIF ENDIF END DO c bugn2 correct bug warm n2 due to change of resolution c ----------------------------------------------------- else if (modif(1:len_trim(modif)) .eq. 'bugn2') then write(*,*) 'Ok we are going to try to solve this bug' write(*,*) 'First extract a profile of tsoil and tsurf' write(*,*) 'that you want to set as reference :' write(*,*) 'tsoil_ref and tsurf_ref ' write(*,*) 'number of points to check ' 544 read(*,*,iostat=ierr) jref1 if(ierr.ne.0) goto 544 write(*,*) 'choice acceptable tsurf range for n2: t1,t2:' 540 read(*,*,iostat=ierr) val1,val2 if(ierr.ne.0) goto 540 ! Check tsurf at all locations DO ig=1,ngridmx IF (qsurf(ig,1).gt.0.001) then IF ((tsurf(ig).lt.val1) .or. & (tsurf(ig).ge.val2)) then write(*,*) 'Pb tsurf point ',ig ! IF low topo et peu/pas de N2: add ch4, CO, N2 do jref=1,jref1 IF ((ig-jref.ge.1).and. & (qsurf(int(ig-jref),igcm_n2).gt.0.001).and. & (tsurf(ig-jref).gt.val1 & .and.tsurf(ig-jref).lt.val2)) then tsurf(ig)=tsurf(int(ig-jref)) DO l=1,nsoilmx tsoil(ig,l)=tsoil(int(ig-jref),l) ENDDO goto 541 ELSEIF ((ig+jref.le.ngridmx).and. & (qsurf(int(ig+jref),igcm_n2).gt.0.001).and. & (tsurf(ig+jref).gt.val1 & .and.tsurf(ig+jref).lt.val2)) then tsurf(ig)=tsurf(int(ig+jref)) DO l=1,nsoilmx tsoil(ig,l)=tsoil(int(ig+jref),l) ENDDO goto 541 ENDIF enddo 541 write(*,*) 'found point with n2' ENDIF ENDIF END DO ! second tour ! DO ig=1,ngridmx ! IF (qsurf(ig,1).gt.0.001) then ! IF ((tsurf(ig).lt.val1) .or. ! & (tsurf(ig).ge.val2)) then ! ! IF low topo et peu/pas de N2: add ch4, CO, N2 ! IF ((ig-1.lt.1).or.(ig+1.gt.ngridmx)) then ! write(*,*) 'pole : doing nothing' ! ELSEIF (qsurf(ig-1,igcm_n2).gt.0.001) then ! tsurf(ig)=tsurf(ig-1) ! DO l=1,nsoilmx ! tsoil(ig,l)=tsoil(ig-1,l) ! ENDDO ! ELSEIF (qsurf(ig+1,igcm_n2).gt.0.001) then ! tsurf(ig)=tsurf(ig+1) ! DO l=1,nsoilmx ! tsoil(ig,l)=tsoil(ig+1,l) ! ENDDO ! ENDIF ! ENDIF ! ENDIF ! END DO c flatnosp : flat topo outside a specific terrain (SP) c ----------------------------------------------------- else if (modif(1:len_trim(modif)) .eq. 'flatnosp') then write(*,*) 'Choice of topography (m) below which we keep ?' 551 read(*,*,iostat=ierr) val if(ierr.ne.0) goto 551 write(*,*) 'gravity g is : ',g geop=g*val write(*,*) 'Choice of minimum amount of N2 ice we keep ?' 552 read(*,*,iostat=ierr) val2 if(ierr.ne.0) goto 552 write(*,*) 'phis=',phis write(*,*) 'phisfi=',phisfi do ig=1,ngridmx IF ( (qsurf(ig,1).le.val2) .and. & (phisfi(ig).gt.geop) ) then write(*,*) 'hello SP',phisfi(ig),ig phisfi(ig)=0. ENDIF end do phisold = phis write(*,*) 're scale the pressure :' r = 8.314511E+0 *1000.E+0/mugaz temp=40. write(*,*) 'r, sale height temperature =',r,temp,g c update new phis and ps CALL gr_fi_dyn(1,ngridmx,iip1,jjp1,phisfi,phis) do j=1,jjp1 do i=1,iip1 ps(i,j) = ps(i,j) * . exp((phisold(i,j)-phis(i,j))/(temp * r)) enddo enddo c periodicity of surface ps in longitude do j=1,jjp1 ps(1,j) = ps(iip1,j) end do c flatregion : flat topo specific to region c ----------------------------------------------------- else if (modif(1:len_trim(modif)) .eq. 'flatregion') then write(*,*) 'Choice band : lat1 and lat2 ?' 553 read(*,*,iostat=ierr) val,val2 if(ierr.ne.0) goto 553 write(*,*) 'Choice band : lon1 and lon2 ?' 554 read(*,*,iostat=ierr) val3,val4 if(ierr.ne.0) goto 554 write(*,*) 'Choice of topography range ?' 555 read(*,*,iostat=ierr) val5,val6 if(ierr.ne.0) goto 555 write(*,*) 'Choice topo ?' 556 read(*,*,iostat=ierr) val7 if(ierr.ne.0) goto 556 write(*,*) 'phis=',phis write(*,*) 'phisfi=',phisfi 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) ) then IF ( (phisfi(ig).ge.val5) .and. & (phisfi(ig).le.val6) ) then write(*,*) 'hello topo',phisfi(ig),ig phisfi(ig)=val7 ENDIF ENDIF end do r = 8.314511E+0 *1000.E+0/mugaz temp=40. phisold = phis c update new phis and ps CALL gr_fi_dyn(1,ngridmx,iip1,jjp1,phisfi,phis) do j=1,jjp1 do i=1,iip1 ps(i,j) = ps(i,j) * . exp((phisold(i,j)-phis(i,j))/(temp * r)) enddo enddo c periodicity of surface ps in longitude do j=1,jjp1 ps(1,j) = ps(iip1,j) end do c changetopo c ----------------------------------------------------- else if (modif(1:len_trim(modif)) .eq. 'changetopo') then write(*,*) 'Choice band : lat1 and lat2 ?' 573 read(*,*,iostat=ierr) val,val2 if(ierr.ne.0) goto 573 write(*,*) 'Choice band : lon1 and lon2 ?' 574 read(*,*,iostat=ierr) val3,val4 if(ierr.ne.0) goto 574 write(*,*) 'Choice of topography range ?' 575 read(*,*,iostat=ierr) val5,val6 if(ierr.ne.0) goto 575 write(*,*) 'Choice change topo factor?' 576 read(*,*,iostat=ierr) val7 if(ierr.ne.0) goto 576 write(*,*) 'Choice change topo add?' 577 read(*,*,iostat=ierr) val8 if(ierr.ne.0) goto 577 write(*,*) 'phis=',phis write(*,*) 'phisfi=',phisfi 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) ) then IF ( (phisfi(ig).ge.val5) .and. & (phisfi(ig).le.val6) ) then write(*,*) 'hello topo',phisfi(ig),ig phisfi(ig)=phisfi(ig)*val7 phisfi(ig)=phisfi(ig)+val8 ENDIF ENDIF end do r = 8.314511E+0 *1000.E+0/mugaz temp=40. phisold = phis c update new phis and ps CALL gr_fi_dyn(1,ngridmx,iip1,jjp1,phisfi,phis) do j=1,jjp1 do i=1,iip1 ps(i,j) = ps(i,j) * . exp((phisold(i,j)-phis(i,j))/(temp * r)) enddo enddo c periodicity of surface ps in longitude do j=1,jjp1 ps(1,j) = ps(iip1,j) end do c corrtopo: corr topo specific bug c ----------------------------------------------------- else if (modif(1:len_trim(modif)) .eq. 'corrtopo') then ! get min max lon print*, latfi*180/pi print*, '***************************' print*, '***************************' print*, '***************************' print*, '***************************' print*, '***************************' print*, lonfi*180/pi print*, 'iip1=',iip1 do ig=2,ngridmx-1 IF (lonfi(ig)*180./pi.eq.-180.) THEN print*, lonfi(ig),lonfi(ig+iip1-2) phisfi(ig)=(phisfi(ig+1)+phisfi(ig+iip1))/2 ENDIF enddo phisold = phis r = 8.314511E+0 *1000.E+0/mugaz temp=40. c update new phis and ps CALL gr_fi_dyn(1,ngridmx,iip1,jjp1,phisfi,phis) do j=1,jjp1 do i=1,iip1 ps(i,j) = ps(i,j) * . exp((phisold(i,j)-phis(i,j))/(temp * r)) enddo enddo c periodicity of surface ps in longitude do j=1,jjp1 ps(1,j) = ps(iip1,j) end do c asymtopo: c ----------------------------------------------------- else if (modif(1:len_trim(modif)) .eq. 'asymtopo') then ! get diff altitude write(*,*) 'Diff N-S altitude in km (pos = N higher) ?' 591 read(*,*,iostat=ierr) val if(ierr.ne.0) goto 591 write(*,*) 'Coeff smoot topo (small = smoother) ?' 592 read(*,*,iostat=ierr) val2 if(ierr.ne.0) goto 592 do ig=1,ngridmx phisfi(ig)=phisfi(ig)+val*1000.*g*tanh(latfi(ig)*180/pi*val2) enddo phisold = phis r = 8.314511E+0 *1000.E+0/mugaz temp=40. c update new phis and ps CALL gr_fi_dyn(1,ngridmx,iip1,jjp1,phisfi,phis) do j=1,jjp1 do i=1,iip1 ps(i,j) = ps(i,j) * . exp((phisold(i,j)-phis(i,j))/(temp * r)) enddo enddo c periodicity of surface ps in longitude do j=1,jjp1 ps(1,j) = ps(iip1,j) end do c seticeNH : set the ices in the SP crater with NH topo c ----------------------------------------------------- else if (modif(1:len_trim(modif)) .eq. 'seticeNH') then open(333,file='./tsoil_180_30',form='formatted') do i=1,24 read(333,*) levsoil(i), tempsoil(i) enddo close(333) open(334,file='./tsurf_180_30',form='formatted') read(334,*) val close(334) write(*,*) 'Tsoil and tsurf ref are:' write(*,*) tempsoil write(*,*) 'tsurf=',val ! def SP: val2=-43. !lat1 val3=51. !lat2 val4=143. ! lon1 val5=-158. ! lon2 val6=-150 ! phisfi min ! Rm all N2 ice outside SP DO ig=1,ngridmx ! IF inside SP IF ( (latfi(ig)*180./pi.ge.val2) .and. & (latfi(ig)*180./pi.le.val3) .and. & ( ((lonfi(ig)*180./pi.ge.-180.) .and. & (lonfi(ig)*180./pi.le.val5)) .or. & ((lonfi(ig)*180./pi.ge.val4) .and. & (lonfi(ig)*180./pi.le.180.))) ) then ! IF low topo et peu/pas de N2: add ch4, CO, N2 IF ((phisfi(ig).le.val6).and. & (qsurf(ig,igcm_n2).le.500)) then qsurf(ig,igcm_n2)=1000. qsurf(ig,igcm_ch4_ice)=1000. qsurf(ig,igcm_co_ice)=1000. tsurf(ig)=val DO l=1,nsoilmx tsoil(ig,l)=tempsoil(l) ENDDO ENDIF ! IF high topo : rm N2 IF ( (qsurf(ig,igcm_n2).ge.20.).and. & (phisfi(ig).ge.val6) ) then qsurf(ig,igcm_n2)=0. ENDIF ! Mais on garde ch4 et CO en prenant la meme quantite ! qu'une autre latitude IF ( (qsurf(ig,igcm_ch4_ice).ge.20.) .and. & (phisfi(ig).ge.val6) ) then jref=int(ig/jjp1)+2 qsurf(ig,igcm_ch4_ice)=qsurf(jref,igcm_ch4_ice) ENDIF IF ( (qsurf(ig,igcm_co_ice).ge.20.) .and. & (phisfi(ig).ge.val6) ) then jref=int(ig/jjp1)+2 qsurf(ig,igcm_co_ice)=qsurf(jref,igcm_co_ice) ENDIF ELSE ! Outside SP IF (qsurf(ig,igcm_n2).ge.20.) then qsurf(ig,igcm_n2)=0. ENDIF IF (qsurf(ig,igcm_ch4_ice).ge.20.) then jref=int(ig/jjp1)+2 qsurf(ig,igcm_ch4_ice)=qsurf(jref,igcm_ch4_ice) ENDIF IF (qsurf(ig,igcm_co_ice).ge.20.) then jref=int(ig/jjp1)+2 qsurf(ig,igcm_co_ice)=qsurf(jref,igcm_co_ice) ENDIF ENDIF END DO c 'nomountain ' c -------------- else if (modif(1:len_trim(modif)) .eq. 'nomountain') then do j=1,jjp1 do i=1,iip1 if (phis(i,j).gt.0.1) then phis(i,j) = 0. ps(i,j)=ps(iim/4,j) ! assuming no topo here endif enddo enddo CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,phis,phisfi) c 'relief' c -------------- else if (modif(1:len_trim(modif)) .eq.'relief') then write(*,*) "add a circular mountain/crater" write(*,*) 'Center: lat lon ?' 707 read(*,*,iostat=ierr) lat0, lon0 if(ierr.ne.0) goto 707 if(lon0.gt.180.) lon0=lon0-360. lat0 = lat0 * pi/180. lon0 = lon0 * pi/180. DO ig=1,ngridmx if(abs(latfi(ig)-lat0).lt.pi/float(jjm) ) then if(abs(lonfi(ig)-lon0).lt.2*pi/float(iim) ) then ig0 = ig end if end if END DO write(*,*) "Reference Point to be used:" write(*,*) 'ig0',ig0 write(*,*) 'lat=',latfi(ig0)*180./pi write(*,*) 'lon=',lonfi(ig0)*180./pi write(*,*) 'radius (km), height (km) negative if crater ?' 508 read(*,*,iostat=ierr) radm, height if(ierr.ne.0) goto 508 write(*,*) 'Sale height temperature (ex:38) ?' 509 read(*,*,iostat=ierr) temp if(ierr.ne.0) goto 509 r = 8.314511E+0 *1000.E+0/mugaz do j=1,jjp1 do i=1,iip1 dist= dist_pluto(lat0,lon0,rlatu(j),rlonv(i)) if (dist.le.radm) then phisprev= phis(i,j) phis(i,j)=phisprev+1000.*g*height*(radm-dist)/radm write(*,*) 'lat=',rlatu(j)*180./pi write(*,*) 'lon=',rlonv(i)*180./pi write(*,*) 'dist', dist write(*,*) 'z(m)=', phis(i,j)/g ps(i,j) = ps(i,j) * . exp((phis(i,j))/(temp * r)) end if end do end do c periodicity of surface ps in longitude do j=1,jjp1 ps(1,j) = ps(iip1,j) end do CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,phis,phisfi) c subsoil_n2 : choice of thermal inertia values for N2 only c ---------------------------------------------------------------- else if (modif(1:len_trim(modif)) .eq. 'subsoil_n2') then write(*,*) 'New value for subsoil thermal inertia ?' 102 read(*,*,iostat=ierr) ith_bb if(ierr.ne.0) goto 102 write(*,*) 'thermal inertia (new value):',ith_bb write(*,*)'At which depth (in m.) does the ice layer begin?' write(*,*)'(currently, the deepest soil layer extends down to:' & ,layer(1),' - ',layer(nsoilmx),')' write(*,*)'write 0 for uniform value for all subsurf levels?' 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 if(val2.eq.0) then write(*,*)'Thermal inertia set for all subsurface layers' ierr=0 else write(*,*)'Depth should be more than ',layer(1) ierr=1 endif endif endif enddo ! of do while ! find the reference index iref the depth corresponds to if(val2.eq.0) then iref=1 write(*,*)'Level selected is first level: ',layer(iref),' m' else do isoil=1,nsoilmx-1 if((val2.gt.layer(isoil)).and.(val2.lt.layer(isoil+1))) & then iref=isoil write(*,*)'Level selected : ',layer(isoil),' m' exit endif enddo endif DO ig=1,ngridmx DO j=iref,nsoilmx c if((qsurf(ig,igcm_ch4_ice).lt.1.).and. c & (qsurf(ig,igcm_co_ice).lt.1.))then if (qsurf(ig,igcm_n2).gt.0.001) then ithfi(ig,j)=ith_bb endif ENDDO ENDDO CALL gr_fi_dyn(nsoilmx,ngridmx,iip1,jjp1,ithfi,ith) c subsoil_ch4 : choice of thermal inertia values for CH4 only c ---------------------------------------------------------------- else if (modif(1:len_trim(modif)) .eq. 'subsoil_ch4') then write(*,*) 'New value for subsoil thermal inertia ?' 103 read(*,*,iostat=ierr) ith_bb if(ierr.ne.0) goto 103 write(*,*) 'thermal inertia (new value):',ith_bb write(*,*)'At which depth (in m.) does the ice layer begin?' write(*,*)'(currently, the deepest soil layer extends down to:' & ,layer(1),' - ',layer(nsoilmx),')' write(*,*)'write 0 for uniform value for all subsurf levels?' 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 if(val2.eq.0) then write(*,*)'Thermal inertia set for all subsurface layers' ierr=0 else write(*,*)'Depth should be more than ',layer(1) ierr=1 endif endif endif enddo ! of do while ! find the reference index iref the depth corresponds to if(val2.eq.0) then iref=1 write(*,*)'Level selected is first level: ',layer(iref),' m' else do isoil=1,nsoilmx-1 if((val2.gt.layer(isoil)).and.(val2.lt.layer(isoil+1))) & then iref=isoil write(*,*)'Level selected : ',layer(isoil),' m' exit endif enddo endif DO ig=1,ngridmx DO j=iref,nsoilmx if (qsurf(ig,igcm_ch4_ice).gt.0.001.and. & qsurf(ig,igcm_n2).lt.0.001) then ithfi(ig,j)=ith_bb endif ENDDO ENDDO CALL gr_fi_dyn(nsoilmx,ngridmx,iip1,jjp1,ithfi,ith) c subsoil_alb : choice of thermal inertia values from albedo c ---------------------------------------------------------------- else if (modif(1:len_trim(modif)) .eq. 'subsoil_alb') then write(*,*) 'Choice range albedo for defining TI' 145 read(*,*,iostat=ierr) val,val2 if(ierr.ne.0) goto 145 write(*,*) 'New value for subsoil thermal inertia ?' 144 read(*,*,iostat=ierr) ith_bb if(ierr.ne.0) goto 144 write(*,*) 'thermal inertia (new value):',ith_bb write(*,*)'At which depth (in m.) does the ice layer begin?' write(*,*)'(currently, the deepest soil layer extends down to:' & ,layer(1),' - ',layer(nsoilmx),')' write(*,*)'write 0 for uniform value for all subsurf levels?' ierr=1 do while (ierr.ne.0) read(*,*,iostat=ierr) val3 if (ierr.eq.0) then ! got a value, but do a sanity check if(val3.gt.layer(nsoilmx)) then write(*,*)'Depth should be less than ',layer(nsoilmx) ierr=1 endif if(val3.lt.layer(1)) then if(val3.eq.0) then write(*,*)'Thermal inertia set for all subsurface layers' ierr=0 else write(*,*)'Depth should be more than ',layer(1) ierr=1 endif endif endif enddo ! of do while ! find the reference index iref the depth corresponds to if(val3.eq.0) then iref=1 write(*,*)'Level selected is first level: ',layer(iref),' m' else do isoil=1,nsoilmx-1 if ((val3.gt.layer(isoil)).and.(val3.lt.layer(isoil+1))) & then iref=isoil+1 write(*,*)'Level selected : ',layer(isoil+1),' m' exit endif enddo endif DO ig=1,ngridmx IF ( (albfi(ig).ge.val) .and. & (albfi(ig).le.val2) ) then DO j=iref,nsoilmx ithfi(ig,j)=ith_bb ENDDO ENDIF ENDDO CALL gr_fi_dyn(nsoilmx,ngridmx,iip1,jjp1,ithfi,ith) c subsoil_all : choice of thermal inertia values (global) c ---------------------------------------------------------------- else if (modif(1:len_trim(modif)) .eq. 'subsoil_all') then write(*,*) 'New value for subsoil thermal inertia ?' 104 read(*,*,iostat=ierr) ith_bb if(ierr.ne.0) goto 104 write(*,*) 'thermal inertia (new value):',ith_bb write(*,*)'At which depth (in m.) does the ice layer begin?' write(*,*)'(currently, the deepest soil layer extends down to:' & ,layer(1),' - ',layer(nsoilmx),')' write(*,*)'write 0 for uniform value for all subsurf levels?' 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 if(val2.eq.0) then write(*,*)'Thermal inertia set for all subsurface layers' ierr=0 else write(*,*)'Depth should be more than ',layer(1) ierr=1 endif endif endif enddo ! of do while ! find the reference index iref the depth corresponds to if(val2.eq.0) then iref=1 write(*,*)'Level selected is first level: ',layer(iref),' m' else do isoil=1,nsoilmx-1 !write(*,*)'isoil, ',isoil,val2 !write(*,*)'lay(i),lay(i+1):',layer(isoil),layer(isoil+1),' m' if ((val2.gt.layer(isoil)).and.(val2.lt.layer(isoil+1))) & then iref=isoil+1 write(*,*)'Level selected : ',layer(isoil+1),' m' exit endif enddo endif DO ig=1,ngridmx DO j=iref,nsoilmx ithfi(ig,j)=ith_bb ENDDO ENDDO CALL gr_fi_dyn(nsoilmx,ngridmx,iip1,jjp1,ithfi,ith) c diurnal_TI : choice of thermal inertia values (global) c ---------------------------------------------------------------- else if (modif(1:len_trim(modif)) .eq. 'diurnal_TI') then write(*,*) 'New value for diurnal thermal inertia ?' 106 read(*,*,iostat=ierr) ith_bb if(ierr.ne.0) goto 106 write(*,*) 'Diurnal thermal inertia (new value):',ith_bb write(*,*)'At which depth (in m.) does the ice layer ends?' write(*,*)'(currently, the soil layer 1 and nsoil are:' & ,layer(1),' - ',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 !write(*,*)'isoil, ',isoil,val2 !write(*,*)'lay(i),lay(i+1):',layer(isoil),layer(isoil+1),' m' if ((val2.gt.layer(isoil)).and.(val2.lt.layer(isoil+1))) & then iref=isoil+1 write(*,*)'Level selected : ',layer(isoil+1),' m' exit endif enddo DO ig=1,ngridmx DO j=1,iref ithfi(ig,j)=ith_bb ENDDO ENDDO CALL gr_fi_dyn(nsoilmx,ngridmx,iip1,jjp1,ithfi,ith) c Choice surface temperature c ----------------------------------------------------- else if (modif(1:len_trim(modif)) .eq. 'initsurf') then write(*,*) 'New value for initial surface temperature ?' 105 read(*,*,iostat=ierr) tsurf_bb if(ierr.ne.0) goto 105 write(*,*) 'new surface temperature (K):',tsurf_bb DO ig=1,ngridmx tsurf(ig)=tsurf_bb ENDDO c Modify surface temperature (for GCM start) c ----------------------------------------------------- else if (modif(1:len_trim(modif)) .eq. 'modtsurf') then write(*,*) 'Choice band : lat1 and lat2 ?' 455 read(*,*,iostat=ierr) val,val2 if(ierr.ne.0) goto 455 write(*,*) 'Choice band : lon1 and lon2 ?' 456 read(*,*,iostat=ierr) val4,val5 if(ierr.ne.0) goto 456 write(*,*) 'Choice topo min max ' 473 read(*,*,iostat=ierr) val9,val10 if(ierr.ne.0) goto 473 write(*,*) 'Choice tsurf min max ' 474 read(*,*,iostat=ierr) val11,val12 if(ierr.ne.0) goto 474 write(*,*) 'Choice multiplicative factor' 457 read(*,*,iostat=ierr) val3 if(ierr.ne.0) goto 457 write(*,*) 'Choice additional tsurf' 458 read(*,*,iostat=ierr) val6 if(ierr.ne.0) goto 458 write(*,*) 'Choice multiplicative factor tsoil' 459 read(*,*,iostat=ierr) val7 if(ierr.ne.0) goto 459 write(*,*) 'Choice additional tsoil' 469 read(*,*,iostat=ierr) val8 if(ierr.ne.0) goto 469 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. & (phisfi(ig).ge.val9) .and. & (phisfi(ig).lt.val10) .and. & (tsurf(ig).ge.val11) .and. & (tsurf(ig).lt.val12) ) then tsurf(ig)=tsurf(ig)*val3 tsurf(ig)=tsurf(ig)+val6 DO l=1,nsoilmx tsoil(ig,l)=tsoil(ig,l)*val7 tsoil(ig,l)=tsoil(ig,l)+val8 ENDDO ENDIF ENDDO c copy latitudes tsurf / tsoil c ----------------------------------------------------- else if (modif(1:len_trim(modif)) .eq. 'copylat') then write(*,*) 'all latitudes : ',rlatu*180./pi write(*,*) 'Choice band to be modified lat1 ?' 465 read(*,*,iostat=ierr) val if(ierr.ne.0) goto 465 write(*,*) 'Choice band copied lat2 ?' 466 read(*,*,iostat=ierr) val2 if(ierr.ne.0) goto 466 write(*,*) 'Choice multiplicative factor' 467 read(*,*,iostat=ierr) val3 if(ierr.ne.0) goto 467 write(*,*) 'Choice additional tsurf' 468 read(*,*,iostat=ierr) val4 if(ierr.ne.0) goto 468 j=1 DO ig=1,ngridmx IF ((latfi(ig)*180./pi.lt.val2+180./iip1) .and. & (latfi(ig)*180./pi.gt.val2-180./iip1)) then randtab(j)=ig j=j+1 ENDIF ENDDO print*,j, ' latitudes found' k=1 DO ig=1,ngridmx IF ((latfi(ig)*180./pi.lt.val+180./iip1) .and. & (latfi(ig)*180./pi.gt.val-180./iip1)) then tsurf(ig)=tsurf(randtab(k))*val3 tsurf(ig)=tsurf(ig)+val4 DO l=1,nsoilmx tsoil(ig,l)=tsoil(randtab(k),l)*val3 tsoil(ig,l)=tsoil(ig,l)+val4 ENDDO k=k+1 ENDIF ENDDO print*,k, ' latitudes copied' IF (j.ne.k) stop c copy longitudes c ----------------------------------------------------- else if (modif(1:len_trim(modif)) .eq. 'copylon') then write(*,*) 'all longitudes : ',rlonu*180./pi write(*,*) 'Choice band to be modified lon1 ?' 475 read(*,*,iostat=ierr) val if(ierr.ne.0) goto 475 write(*,*) 'Choice band copied lon2 ?' 476 read(*,*,iostat=ierr) val2 if(ierr.ne.0) goto 476 j=1 DO ig=2,ngridmx-1 IF ((lonfi(ig)*180./pi.lt.183.) .and. & (lonfi(ig)*180./pi.gt.180.)) then randtab(j)=ig j=j+1 print*, 'lon = ',lonfi(ig) ENDIF ENDDO print*,j, ' longitudes found' k=1 DO ig=2,ngridmx-1 IF ((lonfi(ig)*180./pi.lt.180.) .and. & (lonfi(ig)*180./pi.gt.179.)) then tsurf(ig)=tsurf(randtab(k)) DO l=1,nsoilmx tsoil(ig,l)=tsoil(randtab(k),l) ENDDO qsurf(ig,1)=qsurf(randtab(k),1) phisfi(ig)=phisfi(randtab(k)) k=k+1 ENDIF ENDDO print*,k, ' longitudes copied' IF (j.ne.k) stop CALL gr_fi_dyn(1,ngridmx,iip1,jjp1,phisfi,phis) c copy latlon c ----------------------------------------------------- else if (modif(1:len_trim(modif)) .eq. 'copylatlon') then write(*,*) 'all longitudes : ',rlonu*180./pi write(*,*) 'Choice lat/lon to be copied ?' 495 read(*,*,iostat=ierr) val,val2 if(ierr.ne.0) goto 495 write(*,*) 'Choice indx lon1 lon2 for band ?' 496 read(*,*,iostat=ierr) val3,val4 if(ierr.ne.0) goto 496 write(*,*) 'Choice latitude indx range where we copy ?' 497 read(*,*,iostat=ierr) val5,val6 if(ierr.ne.0) goto 497 ! choice coordinate DO ig=2,ngridmx-1 IF ( (lonfi(ig)*180./pi.gt.val2-1) .and. & (lonfi(ig)*180./pi.lt.val2+1) .and. & (latfi(ig)*180./pi.lt.val+1) .and. & (latfi(ig)*180./pi.gt.val-1) ) then ig0=ig print*,'lat/lon=',latfi(ig)*180./pi,lonfi(ig)*180./pi,ig0 ENDIF ENDDO DO ig=2,ngridmx-1 IF ( (lonfi(ig)*180./pi.lt.val4) .and. & (lonfi(ig)*180./pi.gt.val3) .and. & (latfi(ig)*180./pi.lt.val6) .and. & (latfi(ig)*180./pi.gt.val5) .and. & (qsurf(ig,1).lt.0.001) ) then tsurf(ig)=tsurf(ig0) print*,'corrd=',latfi(ig)*180./pi,lonfi(ig)*180./pi DO l=1,nsoilmx tsoil(ig,l)=tsoil(ig0,l) ENDDO ENDIF ENDDO c Choice surface temperature depending on N2 ice distribution c ----------------------------------------------------- else if (modif(1:len_trim(modif)) .eq. 'tsurfice') then write(*,*) 'Initial soil and surface temp below n2 ?' 107 read(*,*,iostat=ierr) tsurf_bb if(ierr.ne.0) goto 107 write(*,*) 'new surface/soil temp N2 (K):',tsurf_bb write(*,*) 'Initial soil and surface temp when no n2 ?' 108 read(*,*,iostat=ierr) tsurf_bb2 if(ierr.ne.0) goto 108 write(*,*) 'new surface/soil temp when no n2 (K):',tsurf_bb2 DO ig=1,ngridmx if (qsurf(ig,igcm_n2).gt.0.001) then tsurf(ig)=tsurf_bb do l=1,nsoilmx tsoil(ig,l) = tsurf_bb end do else if (tsurf_bb2.ne.0) then tsurf(ig)=tsurf_bb2 do l=1,nsoilmx tsoil(ig,l) = tsurf_bb2 end do endif ENDDO c read an albedo map c ----------------------------------------------------- else if (modif(1:len_trim(modif)) .eq. 'albedomap') then ! Get field 2D fichnom = 'albedo.nc' ierr = NF_OPEN (fichnom, NF_NOWRITE,nid_fi_input) IF (ierr.NE.NF_NOERR) THEN write(6,*)' Problem opening albedo file:',fichnom write(6,*)' ierr = ', ierr CALL ABORT ENDIF ierr = NF_INQ_VARID (nid_fi_input, trim("albedodat"), & nvarid_input) IF (ierr .NE. NF_NOERR) THEN PRINT*, "Could not find asked variable in albedo.nc" CALL abort ENDIF #ifdef NC_DOUBLE ierr = NF_GET_VAR_DOUBLE(nid_fi_input,nvarid_input, & field_input) #else ierr = NF_GET_VAR_REAL(nid_fi_input,nvarid_input,field_input) #endif IF (ierr .NE. NF_NOERR) THEN PRINT*, "Could not get asked variable in albedo.nc" CALL abort ELSE PRINT*, "Got variable in albedo.nc" ENDIF DO ig=1,ngridmx albfi(ig)=field_input(ig) ENDDO c slope : add slope on all longitudes c ----------------------------------------------------- else if (modif(1:len_trim(modif)) .eq. 'slope') then write(*,*) 'choice latitude alt minimum & altitude value (m):' 603 read(*,*,iostat=ierr) val1,val2 if(ierr.ne.0) goto 603 write(*,*) 'choice latitude alt maximum & altitude value (m):' 604 read(*,*,iostat=ierr) val3,val4 if(ierr.ne.0) goto 604 write(*,*) 're scale the pressure :' r = 8.314511E+0 *1000.E+0/mugaz temp=40. write(*,*) 'r, sale height temperature =',r,temp do j=1,jjp1 do i=1,iip1 phisprev= phis(i,j) IF ( ( (rlatu(j)*180./pi.ge.val1) .and. & (rlatu(j)*180./pi.le.val3) ) .or. & ( (rlatu(j)*180./pi.le.val1) .and. & (rlatu(j)*180./pi.ge.val3) )) then val5=abs(val2-val4)/abs(val1-val3)* & abs(val1-rlatu(j)*180./pi)+val2 phis(i,j)=val5*g ps(i,j) = ps(i,j) * . exp((phisprev-phis(i,j))/(temp * r)) ENDIF end do end do c periodicity of surface ps in longitude do j=1,jjp1 ps(1,j) = ps(iip1,j) end do CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,phis,phisfi) c digsp : change topography of SP c ----------------------------------------------------- else if (modif(1:len_trim(modif)) .eq. 'digsp') then write(*,*) 'choice latitudes:' 517 read(*,*,iostat=ierr) val1,val2 if(ierr.ne.0) goto 517 write(*,*) 'choice longitudes:' 518 read(*,*,iostat=ierr) val3,val4 if(ierr.ne.0) goto 518 write(*,*) 'choice phi limite' 519 read(*,*,iostat=ierr) val5 if(ierr.ne.0) goto 519 write(*,*) 'choice change alt (m)' 520 read(*,*,iostat=ierr) val6 if(ierr.ne.0) goto 520 write(*,*) 're scale the pressure :' r = 8.314511E+0 *1000.E+0/mugaz temp=40. write(*,*) 'r, sale height temperature =',r,temp do j=1,jjp1 do i=1,iip1 phisprev= phis(i,j) IF ( ( (rlatu(j)*180./pi.ge.val1) .and. & (rlatu(j)*180./pi.le.val2) ) .and. & ( (rlonv(j)*180./pi.ge.val3) .and. & (rlonv(j)*180./pi.le.val4) ) .and. & (phis(i,j).le.val5)) then phis(i,j)=phis(i,j)+val6*g ps(i,j) = ps(i,j) * . exp((phisprev-phis(i,j))/(temp * r)) ENDIF end do end do c periodicity of surface ps in longitude do j=1,jjp1 ps(1,j) = ps(iip1,j) end do CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,phis,phisfi) c copy field 2D c ----------------------------------------------------- else if (modif(1:len_trim(modif)) .eq. 'copytsoil') then ! Get field 2D fichnom = 'startfi_input.nc' ierr = NF_OPEN (fichnom, NF_NOWRITE,nid_fi_input) IF (ierr.NE.NF_NOERR) THEN write(6,*)' Problem opening file:',fichnom write(6,*)' ierr = ', ierr CALL ABORT ENDIF ! Choice variable to be copied c write(*,*) 'Choice variable to be copied ?' c515 read(*,fmt='(a20)',iostat=ierr) modif c if(ierr.ne.0) goto 515 c write(*,*) 'variable asked :',trim(modif) ierr = NF_INQ_VARID (nid_fi_input, trim("tsurf"),nvarid_input) IF (ierr .NE. NF_NOERR) THEN PRINT*, "Could not find asked variable in startfi_input.nc" CALL abort ENDIF ierr = NF_INQ_VARID (nid_fi_input, trim("tsoil"), & nvarid_inputs) IF (ierr .NE. NF_NOERR) THEN PRINT*, "Could not find asked variable s in startfi_input.nc" CALL abort ENDIF #ifdef NC_DOUBLE ierr = NF_GET_VAR_DOUBLE(nid_fi_input,nvarid_input, & field_input) #else ierr = NF_GET_VAR_REAL(nid_fi_input,nvarid_input,field_input) #endif IF (ierr .NE. NF_NOERR) THEN PRINT*, "Could not get asked variable in startfi_input.nc" CALL abort ELSE PRINT*, "Got variable in startfi_input.nc" ENDIF #ifdef NC_DOUBLE ierr=NF_GET_VAR_DOUBLE(nid_fi_input,nvarid_inputs, & field_inputs) #else ierr=NF_GET_VAR_REAL(nid_fi_input,nvarid_inputs,field_inputs) #endif IF (ierr .NE. NF_NOERR) THEN PRINT*, "Could not get asked variable s in startfi_input.nc" CALL abort ELSE PRINT*, "Got variable s in startfi_input.nc" ENDIF write(*,*) 'Choice domain lon1 lon2 ?' 525 read(*,*,iostat=ierr) val,val2 if(ierr.ne.0) goto 525 write(*,*) 'Choice domain lat1 lat2 ?' 526 read(*,*,iostat=ierr) val3,val4 if(ierr.ne.0) goto 526 write(*,*) 'No change if N2 ice (0) / Change (1) ?' 527 read(*,*,iostat=ierr) val5 if(ierr.ne.0) goto 527 ! Loop DO ig=1,ngridmx IF ( (lonfi(ig)*180./pi.ge.val) .and. & (lonfi(ig)*180./pi.le.val2) .and. & (latfi(ig)*180./pi.ge.val3) .and. & (latfi(ig)*180./pi.le.val4) ) then if (qsurf(ig,igcm_n2).lt.0.001.or.val5.eq.1) then tsurf(ig)=field_input(ig) do l=1,nsoilmx !tsoil(ig,l) = field_input(ig) tsoil(ig,l) = field_inputs(ig,l) end do endif ENDIF ENDDO c modify ice distribution depending on albedo c ----------------------------------------------------- else if (modif(1:len_trim(modif)) .eq. 'mod_distrib_ice') then write(*,*) 'Choice range albedo for CH4 ice' 220 read(*,*,iostat=ierr) val,val2 if(ierr.ne.0) goto 220 write(*,*) 'Choice range albedo for N2 ice' 221 read(*,*,iostat=ierr) val3,val4 if(ierr.ne.0) goto 221 write(*,*) 'Choice extra mass of CH4 ice' 222 read(*,*,iostat=ierr) val5 if(ierr.ne.0) goto 222 write(*,*) 'Choice extra mass of N2 ice' 223 read(*,*,iostat=ierr) val6 if(ierr.ne.0) goto 223 DO ig=1,ngridmx IF ( (albfi(ig).ge.val) .and. & (albfi(ig).le.val2) ) then qsurf(ig,igcm_ch4_ice)=qsurf(ig,igcm_ch4_ice)+val5 ENDIF IF ( (albfi(ig).ge.val3) .and. & (albfi(ig).le.val4) ) then qsurf(ig,igcm_n2)=qsurf(ig,igcm_n2)+val6 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. else if (map_pluto_dat(i,j).eq.4) then CH4_pluto_dat(i,j)=100000. else if (map_pluto_dat(i,j).eq.0) then alb_dat(i,j)=0.3 else if (map_pluto_dat(i,j).eq.6) then alb_dat(i,j)=0.6 else if (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 radp radp=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 = radp * c return end