C====================================================================== PROGRAM newstart c======================================================================= c c c Auteur: Christophe Hourdin/Francois Forget/Yann Wanherdrick c ------ c c Objet: Create or modify the initial state for the LMD Mars GCM c ----- (fichiers NetCDF start et startfi) c c c======================================================================= ! to use 'getin' USE ioipsl_getincom USE planet_h USE comsoil_h use datafile_mod implicit none #include "dimensions.h" #include "dimphys.h" #include "surfdat.h" !#include "comsoil.h" #include "paramet.h" #include "comconst.h" #include "comvert.h" #include "comgeom2.h" #include "control.h" #include "logic.h" #include "description.h" #include "ener.h" #include "temps.h" #include "lmdstd.h" #include "comdissnew.h" #include "clesph0.h" #include "serre.h" #include "netcdf.inc" #include "advtrac.h" #include "tracer.h" #include "callkeys.h" c======================================================================= c Declarations c======================================================================= c Variables dimension du fichier "start_archive" c------------------------------------ CHARACTER relief*3 INTEGER tapeout ! unit numbers for (standard) outputs parameter (tapeout=6) integer tapeerr ! unit number for error message parameter (tapeerr=0) c Variables pour les lectures NetCDF des fichiers "start_archive" c-------------------------------------------------- INTEGER nid_dyn, nid_fi,nid,nvarid, nid_fi_input INTEGER nvarid_input,nvarid_inputs 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 q(iip1,jjp1,llm,nqmx) ! 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 dv(ip1jm,llm),du(ip1jmp1,llm) ! REAL dh(ip1jmp1,llm),dp(ip1jmp1) 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 tsoil(ngridmx,nsoilmx) ! soil temperature REAL emis(ngridmx) ! surface emissivity REAL qsurf(ngridmx,nqmx) REAL qsurf_tmp(ngridmx,nqmx) REAL q2(ngridmx,nlayermx+1) real alb(iip1,jjp1),albfi(ngridmx) ! albedos real ith(iip1,jjp1,nsoilmx),ithfi(ngridmx,nsoilmx) ! 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 field_inputs(ngridmx,nsoilmx) INTEGER i,j,l,isoil,ig,idum,k integer ierr c Variables on the new grid along scalar points c------------------------------------------------------ ! REAL p(iip1,jjp1) REAL t(iip1,jjp1,llm) REAL tset(iip1,jjp1,llm) real phisold_newgrid(iip1,jjp1) real phisold(iip1,jjp1) ! REAL tanh(ip1jmp1) 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) REAL :: beta(iip1,jjp1,llm) ! REAL dteta(ip1jmp1,llm) 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,choice integer randtab(ngridmx) character*80 fichnom integer Lmodif,iq character modif*20 real z_reel(iip1,jjp1) real tempsoil(22),levsoil(22) ! real z_reel(ngridmx) real ith_bb,Tiso,tsurf_bb,tsurf_bb2,geop real ptoto,pcap,patm,airetot,ptotn,patmn character(len=50) :: surfacefile ! "surface.nc" (or equivalent file) ! real ssum 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 INTEGER :: nq,numvanle character(len=20) :: txt ! to store some text ! MONS data: ! real :: MONS_Hdn(iip1,jjp1) ! Hdn: %WEH=Mass fraction of H2O ! real :: MONS_d21(iip1,jjp1) ! ice table "depth" (in kg/m2) ! ! coefficient to apply to convert d21 to 'true' depth (m) ! real :: MONS_coeff ! real :: MONS_coeffS ! coeff for southern hemisphere ! real :: MONS_coeffN ! coeff for northern hemisphere !! real,parameter :: icedepthmin=1.e-3 ! Ice begins at most at that depth c Special Pluto Map from Lellouch & Grundy c------------------------------------------ c integer,parameter :: im_plu=180 ! TB15 coeur c integer,parameter :: jm_plu=89 integer,parameter :: im_plu=360 ! TB15 from topo integer,parameter :: jm_plu=180 c integer,parameter :: im_plu=60 !60 TB15 Lell c integer,parameter :: jm_plu=30 !30 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 c--------------------------------------- ! INTEGER :: visuid ! real :: time_step,t_ops,t_wrt ! CHARACTER*80 :: visu_file preff = 2. ! Pluto ! 610. ! for Mars, instead of 101325. (Earth) pa = 0.2 ! Pluto ! 20 ! for Mars, instead of 500 (Earth) pi=2.*asin(1.) 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 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 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: in nc: 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 if (preff.eq.0) then preff=2. ! Pluto pa=0.2 ! Pluto endif 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?" read(*,*) preffnew write(*,*)"New value of reference pressure pa for purely" write(*,*)"pressure levels (for hybrid coordinates)?" read(*,*) panew preff=preffnew pa=panew 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 call tabfi (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 call tabfi (nid_fi,Lmodif,tab0,day_ini,lllm,p_rad, . p_omeg,p_g,p_cpp,p_mugaz,p_daysec,time) endif rad = p_rad omeg = p_omeg g = p_g cpp = p_cpp mugaz = p_mugaz daysec = p_daysec kappa= 8.314511E+0 *1000.E+0/(mugaz*cpp) c======================================================================= c INITIALISATIONS DIVERSES c======================================================================= ! Load tracer names: call iniadvtrac(nq,numvanle) ! tnom(:) now contains tracer names ! Initialize global tracer indexes (stored in tracer.h) call initracer() ! Load parameters from run.def file CALL defrun_new( 99, .TRUE. ) CALL iniconst CALL inigeom idum=-1 idum=0 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) c======================================================================= c lecture topographie, albedo, inertie thermique, relief sous-maille c======================================================================= if (choix_1.ne.1) then ! pour ne pas avoir besoin du fichier ! surface.dat dans le cas des start ! get datadir OPEN(99,file='callphys.def',status='old',form='formatted' . ,iostat=ierr) CLOSE(99) ! first call to getin will open the file c IF(ierr.NE.0) THEN ! if file callphys.def is found WRITE(tapeout,*) "Problem in Newstart: where is callphys.def?" WRITE(tapeout,*) "" STOP ENDIF call getin("datadir",datadir) write(*,*) '' write(*,*) 'Please choose:' write(*,*) '1- use of surface file netcdf ' write(*,*) '2- define topography manually' write(*,*) '3- skip this part (warning : flat topo will be set)' read(*,*) choice if (choice.eq.1) then write(*,*) 'NC file' CALL datareadnc(relief,phis,alb,surfith, & zmeaS,zstdS,zsigS,zgamS,ztheS) CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,phis,phisfi) else if (choice.eq.2) then ! When not using a "surface.nc" file : define topo manually 302 write(*,*) 'Do you want a flat terrain? (y/n)' read(*,fmt='(a)') yes if (yes.eq.'y') then phis(:,:)=0 else !! Creation topography : Choice longitudes val1 and val2 write(*,*)'Choice range longitude (in deg.), and topography' ierr=1 do while (ierr.ne.0) read(*,*,iostat=ierr) val1,val2 write(*,*)'Values are:',val1,val2 if (ierr.eq.0) then ! got a value ! do a sanity check if((val1.lt.-180.).or.(val1.gt.180.)) then write(*,*)'Longitude should be between -180 and 180 deg!' ierr=1 else if((val2.lt.-180.).or.(val2.gt.180.)) then write(*,*)'Longitude should be between -180 and 180 deg!' ierr=1 else if(val1.ge.val2) then write(*,*)'Longitude 1 should be lower than longitude 2 !' ierr=1 else ! find corresponding iref (nearest longitude) ! note: rlonv(:) contains increasing values of longitude ! starting from -PI to PI do i=1,iip1 ! get val1 write(*,*) 'val1:',rlonv(i)*180./pi,rlonv(i+1)*180./pi if ((rlonv(i)*180./pi.le.val1).and. & (rlonv(i+1)*180./pi.ge.val1)) then ! find which grid point is nearest to val: if (abs(rlonv(i)*180./pi-val1).le. & abs((rlonv(i+1)*180./pi-val1))) then iref1=i else iref1=i+1 endif write(*,*)'Will use nearest grid longitude 1:', & rlonv(iref1)*180./pi endif enddo ! of do j=1,jjp1 do i=1,iip1 ! get val2 if ((rlonv(i)*180./pi.le.val2).and. & (rlonv(i+1)*180./pi.ge.val2)) then ! find which grid point is nearest to val: if (abs(rlonv(i)*180./pi-val2).le. & abs((rlonv(i+1)*180./pi-val2))) then iref2=i else iref2=i+1 endif write(*,*)'Will use nearest grid longitude 2:', & rlonv(iref2)*180./pi endif enddo ! of do j=1,jjp1 endif ! of if((val.lt.0.).or.(val.gt.90)) endif ! of if((val.lt.0.).or.(val.gt.90)) enddo ! of do while !! Choice latitudes val1 and val2 write(*,*)'Choice range latitudes (in deg.), and topography' ierr=1 do while (ierr.ne.0) read(*,*,iostat=ierr) val1,val2 write(*,*)'Values are:',val1,val2 if (ierr.eq.0) then ! got a value ! do a sanity check if((val1.lt.-90.).or.(val1.gt.90.)) then write(*,*)'Latitude should be between -90 and 90 deg!' ierr=1 else if((val2.lt.-90.).or.(val2.gt.90.)) then write(*,*)'Latitude should be between -90 and 90 deg!' ierr=1 else if(val1.ge.val2) then write(*,*)'Latitude 1 should be higuer than latitude 2 !' 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 ! get val1 write(*,*) 'val1:',rlatu(j)*180./pi,rlatu(j+1)*180./pi if ((rlatu(j)*180./pi.ge.val1).and. & (rlatu(j+1)*180./pi.le.val1)) then ! find which grid point is nearest to val: if (abs(rlatu(j)*180./pi-val1).le. & abs((rlatu(j+1)*180./pi-val1))) then jref1=j else jref1=j+1 endif write(*,*)'Will use nearest grid latitude:', & rlatu(jref1)*180./pi endif enddo ! of do j=1,jjp1 do j=1,jjp1 ! get val2 write(*,*) 'val2:',rlatu(j)*180./pi,rlatu(j+1)*180./pi if ((rlatu(j)*180./pi.ge.val2).and. & (rlatu(j+1)*180./pi.le.val2)) then ! find which grid point is nearest to val: if (abs(rlatu(j)*180./pi-val2).le. & abs((rlatu(j+1)*180./pi-val2))) then jref2=j else jref2=j+1 endif write(*,*)'Will use nearest grid latitude:', & rlatu(jref2)*180./pi endif enddo ! of do j=1,jjp1 endif ! of if((val.lt.0.).or.(val.gt.90)) endif ! enddo ! of do while write(*,*) 'Choice of topography (m) (eg: 1000 or -1000) ?' 301 read(*,*,iostat=ierr) val if(ierr.ne.0) goto 301 write(*,*) 'gravity g is : ',g geop=g*val c phis(:,:)=0 phis(iref1:iref2,jref2:jref1)=geop write(*,*) 'phis=',phis write(*,*) 'iip1=',iip1,'jjp1=',jjp1 write(*,*) 'iref1=',iref1,'iref2=',iref2 write(*,*) 'jref1=',jref1,'jref2=',jref2 write(*,*) 'Do you want another topography choice ?' read(*,fmt='(a)') yes if (yes.eq.'y') goto 302 endif ! 'yes flat terrain' zmeaS(:,:)=0 zstdS(:,:)=0 zsigS(:,:)=0 zgamS(:,:)=0 ztheS(:,:)=0 yes=' ' do while ((yes.ne.'y').and.(yes.ne.'n')) write(*,*)'Do you want to change soil albedo and IT (y/n)?' read(*,fmt='(a)') yes end do if (yes.eq.'y') then 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) CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,surfith,surfithfi) CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,alb,albfi) endif !yes for alb and IT changes CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,phis,phisfi) else write(*,*)'OK : skipping topography change' endif ! of if (choice) 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.ne.1) 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(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(:,:) phis(:,:)=phisold_newgrid(:,:) CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,phis,phisfi) ierr= NF_CLOSE(nid) else if (choix_1.eq.1) then ! c'est l'appel a tabfi de phyeta0 qui ! permet de changer les valeurs du ! tab_cntrl Lmodif=1 tab0=0 Lmodif=1 ! Lmodif set to 1 to allow modifications in phyeta0 write(*,*) 'Reading file START' fichnom = 'start.nc' CALL dynetat0(fichnom,nqmx,vcov,ucov,teta,q,masse, . ps,phis,time) write(*,*) 'Reading file STARTFI' fichnom = 'startfi.nc' CALL phyetat0 (fichnom,tab0,Lmodif,nsoilmx,nqmx, . day_ini,time, . tsurf,tsoil,emis,q2,qsurf) ! copy albedo and soil thermal inertia 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 neede later on if reinitializing soil thermal inertia surfithfi(i)=ithfi(i,1) enddo else CALL exit(1) endif dtvr = daysec/FLOAT(day_step) dtphys = dtvr * FLOAT(iphysiq) !preff=preffnew !pa=panew 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(*,*) 'qname : change tracer name' write(*,*) 'q=0 : ALL tracer = zero' write(*,*) 'q=x : give a specific uniform value to one tracer' write(*,*) 'q=kx : initial value of tracer x coeff k' write(*,*) 'qnogcm: initialise 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' c write(*,*) 'emis : change surface emissivity' write(*,*) 'therm_ini_s: Set soil thermal inertia to reference sur &face values' write(*,*) 'inert3d: give uniform value of thermal inertia' write(*,*) 'subsoilice_n: Put deep underground ice layer in northe &rn hemisphere' write(*,*) 'subsoilice_s: Put deep underground ice layer in southe &rn hemisphere' c write(*,*) 'surfalb: change bare groud albedo in startfi' write(*,*) 'reservoir: increase/decrease reservoir of ice' write(*,*) 'reservoir_SP: increase/decrease reservoir ice 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 latitude/longitude' write(*,*) 'copylon: copy tsurf tsoil latitude, n2surf and phisfi' write(*,*) 'tsurfice: surface and soil temperature depending &on if N2 ice is present or not' write(*,*) 'icarus: special option to change surface/soil &temperatures over a latitudinal band' 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 (in SP or outside)' 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 with N2 ice 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 longitude (specific triton)' write(*,*) 'digsp: specific to dig SP' write(*,*) 'bugn2:correct bug of warm n2 due to high resolution' write(*,*) 'bugch4:correct bug of warm ch4 due to high resolution' 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 mass of N2 ice' 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 startfi_input.nc' write(*,*) 'albedomap: read in an albedomap albedo.nc' write(*,*) 'mod_distrib_ice: modify ice distribution from albedo' write(*,*) write(*,*) 'Please note that emis and albedo are set in physiq' write(*,*) write(*,*) 'Change to perform ?' write(*,*) ' (enter keyword or return to end)' write(*,*) read(*,fmt='(a20)') modif if (modif(1:1) .eq. ' ') then write(*,*) 'weird space, exiting' exit ! exit loop on changes endif write(*,*) write(*,*) trim(modif) , ' : ' c TB16 slope : add slope on all longitudes c ----------------------------------------------------- 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 TB17 digsp : change topography of SP c ----------------------------------------------------- elseif (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 TB16 : 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 TB16 topo_sp : change topo of SP c ----------------------------------------------------- elseif (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 TB16 fill_sp inv: fill inverted SP with N2 ice and adjust phisfi (flat SP) c ----------------------------------------------------- elseif (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 !write(*,*) 'TB16 : phisfi',phisfi write(*,*) 'TB16 : qsurf',qsurf(:,igcm_n2) c update new phis and ps CALL gr_fi_dyn(1,ngridmx,iip1,jjp1,phisfi,phis) !write(*,*) 'TB16 : phis',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 TB18 adjust phisfi according to the amount of N2 ice c ----------------------------------------------------- elseif (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 TB18 correct phisfi c ----------------------------------------------------- elseif (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 TB18 correct phisfi c ----------------------------------------------------- elseif (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 TB19 No flat topo c ----------------------------------------------------- elseif (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 TB16 fill_sp : fill SP with N2 ice and adjust phisfi (flat SP) c ----------------------------------------------------- elseif (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 !write(*,*) 'TB16 : phisfi',phisfi write(*,*) 'TB16 : qsurf',qsurf(:,igcm_n2) c update new phis and ps CALL gr_fi_dyn(1,ngridmx,iip1,jjp1,phisfi,phis) !write(*,*) 'TB16 : phis',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 TB17 bugch4 correct bug warm ch4 due to change of resolution c ----------------------------------------------------- elseif (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) val4 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 val5=1,val4 IF ((ig-val5.ge.1).and.qsurf(ig,igcm_n2).eq.0..and. & (qsurf(int(ig-val5),igcm_ch4_ice).gt.0.001).and. & (tsurf(ig-val5).gt.val1 & .and.tsurf(ig-val5).lt.val2)) then tsurf(ig)=tsurf(int(ig-val5)) DO l=1,nsoilmx tsoil(ig,l)=tsoil(int(ig-val5),l) ENDDO goto 548 ELSEIF ((ig+val5.le.ngridmx).and. & qsurf(ig,igcm_n2).eq.0..and. & (qsurf(int(ig+val5),igcm_ch4_ice).gt.0.001).and. & (tsurf(ig+val5).gt.val1 & .and.tsurf(ig+val5).lt.val2)) then tsurf(ig)=tsurf(int(ig+val5)) DO l=1,nsoilmx tsoil(ig,l)=tsoil(int(ig+val5),l) ENDDO goto 548 ENDIF enddo 548 write(*,*) 'found point with ch4' ENDIF ENDIF END DO c TB17 bugn2 correct bug warm n2 due to change of resolution c ----------------------------------------------------- elseif (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) val4 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 val5=1,val4 IF ((ig-val5.ge.1).and. & (qsurf(int(ig-val5),igcm_n2).gt.0.001).and. & (tsurf(ig-val5).gt.val1 & .and.tsurf(ig-val5).lt.val2)) then tsurf(ig)=tsurf(int(ig-val5)) DO l=1,nsoilmx tsoil(ig,l)=tsoil(int(ig-val5),l) ENDDO goto 541 ELSEIF ((ig+val5.le.ngridmx).and. & (qsurf(int(ig+val5),igcm_n2).gt.0.001).and. & (tsurf(ig+val5).gt.val1 & .and.tsurf(ig+val5).lt.val2)) then tsurf(ig)=tsurf(int(ig+val5)) DO l=1,nsoilmx tsoil(ig,l)=tsoil(int(ig+val5),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 TB17 flatnosp : flat topo outside a specific terrain (SP) c ----------------------------------------------------- elseif (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 TB18 flatregion : flat topo specific to region c ----------------------------------------------------- elseif (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 TB18 changetopo c ----------------------------------------------------- elseif (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 TB18 corrtopo: corr topo specific bug c ----------------------------------------------------- elseif (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 TB20 asymtopo: c ----------------------------------------------------- elseif (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 TB16 seticeNH : set the ices in the SP crater with NH topo c ----------------------------------------------------- elseif (modif(1:len_trim(modif)) .eq. 'seticeNH') then open(333,file='./tsoil_180_30',form='formatted') do i=1,22 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 -------------- elseif (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 -------------- elseif (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 'flat : no topography ("aquaplanet")' c ------------------------------------- elseif (modif(1:len_trim(modif)) .eq. 'flat') then ! set topo to zero CALL initial0(ip1jmp1,z_reel) CALL multscal(ip1jmp1,z_reel,g,phis) CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,phis,phisfi) write(*,*) 'topography set to zero.' c Choice for 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 write(*,*) 'New value for ps (Pa) ?' 201 read(*,*,iostat=ierr) patm if(ierr.ne.0) goto 201 write(*,*) write(*,*) ' new ps everywhere (Pa) = ', patm write(*,*) do j=1,jjp1 do i=1,iip1 ps(i,j)=patm enddo enddo endif c TB15 : 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 TB16 : 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 TB23 : 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 TB16 : 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 TB20 : 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 TB15 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 TB17 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 TB17 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 TB17 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 TB20 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 TB15 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 TB23 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 TB23 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 TB20 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 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, nqmx 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,nqmx write(*,'(i3,a3,a20)')iq,' : ',trim(tnom(iq)) enddo write(*,'(a35,i3)') & '(enter tracer number; between 1 and ',nqmx write(*,*)' or any other value to quit this option)' read(*,*) iq if ((iq.ge.1).and.(iq.le.nqmx)) then write(*,*)'Change tracer name ',trim(tnom(iq)),' to ?' read(*,*) txt tnom(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.nqmx)) enddo ! of do while (yes.ne.'y') c q=0 : set tracers to zero c ------------------------- else if (modif(1:len_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, nqmx 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, nqmx DO ig=1,ngridmx qsurf(ig,iq)=0. ENDDO ENDDO c q=x : initialise tracer manually c -------------------------------- else if (modif(1:len_trim(modif)).eq.'q=x') then write(*,*) 'Which tracer do you want to modify ?' do iq=1,nqmx write(*,*)iq,' : ',trim(tnom(iq)),' : ',q(1,1,1,iq) enddo write(*,*) '(choose between 1 and ',nqmx,')' read(*,*) iq write(*,*)'mixing ratio of tracer ',trim(tnom(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 yes=' ' do while ((yes.ne.'y').and.(yes.ne.'n')) write(*,*)'Do you want to change surface value (y/n)?' read(*,fmt='(a)') yes end do if (yes.eq.'y') then write(*,*) 'SURFACE value of tracer ',trim(tnom(iq)), & ' ? (kg/m2)' read(*,*) val DO ig=1,ngridmx qsurf(ig,iq)= val ENDDO endif ! Very specific case ! ---------------------------------------- c DO ig=1,ngridmx c if (ig.lt.ngridmx*11/32+1) then c qsurf(ig,iq)=300 c else if (ig.gt.ngridmx*26/32+1) then c qsurf(ig,iq)=val c else c qsurf(ig,iq)=0 c endif c ENDDO c endif c q=kx : initialise tracer manually x coeff k c -------------------------------- else if (modif(1:len_trim(modif)).eq.'q=kx') then write(*,*) 'Which tracer do you want to modify ?' do iq=1,nqmx write(*,*)iq,' : ',trim(tnom(iq)),' : ',q(1,1,1,iq) enddo write(*,*) '(choose between 1 and ',nqmx,')' read(*,*) iq write(*,*)'coefficient for mmr of tracer ',trim(tnom(iq)), & ' ? (kg/kg)' read(*,*) val DO l=1,llm DO j=1,jjp1 DO i=1,iip1 q(i,j,l,iq)=q(i,j,l,iq)*val ENDDO ENDDO ENDDO c q=kx_loc : initialise tracer manually x coeff k depending on location c -------------------------------- else if (modif(1:len_trim(modif)).eq.'q=kx_loc') then write(*,*) 'Which tracer do you want to modify ?' do iq=1,nqmx write(*,*)iq,' : ',trim(tnom(iq)),' : ',q(1,1,1,iq) enddo write(*,*) '(choose between 1 and ',nqmx,')' read(*,*) iq write(*,*)'coefficient for mmr of tracer ',trim(tnom(iq)), & ' ? (kg/kg)' read(*,*) val write(*,*)'Choice range latitudes ' ierr=1 do while (ierr.ne.0) read(*,*,iostat=ierr) val1,val2 write(*,*)'Values are:',val1,val2 enddo write(*,*)'Choice range longitudes ' ierr=1 do while (ierr.ne.0) read(*,*,iostat=ierr) val3,val4 write(*,*)'Values are:',val3,val4 enddo DO l=1,llm DO j=1,jjp1 DO i=1,iip1 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) ) ) then q(i,j,l,iq)=q(i,j,l,iq)*val ENDIF ENDDO ENDDO ENDDO 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,nqmx write(*,*)iq,' : ',trim(tnom(iq)),' : ',q(1,1,1,iq) enddo write(*,*) '(choose between 1 and ',nqmx,')' 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 TB14 surface emissivity. Removed. Check v7.1 and before c ------------------------------------------------- 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 TB16 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 TB16 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 TB16 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 TB17 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) val7 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(val7) DO l=1,nsoilmx tsoil(ig,l)=tsoil(val7,l) ENDDO ENDIF ENDDO c TB17 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 print*, 'TB18 igref=',igref print*, 'TB18 latref=',latfi(igref)*180./pi print*, 'TB18 lonref=',lonfi(igref)*180./pi 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 print*, 'TB18 igref2=',igref2 print*, 'TB18 latref2=',latfi(igref2)*180./pi print*, 'TB18 lonref2=',lonfi(igref2)*180./pi ! 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) print*, 'TB18 ig c1=',ig print*, 'TB18 ig c1 lat=',latfi(ig)*180./pi print*, 'TB18 ig c1 lon=',lonfi(ig)*180./pi print*, 'TB18 d1=',array_dist(i) print*, 'TB18 a1=',array_angle(i) 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) print*, 'TB18 ig c2=',ig print*, 'TB18 ig c2 lat=',latfi(ig)*180./pi print*, 'TB18 ig c2 lon=',lonfi(ig)*180./pi print*, 'TB18 d2=',val6 print*, 'TB18 a2=',angle ! 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 print*, 'TB18 igref3=',igref3 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 TB18 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 TB16 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 TB16 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 TB16 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 TB16 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 TB16 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 TB17 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 TB17 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 ! choice=2546 !else ! val10=modulo((ig-1),96) ! choice=ig+int(96/2)-val10 ! !choice=int(1+96*(val9-1)+val10) ! write(*,*) 'TB17: choice:',choice,val9,val10 !endif !tsurf(ig)=tsurf(choice) !DO l=1,nsoilmx ! tsoil(ig,l)=tsoil(choice,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 TB17 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) val9 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 ! choice=2546 !else !val10=modulo((ig-1),96) !choice=ig+int(96/2)-val10 !choice=int(1+96*(val9-1)+val10) if (val9.ge.1) then choice=val9 tsurf(ig)=tsurf(choice) DO l=1,nsoilmx tsoil(ig,l)=tsoil(choice,l) ENDDO else if (val9.eq.0) then choice=int(ig/96.)*96+84 print*, 'choice=',choice tsurf(ig)=tsurf(choice) DO l=1,nsoilmx tsoil(ig,l)=tsoil(choice,l) ENDDO endif ENDIF ENDDO c TB16 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 TB17 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 TB18 checktsoil : changing tsoil if no n2 c ---------------------------------------------------------- else if (modif(1:len_trim(modif)) .eq. 'checktsoil') then write(*,*) 'selecting area where tsoil to be check ' write(*,*) 'Choice band : lat1 and lat2 ?' 511 read(*,*,iostat=ierr) val,val1 if(ierr.ne.0) goto 511 write(*,*) 'Choice band : lon1 and lon2 ?' 512 read(*,*,iostat=ierr) val2,val3 if(ierr.ne.0) goto 512 write(*,*) 'Choice temp min max in which change is made' 513 read(*,*,iostat=ierr) val4,val5 if(ierr.ne.0) goto 513 write(*,*) 'Choice phisfi min max where change n2' 514 read(*,*,iostat=ierr) val6,val11 if(ierr.ne.0) goto 514 DO ig=1,ngridmx IF ( (latfi(ig)*180./pi.ge.val) .and. & (latfi(ig)*180./pi.le.val1) .and. & (lonfi(ig)*180./pi.ge.val2) .and. & (lonfi(ig)*180./pi.le.val3) .and. & (((tsurf(ig).ge.val4) .and. & (tsurf(ig).le.val5)) .or. & ((tsoil(ig,nsoilmx).ge.val4) .and. & (tsoil(ig,nsoilmx).le.val5))) .and. & (phisfi(ig).ge.val6) .and. & (phisfi(ig).le.val11) .and. & (qsurf(ig,igcm_n2).le.0.001) ) then ! DO i=1,ngridmx ! IF ( (latfi(i).eq.latfi(ig)) .and. ! & (lonfi(i).eq.0.) ) then ! tsurf(ig)=50. !qsurf(ig,igcm_ch4_ice)=qsurf(i,igcm_ch4_ice) ! DO l=1,nsoilmx tsoil(ig,l)=50. !tsoil(i,l) ENDDO !ENDIF !ENDDO ENDIF ENDDO c TB20 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 TB16 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 TB16 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) ! Check: ! do i=1,iip1 ! do j=1,jjp1 ! do isoil=1,nsoilmx ! write(77,*) i,j,isoil," ",ith(i,j,isoil) ! enddo ! enddo ! enddo ! TB14 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 ! endif ! write(*,*)'iref:',iref,' jref:',jref ! write(*,*)'layer',layer ! write(*,*)'mlayer',mlayer ! 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) ! do i=1,nsoilmx ! write(*,*)'i:',i,'ithfi(1,i):',ithfi(1,i) ! enddo ! 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 'surfalb' : TB15 change albedo in start : bare ground c removed. Check v7.1 and before c ---------------------------------------------------------- c ---------------------------------------------------------- c 'reservoir_SP' : TB16 increase or decrase ices reservoir in SP by factor c ---------------------------------------------------------- else if (modif(1:len_trim(modif)).eq.'reservoir_SP') then write(*,*) "Increase/Decrease N2 reservoir by factor :" read(*,*) val7 write(*,*) "Increase/Decrease CH4 reservoir by factor :" read(*,*) val8 write(*,*) "Increase/Decrease CO reservoir by factor :" read(*,*) val9 ! Definition SP: val3=-33. !lat1 val4=50. !lat2 val5=140. ! lon1 val6=-155. ! lon2 choice=-50. ! min phisfi write(*,*) 'def SP :' write(*,*) 'lat :',val3,val4 write(*,*) 'lon :',val5,'180 / -180',val6 write(*,*) 'choice phisfi min :',choice DO ig=1,ngridmx IF ( (latfi(ig)*180./pi.ge.val3) .and. & (latfi(ig)*180./pi.le.val4) .and. & ( ((lonfi(ig)*180./pi.ge.-180.) .and. & (lonfi(ig)*180./pi.le.val6)) .or. & ((lonfi(ig)*180./pi.ge.val5) .and. & (lonfi(ig)*180./pi.le.180.))) ) then IF ((phisfi(ig).le.choice)) then !.and. c & (qsurf(ig,igcm_n2).ge.50)) then qsurf(ig,igcm_n2)=qsurf(ig,igcm_n2)*val7 qsurf(ig,igcm_ch4_ice)=qsurf(ig,igcm_ch4_ice)*val8 qsurf(ig,igcm_co_ice)=qsurf(ig,igcm_co_ice)*val9 ENDIF ENDIF ENDDO c ---------------------------------------------------------- c 'reservoir' : TB15 increase or decrase ices reservoir by factor c ---------------------------------------------------------- else if (modif(1:len_trim(modif)).eq.'reservoir') then write(*,*) "Increase/Decrease N2 reservoir by factor :" read(*,*) val write(*,*) "Increase/Decrease CH4 reservoir by factor :" read(*,*) val2 write(*,*) "Increase/Decrease CO reservoir by factor :" read(*,*) val3 DO ig=1,ngridmx qsurf(ig,igcm_n2)=qsurf(ig,igcm_n2)*val qsurf(ig,igcm_ch4_ice)=qsurf(ig,igcm_ch4_ice)*val2 qsurf(ig,igcm_co_ice)=qsurf(ig,igcm_co_ice)*val3 ENDDO c -------------------------------------------------------- c 'plutomap' : initialize pluto ices distribution c -------------------------------------------------------- else if (modif(1:len_trim(modif)).eq.'plutomap') then write(*,*) 'pluto_map.dat has to be in your simulation file !!' open(27,file='pluto_map.dat') N2_pluto_dat(:,:) =0. CH4_pluto_dat(:,:) =0. CO_pluto_dat(:,:) =0. ! Dimension file pluto_map IF (jm_plu.ne.180) then write(*,*) 'Newstart : you should set jm_plu to 180' stop ENDIF ! Read values do j=1,jm_plu+1 read(27,*) (map_pluto_dat(i,j),i=1,im_plu) do i=1,im_plu !!!!!! BTD_v2 if (map_pluto_dat(i,j).eq.3) then N2_pluto_dat(i,j)=100000. elseif (map_pluto_dat(i,j).eq.4) then CH4_pluto_dat(i,j)=100000. elseif (map_pluto_dat(i,j).eq.0) then alb_dat(i,j)=0.3 elseif (map_pluto_dat(i,j).eq.6) then alb_dat(i,j)=0.6 elseif (map_pluto_dat(i,j).eq.7) then alb_dat(i,j)=0.1 endif !!!!!! BTD !if (map_pluto_dat(i,j).eq.3) then ! CH4_pluto_dat(i,j)=100000. !endif end do end do close (27) ! Interpolate !! latitude and longitude in REindexed pluto map : latv_plu(1)=90. *pi/180. do j=2,jm_plu latv_plu(j)= (90-(j-1 -0.5)*180./(jm_plu-1))*pi/180. end do latv_plu(jm_plu+1)= -90. *pi/180. do i=1,im_plu+1 lonu_plu(i)=(-180+ (i-1)*360./(im_plu) ) *pi/180. enddo !Reindexation to shift the longitude grid like a LMD GCM grid... do j=1,jm_plu+1 N2_pluto_rein(1,j)=(N2_pluto_dat(1,j)+ & N2_pluto_dat(im_plu,j))/2 CH4_pluto_rein(1,j)=(CH4_pluto_dat(1,j)+ & CH4_pluto_dat(im_plu,j))/2 CO_pluto_rein(1,j)=(CO_pluto_dat(1,j)+ & CO_pluto_dat(im_plu,j))/2 alb_rein(1,j)=(alb_dat(1,j)+ & alb_dat(im_plu,j))/2 do i=2,im_plu N2_pluto_rein(i,j)=(N2_pluto_dat(i-1,j)+ & N2_pluto_dat(i,j))/2 CH4_pluto_rein(i,j)=(CH4_pluto_dat(i-1,j)+ & CH4_pluto_dat(i,j))/2 CO_pluto_rein(i,j)=(CO_pluto_dat(i-1,j)+ & CO_pluto_dat(i,j))/2 alb_rein(i,j)=(alb_dat(i-1,j)+ & alb_dat(i,j))/2 end do N2_pluto_rein(im_plu+1,j) = N2_pluto_rein(1,j) CH4_pluto_rein(im_plu+1,j) = CH4_pluto_rein(1,j) CO_pluto_rein(im_plu+1,j) = CO_pluto_rein(1,j) alb_rein(im_plu+1,j) = alb_rein(1,j) end do call interp_horiz(N2_pluto_rein,N2_pluto_gcm, & im_plu,jm_plu,iim,jjm,1,lonu_plu,latv_plu,rlonu,rlatv) call interp_horiz(CH4_pluto_rein,CH4_pluto_gcm, & im_plu,jm_plu,iim,jjm,1,lonu_plu,latv_plu,rlonu,rlatv) call interp_horiz(CO_pluto_rein,CO_pluto_gcm, & im_plu,jm_plu,iim,jjm,1,lonu_plu,latv_plu,rlonu,rlatv) call interp_horiz(alb_rein,alb_gcm, & im_plu,jm_plu,iim,jjm,1,lonu_plu,latv_plu,rlonu,rlatv) ! Switch dyn grid to fi grid qsurf_tmp(:,:) =0. CALL gr_dyn_fi(1,iip1,jjp1,ngridmx, & N2_pluto_gcm,qsurf_tmp(1,igcm_n2)) if (igcm_ch4_ice.ne.0) then CALL gr_dyn_fi(1,iip1,jjp1,ngridmx, & CH4_pluto_gcm,qsurf_tmp(1,igcm_ch4_ice)) endif if (igcm_co_ice.ne.0) then CALL gr_dyn_fi(1,iip1,jjp1,ngridmx, & CO_pluto_gcm,qsurf_tmp(1,igcm_co_ice)) endif alb=alb_gcm CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,alb_gcm,albfi) !print*, 'N2dat=',N2_pluto_gcm !print*, 'COdat=',CO_pluto_gcm !print*, 'CH4dat=',CH4_pluto_gcm print*, 'N2dat=',qsurf_tmp(:,igcm_n2) print*, 'COdat=',qsurf_tmp(:,igcm_co_ice) print*, 'CH4dat=',qsurf_tmp(:,igcm_ch4_ice) ! Fill DO ig=1,ngridmx qsurf(ig,igcm_n2)=qsurf(ig,igcm_n2)+qsurf_tmp(ig,igcm_n2) if (igcm_ch4_ice.ne.0) then qsurf(ig,igcm_ch4_ice)=qsurf(ig,igcm_ch4_ice)+ & qsurf_tmp(ig,igcm_ch4_ice) endif if (igcm_co_ice.ne.0) then qsurf(ig,igcm_co_ice)=qsurf(ig,igcm_co_ice)+ & qsurf_tmp(ig,igcm_co_ice) endif ENDDO ! ***************** end if ! of if (modif(1:len_trim(modif)) .eq. '...') elseif ... 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 = 8.314511E+0 *1000.E+0/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 periodicity of surface ps in longitude do j=1,jjp1 ps(1,j) = ps(iip1,j) end do end if c======================================================================= c======================================================================= 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, beta, 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 if (choix_1.eq.0) then CALL massdair( p3d, masse ) 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,:) ! get option fast or true to skip inidissip if fast model used OPEN(99,file='callphys.def',status='old',form='formatted' . ,iostat=ierr) CLOSE(99) ! first call to getin will open the file c IF(ierr.NE.0) THEN ! if file callphys.def is found WRITE(tapeout,*) "Problem in Newstart: where is callphys.def?" WRITE(tapeout,*) "" STOP ENDIF WRITE(tapeout,*) " " WRITE(tapeout,*) "Option fast (nogcm):" call getin("fast",fast) WRITE(tapeout,*)" fast = ",fast WRITE(tapeout,*) " " IF (.NOT.fast) THEN CALL inidissip ( lstardis, nitergdiv, nitergrot, niterh, * tetagdiv, tetagrot , tetatemp ) ENDIF itau=0 if (choix_1.eq.0) then day_ini=int(date) endif 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 ) write(*,*) 'TB17 : ',preff,pa c CALL caldyn c $ ( itau,ucov,vcov,teta,ps,masse,pk,pkf,phis , c $ phi,conser,du,dv,dteta,dp,w, pbaru,pbarv, day_ini ) CALL dynredem0("restart.nc",day_ini,anneeref,phis,nqmx) CALL dynredem1("restart.nc",0.0,vcov,ucov,teta,q,nqmx,masse,ps) C Ecriture etat initial physique call physdem1("restartfi.nc",lonfi,latfi,nsoilmx,nqmx, . dtphys,float(day_ini), . time,tsurf,tsoil,emis,q2,qsurf, . airefi,albfi,ithfi,zmea,zstd,zsig,zgam,zthe) write(*,*) '***************************************' write(*,*) 'FINAL VALUES IN RESTART: ' write(*,*) 'check consistency' write(*,*) 'to change values in restart use start2archive' write(*,*) '***************************************' write(*,*) 'radius= ',rad write(*,*) 'omeg= ',omeg write(*,*) 'gravity= ',g write(*,*) 'cpp= ', cpp write(*,*) 'mugaz= ',mugaz write(*,*) 'daysec=',daysec write(*,*) 'kappa=', 8.314511E+0 *1000.E+0/(mugaz*cpp) 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//) end !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! function dist_pluto(lat1,lon1,lat2,lon2) implicit none real dist_pluto real lat1,lon1,lat2,lon2 real dlon, dlat real a,c real rad rad=1190. ! km dlon = lon2 - lon1 dlat = lat2 - lat1 a = (sin(dlat/2))**2 + cos(lat1) * cos(lat2) * (sin(dlon/2))**2 c = 2 * atan2( sqrt(a), sqrt(1-a) ) dist_pluto = rad * c return end