Index: trunk/LMDZ.PLUTO/libf/dynphy_lonlat/phypluto/newstart.F =================================================================== --- trunk/LMDZ.PLUTO/libf/dynphy_lonlat/phypluto/newstart.F (revision 3275) +++ trunk/LMDZ.PLUTO/libf/dynphy_lonlat/phypluto/newstart.F (revision 3371) @@ -6,11 +6,10 @@ c Auteur: Christophe Hourdin/Francois Forget/Yann Wanherdrick c ------ -c Derniere modif : 12/03 -c +c Adapted to Pluto: Tanguy Bertrand +c Derniere modif : 06/2024 c c Objet: Create or modify the initial state for the LMD Mars GCM c ----- (fichiers NetCDF start et startfi) c -c c======================================================================= @@ -19,10 +18,9 @@ & is_master use infotrac, only: infotrac_init, nqtot, tname - USE tracer_h, ONLY: igcm_n2 + USE tracer_h, ONLY: igcm_n2,igcm_ch4_ice,igcm_co_ice,igcm_haze, + & igcm_prec_haze,igcm_co_gas,igcm_ch4_gas,noms USE comsoil_h, ONLY: nsoilmx, layer, mlayer, inertiedat USE surfdat_h, ONLY: phisfi, albedodat, & zmea, zstd, zsig, zgam, zthe -! TB24 USE nonoro_gwd_ran_mod, ONLY: du_nonoro_gwd, dv_nonoro_gwd, -! & east_gwstress, west_gwstress use datafile_mod, only: datadir, surfdir use ioipsl_getin_p_mod, only: getin_p @@ -30,6 +28,4 @@ use phyredem, only: physdem0, physdem1 use iostart, only: open_startphy -! use slab_ice_h, only:noceanmx -! USE ocean_slab_mod, ONLY: nslay use filtreg_mod, only: inifilr USE mod_const_mpi, ONLY: COMM_LMDZ @@ -66,5 +62,4 @@ CHARACTER relief*3 - c Variables pour les lectures NetCDF des fichiers "start_archive" c-------------------------------------------------- @@ -90,6 +85,4 @@ 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) @@ -106,26 +99,25 @@ REAL tsurf(ngridmx) ! surface temperature REAL,ALLOCATABLE :: tsoil(:,:) ! soil temperature -! REAL n2ice(ngridmx) ! N2 ice layer REAL emis(ngridmx) ! surface emissivity real emisread ! added by RW REAL,ALLOCATABLE :: qsurf(:,:) + REAL,ALLOCATABLE :: qsurf_tmp(:,:) REAL q2(ngridmx,llm+1) -! REAL rnaturfi(ngridmx) real alb(iip1,jjp1),albfi(ngridmx) ! albedos real,ALLOCATABLE :: ith(:,:,:),ithfi(:,:) ! thermal inertia (3D) real surfith(iip1,jjp1),surfithfi(ngridmx) ! surface thermal inertia (2D) -! REAL latfi(ngridmx),lonfi(ngridmx),airefi(ngridmx) - - INTEGER i,j,l,isoil,ig,idum + REAL :: latfi(ngridmx),lonfi(ngridmx),airefi(ngridmx) + + INTEGER i,j,l,isoil,ig,idum,k real mugaz ! molar mass of the atmosphere - integer ierr,iref + 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 :: teta(iip1, jjp1, llm) REAL :: pk(iip1,jjp1,llm) @@ -135,5 +127,4 @@ REAL :: xpn,xps,xppn(iim),xpps(iim) REAL :: p3d(iip1, jjp1, llm+1) -! REAL dteta(ip1jmp1,llm) c Variable de l'ancienne grille @@ -145,5 +136,5 @@ c variables diverses c------------------- - real choix_1,pp + real choix_1,pp,choice character*80 fichnom character*250 filestring @@ -153,9 +144,11 @@ real tsud,albsud,alb_bb,ith_bb,Tiso,Tabove real ptoto,pcap,patm,airetot,ptotn,patmn,psea -! real ssum character*1 yes logical :: flagtset=.false. , flagps0=.false. - real val, val2, val3, val4 ! to store temporary variables + 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 @@ -174,17 +167,43 @@ real fact2 +c Special Pluto Map from Lellouch & Grundy +c------------------------------------------ + integer,parameter :: im_plu=360 ! from topo + integer,parameter :: jm_plu=180 + real latv_plu(jm_plu+1),lonu_plu(im_plu+1) + real map_pluto_dat(im_plu,jm_plu+1) + + real N2_pluto_dat(im_plu,jm_plu+1) + real CH4_pluto_dat(im_plu,jm_plu+1) + real CO_pluto_dat(im_plu,jm_plu+1) + real alb_dat(im_plu,jm_plu+1) + + real N2_pluto_rein(im_plu+1,jm_plu+1) + real CH4_pluto_rein(im_plu+1,jm_plu+1) + real CO_pluto_rein(im_plu+1,jm_plu+1) + real alb_rein(im_plu+1,jm_plu+1) + real N2_pluto_gcm(iip1,jjp1) + real CH4_pluto_gcm(iip1,jjp1) + real CO_pluto_gcm(iip1,jjp1) + real alb_gcm(iip1,jjp1) + +c Special Topo Map mountain + real lat0, lon0 + integer ig0 + real dist_pluto + real radm,height, phisprev, temp + real preffnew,panew +c Special copy of terrain + real actualmin,angle + integer array_ind(ngridmx) + real array_dist(ngridmx) + real array_angle(ngridmx) c sortie visu pour les champs dynamiques -c--------------------------------------- -! INTEGER :: visuid -! real :: time_step,t_ops,t_wrt -! CHARACTER*80 :: visu_file - !TB24 CALL conf_gcm( 99, .TRUE. ) - cpp = 0. - preff = 0. - pa = 0. ! to ensure disaster rather than minor error if we don`t + preff = 2. + pa = 0.2 ! to ensure disaster rather than minor error if we don`t ! make deliberate choice of these values elsewhere. @@ -204,4 +223,5 @@ allocate(q(iip1,jjp1,llm,nqtot)) allocate(qsurf(ngridmx,nqtot)) + allocate(qsurf_tmp(ngridmx,nqtot)) ! get value of nsoilmx and allocate corresponding arrays @@ -272,5 +292,4 @@ endif - c======================================================================= @@ -354,6 +373,4 @@ . nqtot,day_ini,time, . tsurf,tsoil,emis,q2,qsurf) !) ! temporary modif by RDW -! . cloudfrac,totalfrac,hice,rnat,pctsrf_sic,tslab,tsea_ice, -! . sea_ice) ! copy albedo and soil thermal inertia on (local) physics grid @@ -432,4 +449,21 @@ endif +c Initialisation coordonnees /aires +c ------------------------------- +! Note: rlatu(:) and rlonv(:) are commons defined in "comgeom.h" +! rlatu() and rlonv() are given in radians + latfi(1)=rlatu(1) + lonfi(1)=0. + DO j=2,jjm + DO i=1,iim + latfi((j-2)*iim+1+i)=rlatu(j) + lonfi((j-2)*iim+1+i)=rlonv(i) + ENDDO + ENDDO + latfi(ngridmx)=rlatu(jjp1) + lonfi(ngridmx)=0. + + CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,aire,airefi) + rad = p_rad omeg = p_omeg @@ -438,10 +472,4 @@ mugaz = p_mugaz daysec = p_daysec - - if (p_omeg .eq. -9999.) then - p_omeg = 8.*atan(1.)/p_daysec - print*,"new value of omega",p_omeg - endif - kappa = 8.314*1000./(p_mugaz*p_cpp) ! added by RDW @@ -536,5 +564,4 @@ endif ! of if (choix_1.eq.0) - c======================================================================= c Lecture des fichiers (start ou start_archive) @@ -549,6 +576,4 @@ & surfith,nid) -! & rnat,pctsrf_sic,tslab,tsea_ice,sea_ice) -! TB24 & du_nonoro_gwd,dv_nonoro_gwd,east_gwstress,west_gwstress) write(*,*) "OK, read start_archive file" ! copy soil thermal inertia @@ -572,19 +597,86 @@ do ! infinite loop on list of changes - write(*,*) - write(*,*) - write(*,*) 'List of possible changes :' - write(*,*) '~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~' - write(*,*) - write(*,*) 'flat : no topography ("aquaplanet")' - write(*,*) 'set_ps_to_preff : used if changing preff with topo' - write(*,*) 'qname : change tracer name' - write(*,*) 't=profile : read temperature profile in profile.in' - write(*,*) 'q=0 : ALL tracer =zero' - write(*,*) 'q=x : give a specific uniform value to one tracer' - write(*,*) 'qs=x : give a uniform value to a surface tracer' - write(*,*) 'q=profile : specify a profile for a tracer' - write(*,*) 'subsoil_all : set seasonal subsurface thermal inertia' - write(*,*) 'diurnal_TI : set diurnal subsurface thermal inertia' + write(*,*) + write(*,*) + write(*,*) 'List of possible changes :' + write(*,*) '~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~' + write(*,*) + write(*,*) 'flat : no topography ("aquaplanet")' + write(*,*) 'set_ps_to_preff : used if changing preff with topo' + write(*,*) 'qname : change tracer name' + write(*,*) 't=profile : read temperature profile in profile.in' + write(*,*) 'q=0 : ALL tracer =zero' + write(*,*) 'q=x : give a specific uniform value to one tracer' + write(*,*) 'qs=x : give a uniform value to a surface tracer' + write(*,*) 'q=profile : specify a profile for a tracer' + write(*,*) 'qnogcm: initial tracer for GCM from nogcm (CH4,CO)' + write(*,*) 'isotherm : Isothermal Temperatures' + write(*,*) 'tempprof : specific Temperature profile' + write(*,*) 'initsoil : initialize soil Temperatures' + write(*,*) 'ptot : change total pressure' + write(*,*) 'therm_ini_s: set soil TI to reference surf. values' + write(*,*) 'inert3d: give uniform value of thermal inertia' + write(*,*) 'subsoilice_n: put deep underground ice in N. hemis' + write(*,*) 'subsoilice_s: put deep underground ice in S. hemis' + write(*,*) 'reservoir: increase/decrease reservoir of ice' + write(*,*) 'reservoir_SP: increase/decrease reservoir in SP' + write(*,*) 'plutomap: initialize surface like pluto from ...' + write(*,*) 'subsoil_n2: diu-sea thermal inertia for N2 only' + write(*,*) 'subsoil_ch4: diu-sea thermal inertia for CH4 only' + write(*,*) 'subsoil_all: diu-sea thermal inertia for all terr' + write(*,*) 'subsoil_alb: diu-sea thermal inertia from albedo' + write(*,*) 'diurnal_TI: diurnal thermal inertia for all terr' + write(*,*) 'initsurf: initial surface temperature (global)' + write(*,*) 'modtsurf: initial surface temperature (global)' + write(*,*) 'copylat: copy tsurf and tsoil latitude' + write(*,*) 'copylatlon: copy tsurf and tsoil lat/lon' + write(*,*) 'copylon: copy tsurf tsoil lat, n2surf and phisfi' + write(*,*) 'tsurfice: temperature depending on N2 ice' + write(*,*) 'icarus: option to change surface/soil temperature' + write(*,*) 'icarus_ch4: special option to add ch4' + write(*,*) 'delfrostch4: delete ch4 frost' + write(*,*) 'modch4: remove/modify amount ch4 frost' + write(*,*) 'modch4n2: modify amount ch4 frost according N2' + write(*,*) 'modco: remove/modify amount co frost' + write(*,*) 'modn2: remove/modify amount n2 ice' + write(*,*) 'modcoch4: remove/modify co ch4 where no n2 ' + write(*,*) 'checktsoil: change tsoil where no n2 ' + write(*,*) 'non2noco: if no n2, no co ' + write(*,*) 'modn2_topo: modify n2 ice topo according to topo' + write(*,*) 'modwheren2: modify n2 where already n2' + write(*,*) 'modn2HD: modify n2 for HD runs' + write(*,*) 'modn2HD_SP: modify n2 for HD runs in SP' + write(*,*) 'globch4co: add/remove global amount ch4co frost' + write(*,*) 'source_ch4: add source ch4' + write(*,*) 'nomountain: remove pluto mountains (!)' + write(*,*) 'relief: add pluto crater or mountain' + write(*,*) 'seticeNH: set ice for initial start with NH topo' + write(*,*) 'topo_sp: change topography of SP' + write(*,*) 'fill_sp: fill sp with N2 ice and adjust phisfi' + write(*,*) 'fill_sp_inv: fill inverted sp and adjust' + write(*,*) 'subsoil_SP: diu-sea thermal inertia for SP ' + write(*,*) 'bladed: add ch4 on bladed terrains' + write(*,*) 'cpbladed: add extension bladed terrains' + write(*,*) 'slope: add slope over all long. (for triton)' + write(*,*) 'digsp: specific to dig SP' + write(*,*) 'bugn2: correct bug of warm n2 due to HighRes' + write(*,*) 'bugch4: correct bug of warm ch4 due to HighRes' + write(*,*) 'flatnosp : topo flat except SP (topo controled)' + write(*,*) 'flatregion: topo flat for specific region' + write(*,*) 'changetopo: change topo' + write(*,*) 'asymtopo: N-S asym topo in tanh' + write(*,*) 'corrtopo: correction topo bug' + write(*,*) 'adjustphi: adjust phisfi according to N2 mass' + write(*,*) 'correctphi: adjust phisfi' + write(*,*) 'correctps: test to correct ps' + write(*,*) 'toponoise: no flat topo' + write(*,*) 'asymtriton: asymetry in longitude for triton' + write(*,*) 'copytsoil: copy 2D tsoil field from nc file' + write(*,*) 'albedomap: read in an albedomap albedo.nc' + write(*,*) 'mod_distrib_ice: modify ice distrib. from albedo' + + write(*,*) + write(*,*) 'Please note that emis and albedo are set in physiq' + write(*,*) write(*,*) @@ -825,79 +917,976 @@ endif - - -c subsoil_all : initialize subsurface thermal inertia -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(*,*)'(here, 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 in subsoil_all:',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) +c qnogcm : initialise tracer from nogcm (CH4, CO) +c -------------------------------- + else if (modif(1:len_trim(modif)).eq.'qnogcm') then + write(*,*) 'Which tracer do you want to modify ?' + write(*,*) 'Should be ch4_gas and co_gas' + do iq=1,nqtot + write(*,*)iq,' : ',trim(noms(iq)),' : ',q(1,1,1,iq) + enddo + write(*,*) '(choose between 1 and ',nqtot,')' + read(*,*) iq + DO l=1,llm + DO j=1,jjp1 + DO i=1,iip1 + q(i,j,l,iq)=q(i,j,1,iq) + ENDDO + ENDDO + ENDDO + +c isothermal temperatures +c ------------------------------------------------ + else if (modif(1:len_trim(modif)) .eq. 'isotherm') then + + write(*,*)'Isothermal temperature of the atmosphere' + write(*,*) 'Value of THIS temperature ? :' + 203 read(*,*,iostat=ierr) Tset(1,1,1) + if(ierr.ne.0) goto 203 + flagtset=.true. + DO l=1,llm + DO j=1,jjp1 + DO i=1,iip1 + Tset(i,j,l)=Tset(1,1,1) + ENDDO + ENDDO + ENDDO + write(*,*) 'atmospheric temp =', Tset(2,2,2) + +c specific temperature profile +c ------------------------------------------------ + else if (modif(1:len_trim(modif)) .eq. 'tempprof') then + + write(*,*) 'phi=' + write(*,*) phi(1,1,:) + write(*,*) 'temperature profile in the atmosphere' + write(*,*) 'Value of THIS temperature ? :' + 204 read(*,*,iostat=ierr) Tset(1,1,1) + if(ierr.ne.0) goto 204 + flagtset=.true. + DO l=1,llm + DO j=1,jjp1 + DO i=1,iip1 + Tset(i,j,l)=Tset(1,1,1) + ENDDO + ENDDO + ENDDO + write(*,*) 'atmospheric temp =', Tset(2,2,2) + +c initsoil: subsurface temperature +c --------------------------------- + else if (modif(1:len_trim(modif)) .eq. 'initsoil') then + + write(*,*)'Isothermal temperature of the subsurface' + write(*,*) 'Value of this temperature ? :' + 303 read(*,*,iostat=ierr) Tiso + if(ierr.ne.0) goto 303 + + do l=1,nsoilmx + do ig=1, ngridmx + tsoil(ig,l) = Tiso + end do + end do + +c icarus : changing tsoil tsurf on latitudinal bands +c ----------------------------------------------------- + else if (modif(1:len_trim(modif)) .eq. 'icarus') then + + write(*,*) 'At which lat lon do you want to extract the + & reference tsurf/tsoil profile ?' + 407 read(*,*,iostat=ierr) val,val2 + if(ierr.ne.0) goto 407 + write(*,*) 'You want lat =',val + write(*,*) 'You want lon =',val2 + dist=0 + ref=1000 + do j=1,ngridmx-1 + dist=sqrt((lonfi(j)*180./pi-val2)**2+ + & (latfi(j)*180./pi-val)**2) + IF (dist.lt.ref) then + ref=dist + jref=j + ENDIF + ENDDO + + write(*,*)'Will use nearest grid point which is:', + & latfi(jref)*180./pi,lonfi(jref)*180./pi + + ! Extraction of the profile + write(*,*) 'Profile is =',tsoil(jref,:) + write(*,*) 'tsurf is =',tsurf(jref) + write(*,*) 'Choice lat for latitudinal band with this profile' + 408 read(*,*,iostat=ierr) val3 + if(ierr.ne.0) goto 408 + write(*,*) 'Choice delta lat for this latitudinal band' + 409 read(*,*,iostat=ierr) val4 + if(ierr.ne.0) goto 409 + write(*,*) 'Choice delta temp (K) for this latitudinal band' + 410 read(*,*,iostat=ierr) val5 + if(ierr.ne.0) goto 410 + write(*,*) 'How much N2 ice should I put on this band (kg/m2)' + 415 read(*,*,iostat=ierr) choice + if(ierr.ne.0) goto 415 + write(*,*) 'Choice lat2' + 411 read(*,*,iostat=ierr) val6 + if(ierr.ne.0) goto 411 + write(*,*) 'Choice delta lat for this latitudinal band' + 412 read(*,*,iostat=ierr) val7 + if(ierr.ne.0) goto 412 + write(*,*) 'Choice delta temp (K) for this latitudinal band' + 413 read(*,*,iostat=ierr) val8 + if(ierr.ne.0) goto 413 + + DO ig=1,ngridmx + IF ( (latfi(ig)*180./pi.ge.(val3-val4/2.)) .and. + & (latfi(ig)*180./pi.le.(val3+val4/2.)) .and. + & (qsurf(ig,igcm_n2).lt.0.001) ) then + tsurf(ig)=tsurf(jref)+val5 + DO l=1,nsoilmx + tsoil(ig,l)=tsoil(jref,l)+val5 + ENDDO + qsurf(ig,igcm_n2)=choice + ENDIF + IF ( (latfi(ig)*180./pi.ge.(val6-val7/2.)) .and. + & (latfi(ig)*180./pi.le.(val6+val7/2.)) .and. + & (qsurf(ig,igcm_n2).lt.0.001) ) then + tsurf(ig)=tsurf(jref)+val8 + DO l=1,nsoilmx + tsoil(ig,l)=tsoil(jref,l)+val8 + ENDDO + ENDIF + ENDDO + +c icarus_ch4 : changing tsoil tsurf and adding ch4 +c ----------------------------------------------------- + else if (modif(1:len_trim(modif)) .eq. 'icarus_ch4') then + + write(*,*) 'At which lat lon do you want to extract the + & reference tsurf/tsoil profile ?' + 416 read(*,*,iostat=ierr) val,val2 + if(ierr.ne.0) goto 416 + write(*,*) 'You want lat =',val + write(*,*) 'You want lon =',val2 + dist=0 + ref=1000 + do j=1,ngridmx-1 + dist=sqrt((lonfi(j)*180./pi-val2)**2+ + & (latfi(j)*180./pi-val)**2) + IF (dist.lt.ref) then + ref=dist + jref=j + ENDIF + ENDDO + + write(*,*)'Will use nearest grid point which is:', + & latfi(jref)*180./pi,lonfi(jref)*180./pi + + ! Extraction of the profile + write(*,*) 'Profile is =',tsoil(jref,:) + write(*,*) 'tsurf is =',tsurf(jref) + write(*,*) 'Choice band : lat1 and lat2 ?' + 417 read(*,*,iostat=ierr) val3,val4 + if(ierr.ne.0) goto 417 + write(*,*) 'Choice temp coefficient for this latitudinal band' + 418 read(*,*,iostat=ierr) val5 + if(ierr.ne.0) goto 418 + write(*,*) 'How much CH4 ice do I put on this band (kg/m2)' + 419 read(*,*,iostat=ierr) choice + if(ierr.ne.0) goto 419 + + DO ig=1,ngridmx + IF ( (latfi(ig)*180./pi.ge.val3) .and. + & (latfi(ig)*180./pi.le.val4) .and. + & (qsurf(ig,igcm_ch4_ice).lt.0.001) ) then + tsurf(ig)=tsurf(jref)*val5 + DO l=1,nsoilmx + tsoil(ig,l)=tsoil(jref,l)*val5 + ENDDO + qsurf(ig,igcm_ch4_ice)=choice + ENDIF + ENDDO + +c globch4co : adding/remove global ch4 co +c ---------------------------------------------------------- + else if (modif(1:len_trim(modif)) .eq. 'globch4co') then + + write(*,*) 'Adding or removing ch4 co ' + write(*,*) 'Choice amount add/less ch4 ice (0 = rm all)' + 431 read(*,*,iostat=ierr) val3 + if(ierr.ne.0) goto 431 + write(*,*) 'Choice amount add/less co ice (0 = rm all)' + 432 read(*,*,iostat=ierr) val4 + if(ierr.ne.0) goto 432 + IF (val3.eq.0.) then + DO ig=1,ngridmx + qsurf(ig,igcm_ch4_ice)=0. + ENDDO + ELSE + DO ig=1,ngridmx + qsurf(ig,igcm_ch4_ice)=qsurf(ig,igcm_ch4_ice)+val3 + ENDDO + ENDIF + IF (val4.eq.0.) then + DO ig=1,ngridmx + qsurf(ig,igcm_co_ice)=0. + ENDDO + ELSE + DO ig=1,ngridmx + qsurf(ig,igcm_co_ice)=qsurf(ig,igcm_co_ice)+val4 + ENDDO + ENDIF + +c bladed : adding/remove ch4 on bladed terrains +c ---------------------------------------------------------- + else if (modif(1:len_trim(modif)).eq.'bladed') then + + write(*,*) 'Adding or removing ch4 on the bladed terrains' + write(*,*) 'Choice band : lat1 and lat2 ?' + 450 read(*,*,iostat=ierr) val,val2 + if(ierr.ne.0) goto 450 + write(*,*) 'Choice band : lon1 and lon2 ?' + 451 read(*,*,iostat=ierr) val3,val4 + if(ierr.ne.0) goto 451 + write(*,*) 'Choice phisfi minimum ?' + 454 read(*,*,iostat=ierr) choice + if(ierr.ne.0) goto 454 + write(*,*) 'Choice multiplicative factor' + 452 read(*,*,iostat=ierr) val5 + if(ierr.ne.0) goto 452 + write(*,*) 'Choice additional ch4' + 453 read(*,*,iostat=ierr) val6 + if(ierr.ne.0) goto 453 + write(*,*) 'Choice index for tsurf tsoil' + 449 read(*,*,iostat=ierr) iref + if(ierr.ne.0) goto 449 + + ! shift + DO ig=1,ngridmx + IF ( (latfi(ig)*180./pi.ge.val) .and. + & (latfi(ig)*180./pi.le.val2) .and. + & (lonfi(ig)*180./pi.ge.val3) .and. + & (lonfi(ig)*180./pi.le.val4) .and. + & (phisfi(ig).gt.choice) ) then + qsurf(ig,igcm_ch4_ice)=qsurf(ig,igcm_ch4_ice)*val5 + qsurf(ig,igcm_ch4_ice)=qsurf(ig,igcm_ch4_ice)+val6 + tsurf(ig)=tsurf(iref) + DO l=1,nsoilmx + tsoil(ig,l)=tsoil(iref,l) + ENDDO + ENDIF + ENDDO + +c cpbladed : add extension bladed terrains +c ---------------------------------------------------------- + else if (modif(1:len_trim(modif)).eq.'cpbladed') then + + write(*,*) 'Choice lat,lon, centre of band to be copied ?' + 560 read(*,*,iostat=ierr) val,val2 + if(ierr.ne.0) goto 560 + write(*,*) 'Choice distance (km) from there defining the area' + 561 read(*,*,iostat=ierr) val3 + if(ierr.ne.0) goto 561 + write(*,*) 'Nb of copies ?' + 562 read(*,*,iostat=ierr) l + if(ierr.ne.0) goto 562 + + ! Get index correponding to central points + ! 1/ Reference point + igref=1 + actualmin=1000. + do ig=1,ngridmx + val6=dist_pluto(latfi(ig),lonfi(ig), + & val*pi/180.,val2*pi/180.) + if (val6.lt.actualmin) then + actualmin=val6 + igref=ig + endif + enddo + + do k=1,l + + write(*,*) 'Choice lat,lon where terrain is copied' + 563 read(*,*,iostat=ierr) val4,val5 + if(ierr.ne.0) goto 563 + + if (val5.gt.180.) then + val5=val5-360. + endif + + ! 2/ Target point + igref2=1 + actualmin=1000. + do ig=1,ngridmx + val6=dist_pluto(latfi(ig),lonfi(ig), + & val4*pi/180.,val5*pi/180.) + if (val6.lt.actualmin) then + actualmin=val6 + igref2=ig + endif + enddo + + ! for each point within A1, get distance and angle + ! save in a table + i=1 + array_ind(:)=0 + array_dist(:)=1000. + array_angle(:)=0. + do ig=1,ngridmx + val6=dist_pluto(latfi(ig),lonfi(ig), + & val*pi/180.,val2*pi/180.) + if (val6.lt.val3) then + array_ind(i)=ig + array_dist(i)=val6 + array_angle(i)=atan2(val-latfi(ig)*180./pi, + & val2-lonfi(ig)*180./pi) + i=i+1 + endif + enddo + + ! for each point within A2, get distance and angle + ! and look where fit with previous table, and set value + + do ig=1,ngridmx + val6=dist_pluto(latfi(ig),lonfi(ig), + & val4*pi/180.,val5*pi/180.) + if (val6.lt.val3) then + ! dist = val6 + ! get angle: + angle=atan2(val4-latfi(ig)*180./pi, + & val5-lonfi(ig)*180./pi) + ! Loop where min + actualmin=1000. + do j=1,i + if ( (array_angle(j).lt.angle+0.52).and. + & (array_angle(j).gt.angle-0.52).and. + & (array_dist(j)-val6.lt.actualmin) ) then + actualmin=array_dist(j)-val6 + igref3=j + endif + enddo + phisfi(ig)=phisfi(array_ind(igref3)) + qsurf(ig,igcm_ch4_ice)= + & qsurf(array_ind(igref3),igcm_ch4_ice) + tsurf(ig)=tsurf(array_ind(igref3)) + endif + enddo + enddo + CALL gr_fi_dyn(1,ngridmx,iip1,jjp1,phisfi,phis) + +c delfrostch4/ delete ch4 frost on a latitudinal band +c ---------------------------------------------------------- + else if (modif(1:len_trim(modif)) .eq. 'delfrostch4') then + + write(*,*) 'removing ch4 latitudinally' + write(*,*) 'Choice band : lat1 and lat2 ?' + read(*,*,iostat=ierr) val,val2 + write(*,*) 'Choice band : lon1 and lon2 ?' + 522 read(*,*,iostat=ierr) val4,val5 + if(ierr.ne.0) goto 522 + write(*,*) 'Choice amount max below whcih ch4 is removed' + 521 read(*,*,iostat=ierr) val3 + if(ierr.ne.0) goto 521 + + DO ig=1,ngridmx + IF ( (latfi(ig)*180./pi.ge.val) .and. + & (latfi(ig)*180./pi.le.val2) .and. + & (lonfi(ig)*180./pi.ge.val4) .and. + & (lonfi(ig)*180./pi.lt.val5) .and. + & (qsurf(ig,igcm_ch4_ice).lt.val3) ) then + qsurf(ig,igcm_ch4_ice)=0. + ENDIF + ENDDO + +c modch4 : adding/remove ch4 frost on a latitudinal band +c ---------------------------------------------------------- + else if (modif(1:len_trim(modif)) .eq. 'modch4') then + + write(*,*) 'Adding or removing ch4 latitudinally' + write(*,*) 'Choice band : lat1 and lat2 ?' + read(*,*,iostat=ierr) val,val2 + write(*,*) 'Choice band : lon1 and lon2 ?' + 422 read(*,*,iostat=ierr) val4,val5 + if(ierr.ne.0) goto 422 + write(*,*) 'Choice multiplicative factor' + 421 read(*,*,iostat=ierr) val3 + if(ierr.ne.0) goto 421 + write(*,*) 'Choice additional ch4' + 423 read(*,*,iostat=ierr) val6 + if(ierr.ne.0) goto 423 + + DO ig=1,ngridmx + IF ( (latfi(ig)*180./pi.ge.val) .and. + & (latfi(ig)*180./pi.le.val2) .and. + & (lonfi(ig)*180./pi.ge.val4) .and. + & (lonfi(ig)*180./pi.lt.val5)) then +! & (qsurf(ig,igcm_n2).gt.50)) then +! & (qsurf(ig,igcm_ch4_ice).lt.10) ) then + qsurf(ig,igcm_ch4_ice)=qsurf(ig,igcm_ch4_ice)*val3 + qsurf(ig,igcm_ch4_ice)=qsurf(ig,igcm_ch4_ice)+val6 + ENDIF + ENDDO + +c modch4n2 : adding/remove ch4 frost on a latitudinal band +c ---------------------------------------------------------- + else if (modif(1:len_trim(modif)) .eq. 'modch4n2') then + + write(*,*) 'Adding or removing ch4 latitudinally' + write(*,*) 'Choice band : lat1 and lat2 ?' + read(*,*,iostat=ierr) val,val2 + write(*,*) 'Choice band : lon1 and lon2 ?' + read(*,*,iostat=ierr) val4,val5 + write(*,*) 'Choice range n2 for modif' + read(*,*,iostat=ierr) val7,val8 + write(*,*) 'Choice multiplicative factor ch4' + read(*,*,iostat=ierr) val3 + write(*,*) 'Choice additional ch4' + read(*,*,iostat=ierr) val6 + + DO ig=1,ngridmx + IF ( (latfi(ig)*180./pi.ge.val) .and. + & (latfi(ig)*180./pi.le.val2) .and. + & (lonfi(ig)*180./pi.ge.val4) .and. + & (lonfi(ig)*180./pi.lt.val5) .and. + & (qsurf(ig,igcm_n2).gt.val7) .and. + & (qsurf(ig,igcm_n2).lt.val8) ) then + qsurf(ig,igcm_ch4_ice)=qsurf(ig,igcm_ch4_ice)*val3 + qsurf(ig,igcm_ch4_ice)=qsurf(ig,igcm_ch4_ice)+val6 + ENDIF + ENDDO + +c non2noco : if no n2 no co +c ---------------------------------------------------------- + else if (modif(1:len_trim(modif)) .eq. 'non2noco') then + DO ig=1,ngridmx + IF ( (qsurf(ig,igcm_n2).lt.0.001)) then + qsurf(ig,igcm_co_ice)=0. + ENDIF + ENDDO + +c modco : adding/remove co frost on a latitudinal band +c ---------------------------------------------------------- + else if (modif(1:len_trim(modif)) .eq. 'modco') then + + write(*,*) 'Adding or removing co where N2 is present ' + write(*,*) 'Choice limit mini n2 pour mettre co?' + 460 read(*,*,iostat=ierr) val7 + if(ierr.ne.0) goto 460 + write(*,*) 'Choice band : lat1 and lat2 ?' + 461 read(*,*,iostat=ierr) val,val2 + if(ierr.ne.0) goto 461 + write(*,*) 'Choice band : lon1 and lon2 ?' + 462 read(*,*,iostat=ierr) val4,val5 + if(ierr.ne.0) goto 462 + write(*,*) 'Choice multiplicative factor' + 463 read(*,*,iostat=ierr) val3 + if(ierr.ne.0) goto 463 + write(*,*) 'Choice additional co' + 464 read(*,*,iostat=ierr) val6 + if(ierr.ne.0) goto 464 + + DO ig=1,ngridmx + IF ( (latfi(ig)*180./pi.ge.val) .and. + & (latfi(ig)*180./pi.le.val2) .and. + & (lonfi(ig)*180./pi.ge.val4) .and. + & (lonfi(ig)*180./pi.lt.val5) .and. + & (qsurf(ig,igcm_n2).lt.val7)) then + qsurf(ig,igcm_co_ice)=qsurf(ig,igcm_co_ice)*val3 + qsurf(ig,igcm_co_ice)=qsurf(ig,igcm_co_ice)+val6 + ENDIF + ENDDO + +c modn2 : adding/remove n2 ice +c ---------------------------------------------------------- + else if (modif(1:len_trim(modif)) .eq. 'modn2') then + + write(*,*) 'Adding or removing n2 ' + write(*,*) 'Choice band : lat1 and lat2 ?' + 425 read(*,*,iostat=ierr) val,val2 + if(ierr.ne.0) goto 425 + write(*,*) 'Choice band : lon1 and lon2 ?' + 426 read(*,*,iostat=ierr) val4,val5 + if(ierr.ne.0) goto 426 + write(*,*) 'Choice multiplicative factor' + 427 read(*,*,iostat=ierr) val3 + if(ierr.ne.0) goto 427 + write(*,*) 'Choice additional n2' + 428 read(*,*,iostat=ierr) val6 + if(ierr.ne.0) goto 428 + + DO ig=1,ngridmx + IF ( (latfi(ig)*180./pi.ge.val) .and. + & (latfi(ig)*180./pi.le.val2) .and. + & (lonfi(ig)*180./pi.ge.val4) .and. + & (lonfi(ig)*180./pi.lt.val5) ) then +c & (qsurf(ig,igcm_n2).lt.50)) then + qsurf(ig,igcm_n2)=qsurf(ig,igcm_n2)*val3 + qsurf(ig,igcm_n2)=qsurf(ig,igcm_n2)+val6 + ENDIF +! IF ((phisfi(ig).gt.-1000.)) then +! qsurf(ig,igcm_n2)=0. +! ENDIF + ENDDO + +c modcoch4 : adding/remove co and ch4 frost where co without n2 +c ---------------------------------------------------------- + else if (modif(1:len_trim(modif)) .eq. 'modcoch4') then + + write(*,*) 'Adding or removing co where N2 is not present ' + write(*,*) 'Choice band : lat1 and lat2 ?' + 491 read(*,*,iostat=ierr) val,val2 + if(ierr.ne.0) goto 491 + write(*,*) 'Choice band : lon1 and lon2 ?' + 492 read(*,*,iostat=ierr) val4,val5 + if(ierr.ne.0) goto 492 + write(*,*) 'Choice multiplicative factor' + 493 read(*,*,iostat=ierr) val3 + if(ierr.ne.0) goto 493 + write(*,*) 'Choice additional co ch4' + 494 read(*,*,iostat=ierr) val6 + if(ierr.ne.0) goto 494 + + DO ig=1,ngridmx + IF ( (latfi(ig)*180./pi.ge.val) .and. + & (latfi(ig)*180./pi.le.val2) .and. + & (lonfi(ig)*180./pi.ge.val4) .and. + & (lonfi(ig)*180./pi.le.val5) .and. + & (qsurf(ig,igcm_n2).lt.10.)) then + qsurf(ig,igcm_co_ice)=qsurf(ig,igcm_co_ice)*val3 + qsurf(ig,igcm_co_ice)=qsurf(ig,igcm_co_ice)+val6 + qsurf(ig,igcm_ch4_ice)=qsurf(ig,igcm_ch4_ice)*val3 + qsurf(ig,igcm_ch4_ice)=qsurf(ig,igcm_ch4_ice)+val6 + ENDIF + ENDDO + +c modn2HD : adding/remove n2 ice for HD runs +c ---------------------------------------------------------- + else if (modif(1:len_trim(modif)) .eq. 'modn2HD') then + + write(*,*) 'Adding or removing n2 ' + write(*,*) 'Choice band : lat1 and lat2 ?' + 480 read(*,*,iostat=ierr) val,val1 + if(ierr.ne.0) goto 480 + write(*,*) 'Choice band : lon1 and lon2 ?' + 481 read(*,*,iostat=ierr) val2,val3 + if(ierr.ne.0) goto 481 + write(*,*) 'Choice amount N2 min max where to modify' + 482 read(*,*,iostat=ierr) val4,val5 + if(ierr.ne.0) goto 482 + write(*,*) 'Choice phisfi min max where change n2' + 483 read(*,*,iostat=ierr) val6,val11 + if(ierr.ne.0) goto 483 + write(*,*) 'Choice multiplicative factor' + 484 read(*,*,iostat=ierr) val7 + if(ierr.ne.0) goto 484 + write(*,*) 'Choice additional n2' + 485 read(*,*,iostat=ierr) val8 + if(ierr.ne.0) goto 485 +! write(*,*) 'Choice ind lon equivalent for change tsurf tsoil' +! 486 read(*,*,iostat=ierr) val9 +! if(ierr.ne.0) goto 486 + + DO ig=1,ngridmx + IF ( (latfi(ig)*180./pi.ge.val) .and. + & (latfi(ig)*180./pi.le.val1) .and. + & (lonfi(ig)*180./pi.ge.val2) .and. + & (lonfi(ig)*180./pi.lt.val3) .and. + & (qsurf(ig,igcm_n2).ge.val4) .and. + & (qsurf(ig,igcm_n2).lt.val5) .and. + & (phisfi(ig).ge.val6) .and. + & (phisfi(ig).le.val11) ) then + + qsurf(ig,igcm_n2)=qsurf(ig,igcm_n2)*val7 + qsurf(ig,igcm_n2)=qsurf(ig,igcm_n2)+val8 + !qsurf(ig,igcm_ch4_ice)=0. + qsurf(ig,igcm_co_ice)=0. + ! index of ig + !if (val9.eq.0.) then + ! iref=2546 + !else + ! val10=modulo((ig-1),96) + ! choice=ig+int(96/2)-val10 + ! !choice=int(1+96*(val9-1)+val10) + !endif + !tsurf(ig)=tsurf(iref) + !DO l=1,nsoilmx + ! tsoil(ig,l)=tsoil(iref,l) + !ENDDO + tsurf(ig)=tsurf(ig)+4 + DO l=1,nsoilmx + tsoil(ig,l)=tsoil(ig,l)+4 + ENDDO + ENDIF +! IF ((phisfi(ig).gt.-1000.)) then +! qsurf(ig,igcm_n2)=0. +! ENDIF + ENDDO + +c modn2HD_SP : adding/remove n2 ice for HD runs +c ---------------------------------------------------------- + else if (modif(1:len_trim(modif)) .eq. 'modn2HD_SP') then + + write(*,*) 'Adding or removing n2 ' + write(*,*) 'Choice band : lat1 and lat2 ?' + 501 read(*,*,iostat=ierr) val,val1 + if(ierr.ne.0) goto 501 + write(*,*) 'Choice band : lon1 and lon2 ?' + 502 read(*,*,iostat=ierr) val2,val3 + if(ierr.ne.0) goto 502 + write(*,*) 'Choice amount N2 min max where to modify' + 503 read(*,*,iostat=ierr) val4,val5 + if(ierr.ne.0) goto 503 + write(*,*) 'Choice phisfi min max where change n2' + 504 read(*,*,iostat=ierr) val6,val11 + if(ierr.ne.0) goto 504 + write(*,*) 'Choice multiplicative factor' + 505 read(*,*,iostat=ierr) val7 + if(ierr.ne.0) goto 505 + write(*,*) 'Choice additional n2' + 506 read(*,*,iostat=ierr) val8 + if(ierr.ne.0) goto 506 + write(*,*) 'Choice ind lon equivalent for change tsurf tsoil' + 507 read(*,*,iostat=ierr) iref + if(ierr.ne.0) goto 507 + + DO ig=1,ngridmx + IF ( (latfi(ig)*180./pi.ge.val) .and. + & (latfi(ig)*180./pi.le.val1) .and. + & (lonfi(ig)*180./pi.ge.val2) .and. + & (lonfi(ig)*180./pi.lt.val3) .and. + & (qsurf(ig,igcm_n2).ge.val4) .and. + & (qsurf(ig,igcm_n2).lt.val5) .and. + & (phisfi(ig).ge.val6) .and. + & (phisfi(ig).le.val11) ) then + + qsurf(ig,igcm_n2)=qsurf(ig,igcm_n2)*val7 + qsurf(ig,igcm_n2)=qsurf(ig,igcm_n2)+val8 + qsurf(ig,igcm_ch4_ice)=0. + qsurf(ig,igcm_co_ice)=0. + ! index of ig + !if (val9.eq.0.) then + ! iref=2546 + !else + !val10=modulo((ig-1),96) + !choice=ig+int(96/2)-val10 + !choice=int(1+96*(val9-1)+val10) + if (iref.ge.1) then + tsurf(ig)=tsurf(iref) + DO l=1,nsoilmx + tsoil(ig,l)=tsoil(iref,l) + ENDDO + else if (iref.eq.0) then + choice=int(ig/96.)*96+84 + print*, 'choice=',choice + tsurf(ig)=tsurf(int(choice)) + DO l=1,nsoilmx + tsoil(ig,l)=tsoil(int(choice),l) + ENDDO + endif + ENDIF + ENDDO + + +c modn2_topo : adding/remove n2 ice +c ---------------------------------------------------------- + else if (modif(1:len_trim(modif)) .eq. 'modn2_topo') then + + write(*,*) 'Adding or removing n2 ' + write(*,*) 'Choice band : lat1 and lat2 ?' + 441 read(*,*,iostat=ierr) val,val2 + if(ierr.ne.0) goto 441 + write(*,*) 'Choice band : lon1 and lon2 ?' + 442 read(*,*,iostat=ierr) val4,val5 + if(ierr.ne.0) goto 442 + write(*,*) 'Choice topo 1 et 2 (phi) where change is active' + 443 read(*,*,iostat=ierr) val7,val8 + if(ierr.ne.0) goto 443 + write(*,*) 'Choice multiplicative factor' + 444 read(*,*,iostat=ierr) val3 + if(ierr.ne.0) goto 444 + write(*,*) 'Choice additional n2' + 445 read(*,*,iostat=ierr) val6 + if(ierr.ne.0) goto 445 + + DO ig=1,ngridmx + IF ( (latfi(ig)*180./pi.ge.val) .and. + & (latfi(ig)*180./pi.le.val2) .and. + & (lonfi(ig)*180./pi.ge.val4) .and. + & (lonfi(ig)*180./pi.lt.val5) .and. + & (phisfi(ig).le.val8) .and. + & (phisfi(ig).ge.val7) ) then + qsurf(ig,igcm_n2)=qsurf(ig,igcm_n2)*val3 + qsurf(ig,igcm_n2)=qsurf(ig,igcm_n2)+val6 + ENDIF + ENDDO + +c modwheren2 : adding/remove n2 ice where already n2 +c ---------------------------------------------------------- + else if (modif(1:len_trim(modif)) .eq. 'modwheren2') then + + write(*,*) 'Choice multiplicative factor' + 471 read(*,*,iostat=ierr) val3 + if(ierr.ne.0) goto 471 + write(*,*) 'Choice additional n2' + 472 read(*,*,iostat=ierr) val6 + if(ierr.ne.0) goto 472 + + DO ig=1,ngridmx + IF ( qsurf(ig,igcm_n2).gt.0.001 ) then + qsurf(ig,igcm_n2)=qsurf(ig,igcm_n2)*val3 + qsurf(ig,igcm_n2)=qsurf(ig,igcm_n2)+val6 + ENDIF + ENDDO + +c checktsoil : changing tsoil if no n2 +c ---------------------------------------------------------- + else if (modif(1:len_trim(modif)) .eq. 'checktsoil') then + + write(*,*) 'selecting area where tsoil to be check ' + write(*,*) 'Choice band : lat1 and lat2 ?' + 511 read(*,*,iostat=ierr) val,val1 + if(ierr.ne.0) goto 511 + write(*,*) 'Choice band : lon1 and lon2 ?' + 512 read(*,*,iostat=ierr) val2,val3 + if(ierr.ne.0) goto 512 + write(*,*) 'Choice temp min max in which change is made' + 513 read(*,*,iostat=ierr) val4,val5 + if(ierr.ne.0) goto 513 + write(*,*) 'Choice phisfi min max where change n2' + 514 read(*,*,iostat=ierr) val6,val11 + if(ierr.ne.0) goto 514 + + DO ig=1,ngridmx + IF ( (latfi(ig)*180./pi.ge.val) .and. + & (latfi(ig)*180./pi.le.val1) .and. + & (lonfi(ig)*180./pi.ge.val2) .and. + & (lonfi(ig)*180./pi.le.val3) .and. + & (((tsurf(ig).ge.val4) .and. + & (tsurf(ig).le.val5)) .or. + & ((tsoil(ig,nsoilmx).ge.val4) .and. + & (tsoil(ig,nsoilmx).le.val5))) .and. + & (phisfi(ig).ge.val6) .and. + & (phisfi(ig).le.val11) .and. + & (qsurf(ig,igcm_n2).le.0.001) ) then + +! DO i=1,ngridmx +! IF ( (latfi(i).eq.latfi(ig)) .and. +! & (lonfi(i).eq.0.) ) then +! + tsurf(ig)=50. + !qsurf(ig,igcm_ch4_ice)=qsurf(i,igcm_ch4_ice) +! + DO l=1,nsoilmx + tsoil(ig,l)=50. !tsoil(i,l) + ENDDO + !ENDIF + !ENDDO + ENDIF + ENDDO + +c asymtriton : changing ice, tsurf and tsoil on triton +c ---------------------------------------------------------- + else if (modif(1:len_trim(modif)) .eq. 'asymtriton') then + + write(*,*) 'selecting area where tsoil to be check ' + write(*,*) 'Choice center latitude of sinuisoid distrib?' + 531 read(*,*,iostat=ierr) val1 + if(ierr.ne.0) goto 531 + write(*,*) 'Choice maginutde in latitude?' + 532 read(*,*,iostat=ierr) val2 + if(ierr.ne.0) goto 532 + write(*,*) 'Choice center longitude?' + 533 read(*,*,iostat=ierr) val3 + if(ierr.ne.0) goto 533 +! write(*,*) 'iip1,jjp1=',iip1,jjp1,ngridmx + + DO ig=1,ngridmx + ! Latitude of the sinusoid: + val11=val1+val2*cos(lonfi(ig)*1.57079633*2/pi-val3*pi/180.) + ! If above line and ice: remove ice + IF ( (latfi(ig)*180./pi.ge.val11) .and. + & (latfi(ig)*180./pi.le.val1+val2) .and. + & (qsurf(ig,igcm_n2).gt.0.) ) then + ! Look for same longitude point where no ice + val5=1. + jref=ig + ! 1 + ! ... iip1 ... x (jjp1-2) 32 x 23 + ! 1 + do while (val5.ge.1..and.jref.gt.iip1) + ! northward point + jref=jref-iip1+1 + ! ice at the point + val5=qsurf(jref,igcm_n2) +! write(*,*) jref, +! & latfi(jref)*180./pi,lonfi(jref)*180/pi,val5 + enddo + if (val5.ge.1) write(*,*) 'NO POINT FOUND WITH NO ICE' + + ! Copy point in the selected area + tsurf(ig)=tsurf(jref) + qsurf(ig,igcm_n2)=qsurf(jref,igcm_n2) + qsurf(ig,igcm_ch4_ice)=qsurf(jref,igcm_ch4_ice) + qsurf(ig,igcm_co_ice)=qsurf(jref,igcm_co_ice) + DO l=1,nsoilmx + tsoil(ig,l)=tsoil(jref,l) + ENDDO + ENDIF + ! If below line and no ice: add ice + IF ( (latfi(ig)*180./pi.le.val11) .and. + & (latfi(ig)*180./pi.le.val1+val2) .and. + & (qsurf(ig,igcm_n2).eq.0.) ) then + ! Look for same longitude point where ice + val5=1. + jref=ig + do while (val5.le.1.and.jref.lt.ngridmx-iip1) + ! southward point + jref=jref+iip1-1 + ! ice at the point + val5=qsurf(jref,igcm_n2) + write(*,*) jref, + & latfi(jref)*180./pi,lonfi(jref)*180/pi,val5 + enddo + if (val5.le.1) write(*,*) 'NO POINT FOUND WITH ICE' + + ! Copy point in the selected area + tsurf(ig)=tsurf(jref) + qsurf(ig,igcm_n2)=qsurf(jref,igcm_n2) + qsurf(ig,igcm_ch4_ice)=qsurf(jref,igcm_ch4_ice) + qsurf(ig,igcm_co_ice)=qsurf(jref,igcm_co_ice) + DO l=1,nsoilmx + tsoil(ig,l)=tsoil(jref,l) + ENDDO + ENDIF + + ENDDO + +c source_ch4 : adding source ch4 +c ---------------------------------------------------------- + else if (modif(1:len_trim(modif)) .eq. 'source_ch4') then + + write(*,*) 'Adding ch4 ice latitudinally ' + write(*,*) 'Choice : lat1 and lat2 ?' + 433 read(*,*,iostat=ierr) val,val2 + if(ierr.ne.0) goto 433 + write(*,*) 'Choice : lon1 and lon2 ?' + 434 read(*,*,iostat=ierr) val3,val4 + if(ierr.ne.0) goto 434 + write(*,*) 'Choice amount (kg/m2)' + 435 read(*,*,iostat=ierr) val5 + if(ierr.ne.0) goto 435 + + DO ig=1,ngridmx + IF ( (latfi(ig)*180./pi.ge.val) .and. + & (latfi(ig)*180./pi.le.val2) .and. + & (lonfi(ig)*180./pi.ge.val3) .and. + & (lonfi(ig)*180./pi.lt.val4) ) then + qsurf(ig,igcm_ch4_ice)=qsurf(ig,igcm_ch4_ice)+val5 + ENDIF + ENDDO + +c source_co : adding source co +c ---------------------------------------------------------- + else if (modif(1:len_trim(modif)) .eq. 'source_co') then + + write(*,*) 'Adding co ice latitudinally ' + write(*,*) 'Choice : lat1 and lat2 ?' + 436 read(*,*,iostat=ierr) val,val2 + if(ierr.ne.0) goto 436 + write(*,*) 'Choice : lon1 and lon2 ?' + 437 read(*,*,iostat=ierr) val3,val4 + if(ierr.ne.0) goto 437 + write(*,*) 'Choice amount (kg/m2)' + 438 read(*,*,iostat=ierr) val5 + if(ierr.ne.0) goto 438 + + DO ig=1,ngridmx + IF ( (latfi(ig)*180./pi.ge.val) .and. + & (latfi(ig)*180./pi.le.val2) .and. + & (lonfi(ig)*180./pi.ge.val3) .and. + & (lonfi(ig)*180./pi.lt.val4) ) then + qsurf(ig,igcm_co_ice)=qsurf(ig,igcm_co_ice)+val5 + ENDIF + ENDDO + +! therm_ini_s: (re)-set soil thermal inertia to reference surface values +! ---------------------------------------------------------------------- + else if (modif(1:len_trim(modif)).eq.'therm_ini_s') then +! write(*,*)"surfithfi(1):",surfithfi(1) + do isoil=1,nsoilmx + inertiedat(1:ngridmx,isoil)=surfithfi(1:ngridmx) + enddo + write(*,*)'OK: Soil thermal inertia has been reset to referenc + &e surface values' + write(*,*)"inertiedat(1,1):",inertiedat(1,1) + ithfi(:,:)=inertiedat(:,:) + ! recast ithfi() onto ith() + call gr_fi_dyn(nsoilmx,ngridmx,iip1,jjp1,ithfi,ith) + +! inert3d: set soil thermal inertia to specific uniform value +! ---------------------------------------------------------------------- + else if (modif(1:len_trim(modif)).eq.'inert3d') then + write(*,*) 'Actual value of surf Thermal inertia at ig=1: ', + & inertiedat(1,1), ' SI' + write(*,*) 'Value of Thermal inertia (SI) ? ' + read(*,*) val + do isoil=1,nsoilmx + do ig=1,ngridmx + inertiedat(ig,isoil)=val + enddo + enddo + write(*,*)'OK: Soil TI set to ',val,' SI everywhere' + ithfi(:,:)=inertiedat(:,:) + ! recast ithfi() onto ith() + call gr_fi_dyn(nsoilmx,ngridmx,iip1,jjp1,ithfi,ith) + +! subsoilice_n: Put deep ice layer in northern hemisphere soil +! ------------------------------------------------------------ + else if (modif(1:len_trim(modif)).eq.'subsoilice_n') then + + write(*,*)'From which latitude (in deg.), up to the north pole, + &should we put subterranean ice?' + ierr=1 + do while (ierr.ne.0) + read(*,*,iostat=ierr) val + if (ierr.eq.0) then ! got a value + ! do a sanity check + if((val.lt.0.).or.(val.gt.90)) then + write(*,*)'Latitude should be between 0 and 90 deg. !!!' ierr=1 - 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+1 - write(*,*)'Level selected : ',layer(isoil+1),' m' - exit - endif - enddo - endif - - DO ig=1,ngridmx - DO j=iref,nsoilmx - ithfi(ig,j)=ith_bb - ENDDO - ENDDO - - CALL gr_fi_dyn(nsoilmx,ngridmx,iip1,jjp1,ithfi,ith) - -c diurnal_TI : choice of thermal inertia values (global) -c ---------------------------------------------------------------- - else if (modif(1:len_trim(modif)) .eq. 'diurnal_TI') then - - write(*,*) 'New value for diurnal thermal inertia ?' - 106 read(*,*,iostat=ierr) ith_bb - if(ierr.ne.0) goto 106 - write(*,*) 'Diurnal thermal inertia (new value):',ith_bb - - write(*,*)'At which depth (in m.) does the ice layer ends?' - write(*,*)'(currently, the soil layer 1 and nsoil are:' - & ,layer(1),' - ',layer(nsoilmx),')' + 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 in diurnal_TI:',val2,'ierr=',ierr +! write(*,*)'val2:',val2,'ierr=',ierr if (ierr.eq.0) then ! got a value, but do a sanity check if(val2.gt.layer(nsoilmx)) then @@ -911,47 +1900,323 @@ endif enddo ! of do while - - ! find the reference index iref the depth corresponds to + + ! find the reference index iref the depth corresponds to +! if (val2.lt.layer(1)) then +! iref=1 +! 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))) + if((val2.gt.layer(isoil)).and.(val2.lt.layer(isoil+1))) & then - iref=isoil+1 - write(*,*)'Level selected : ',layer(isoil+1),' m' + iref=isoil exit + endif + enddo + + ! thermal inertia of the ice: + ierr=1 + do while (ierr.ne.0) + write(*,*)'What is the value of subterranean ice thermal inert + &ia? (e.g.: 2000)' + read(*,*,iostat=ierr)iceith + enddo ! of do while + + ! recast ithfi() onto ith() + call gr_fi_dyn(nsoilmx,ngridmx,iip1,jjp1,ithfi,ith) + + do j=1,jref +! write(*,*)'j:',j,'rlatu(j)*180./pi:',rlatu(j)*180./pi + do i=1,iip1 ! loop on longitudes + ! Build "equivalent" thermal inertia for the mixed layer + ith(i,j,iref+1)=sqrt((layer(iref+1)-layer(iref))/ + & (((val2-layer(iref))/(ith(i,j,iref)**2))+ + & ((layer(iref+1)-val2)/(iceith)**2))) + ! Set thermal inertia of lower layers + do isoil=iref+2,nsoilmx + ith(i,j,isoil)=iceith ! ice + enddo + enddo ! of do i=1,iip1 + enddo ! of do j=1,jjp1 + + + CALL gr_dyn_fi(nsoilmx,iip1,jjp1,ngridmx,ith,ithfi) + +! subsoilice_s: Put deep ice layer in southern hemisphere soil +! ------------------------------------------------------------ + else if (modif(1:len_trim(modif)).eq.'subsoilice_s') then + + write(*,*)'From which latitude (in deg.), down to the south pol + &e, should we put subterranean ice?' + ierr=1 + do while (ierr.ne.0) + read(*,*,iostat=ierr) val + if (ierr.eq.0) then ! got a value + ! do a sanity check + if((val.gt.0.).or.(val.lt.-90)) then + write(*,*)'Latitude should be between 0 and -90 deg. !!!' + ierr=1 + else ! find corresponding jref (nearest latitude) + ! note: rlatu(:) contains decreasing values of latitude + ! starting from PI/2 to -PI/2 + do j=1,jjp1 + if ((rlatu(j)*180./pi.ge.val).and. + & (rlatu(j+1)*180./pi.le.val)) then + ! find which grid point is nearest to val: + if (abs(rlatu(j)*180./pi-val).le. + & abs((rlatu(j+1)*180./pi-val))) then + jref=j + else + jref=j+1 + endif + + write(*,*)'Will use nearest grid latitude which is:', + & rlatu(jref)*180./pi + endif + enddo ! of do j=1,jjp1 + endif ! of if((val.lt.0.).or.(val.gt.90)) + endif !of if (ierr.eq.0) + enddo ! of do while + + ! Build layers() (as in soil_settings.F) + val2=sqrt(mlayer(0)*mlayer(1)) + val3=mlayer(1)/mlayer(0) + do isoil=1,nsoilmx + layer(isoil)=val2*(val3**(isoil-1)) + enddo + + write(*,*)'At which depth (in m.) does the ice layer begin?' + write(*,*)'(currently, the deepest soil layer extends down to:' + & ,layer(nsoilmx),')' + ierr=1 + do while (ierr.ne.0) + read(*,*,iostat=ierr) val2 +! write(*,*)'val2:',val2,'ierr=',ierr + if (ierr.eq.0) then ! got a value, but do a sanity check + if(val2.gt.layer(nsoilmx)) then + write(*,*)'Depth should be less than ',layer(nsoilmx) + ierr=1 endif + if(val2.lt.layer(1)) then + write(*,*)'Depth should be more than ',layer(1) + ierr=1 + endif + endif + enddo ! of do while + + ! find the reference index iref the depth corresponds to + do isoil=1,nsoilmx-1 + if((val2.gt.layer(isoil)).and.(val2.lt.layer(isoil+1))) + & then + iref=isoil + exit + endif + enddo + + ! thermal inertia of the ice: + ierr=1 + do while (ierr.ne.0) + write(*,*)'What is the value of subterranean ice thermal inert + &ia? (e.g.: 2000)' + read(*,*,iostat=ierr)iceith + enddo ! of do while + + ! recast ithfi() onto ith() + call gr_fi_dyn(nsoilmx,ngridmx,iip1,jjp1,ithfi,ith) + + do j=jref,jjp1 +! write(*,*)'j:',j,'rlatu(j)*180./pi:',rlatu(j)*180./pi + do i=1,iip1 ! loop on longitudes + ! Build "equivalent" thermal inertia for the mixed layer + ith(i,j,iref+1)=sqrt((layer(iref+1)-layer(iref))/ + & (((val2-layer(iref))/(ith(i,j,iref)**2))+ + & ((layer(iref+1)-val2)/(iceith)**2))) + ! Set thermal inertia of lower layers + do isoil=iref+2,nsoilmx + ith(i,j,isoil)=iceith ! ice + enddo + enddo ! of do i=1,iip1 + enddo ! of do j=jref,jjp1 + + CALL gr_dyn_fi(nsoilmx,iip1,jjp1,ngridmx,ith,ithfi) + +c ---------------------------------------------------------- +c 'reservoir_SP' : increase or decrase ices reservoir in SP by factor +c ---------------------------------------------------------- + else if (modif(1:len_trim(modif)).eq.'reservoir_SP') then + write(*,*) "Increase/Decrease N2 reservoir by factor :" + read(*,*) val7 + write(*,*) "Increase/Decrease CH4 reservoir by factor :" + read(*,*) val8 + write(*,*) "Increase/Decrease CO reservoir by factor :" + read(*,*) val9 + + ! Definition SP: + val3=-33. !lat1 + val4=50. !lat2 + val5=140. ! lon1 + val6=-155. ! lon2 + choice=-50. ! min phisfi + write(*,*) 'def SP :' + write(*,*) 'lat :',val3,val4 + write(*,*) 'lon :',val5,'180 / -180',val6 + write(*,*) 'choice phisfi min :',choice + + DO ig=1,ngridmx + + IF ( (latfi(ig)*180./pi.ge.val3) .and. + & (latfi(ig)*180./pi.le.val4) .and. + & ( ((lonfi(ig)*180./pi.ge.-180.) .and. + & (lonfi(ig)*180./pi.le.val6)) .or. + & ((lonfi(ig)*180./pi.ge.val5) .and. + & (lonfi(ig)*180./pi.le.180.))) ) then + + IF ((phisfi(ig).le.choice)) then !.and. +c & (qsurf(ig,igcm_n2).ge.50)) then + + qsurf(ig,igcm_n2)=qsurf(ig,igcm_n2)*val7 + qsurf(ig,igcm_ch4_ice)=qsurf(ig,igcm_ch4_ice)*val8 + qsurf(ig,igcm_co_ice)=qsurf(ig,igcm_co_ice)*val9 + ENDIF + ENDIF + ENDDO + +c ---------------------------------------------------------- +c 'reservoir' : increase or decrase ices reservoir by factor +c ---------------------------------------------------------- + else if (modif(1:len_trim(modif)).eq.'reservoir') then + write(*,*) "Increase/Decrease N2 reservoir by factor :" + read(*,*) val + write(*,*) "Increase/Decrease CH4 reservoir by factor :" + read(*,*) val2 + write(*,*) "Increase/Decrease CO reservoir by factor :" + read(*,*) val3 + + DO ig=1,ngridmx + qsurf(ig,igcm_n2)=qsurf(ig,igcm_n2)*val + qsurf(ig,igcm_ch4_ice)=qsurf(ig,igcm_ch4_ice)*val2 + qsurf(ig,igcm_co_ice)=qsurf(ig,igcm_co_ice)*val3 + ENDDO + +c -------------------------------------------------------- +c 'plutomap' : initialize pluto ices distribution +c -------------------------------------------------------- + else if (modif(1:len_trim(modif)).eq.'plutomap') then + + write(*,*) 'pluto_map.dat has to be in your simulation file !!' + open(27,file='pluto_map.dat') + N2_pluto_dat(:,:) =0. + CH4_pluto_dat(:,:) =0. + CO_pluto_dat(:,:) =0. + + ! Dimension file pluto_map + IF (jm_plu.ne.180) then + write(*,*) 'Newstart : you should set jm_plu to 180' + stop + ENDIF + ! Read values + do j=1,jm_plu+1 + read(27,*) (map_pluto_dat(i,j),i=1,im_plu) + do i=1,im_plu + + !!!!!! BTD_v2 + if (map_pluto_dat(i,j).eq.3) then + N2_pluto_dat(i,j)=100000. + elseif (map_pluto_dat(i,j).eq.4) then + CH4_pluto_dat(i,j)=100000. + elseif (map_pluto_dat(i,j).eq.0) then + alb_dat(i,j)=0.3 + elseif (map_pluto_dat(i,j).eq.6) then + alb_dat(i,j)=0.6 + elseif (map_pluto_dat(i,j).eq.7) then + alb_dat(i,j)=0.1 + endif + + !!!!!! BTD + !if (map_pluto_dat(i,j).eq.3) then + ! CH4_pluto_dat(i,j)=100000. + !endif + end do + end do + close (27) + ! Interpolate + + !! latitude and longitude in REindexed pluto map : + latv_plu(1)=90. *pi/180. + do j=2,jm_plu + latv_plu(j)= (90-(j-1 -0.5)*180./(jm_plu-1))*pi/180. + end do + latv_plu(jm_plu+1)= -90. *pi/180. + + do i=1,im_plu+1 + lonu_plu(i)=(-180+ (i-1)*360./(im_plu) ) *pi/180. enddo + !Reindexation to shift the longitude grid like a LMD GCM grid... + do j=1,jm_plu+1 + N2_pluto_rein(1,j)=(N2_pluto_dat(1,j)+ + & N2_pluto_dat(im_plu,j))/2 + CH4_pluto_rein(1,j)=(CH4_pluto_dat(1,j)+ + & CH4_pluto_dat(im_plu,j))/2 + CO_pluto_rein(1,j)=(CO_pluto_dat(1,j)+ + & CO_pluto_dat(im_plu,j))/2 + alb_rein(1,j)=(alb_dat(1,j)+ + & alb_dat(im_plu,j))/2 + do i=2,im_plu + N2_pluto_rein(i,j)=(N2_pluto_dat(i-1,j)+ + & N2_pluto_dat(i,j))/2 + CH4_pluto_rein(i,j)=(CH4_pluto_dat(i-1,j)+ + & CH4_pluto_dat(i,j))/2 + CO_pluto_rein(i,j)=(CO_pluto_dat(i-1,j)+ + & CO_pluto_dat(i,j))/2 + alb_rein(i,j)=(alb_dat(i-1,j)+ + & alb_dat(i,j))/2 + end do + N2_pluto_rein(im_plu+1,j) = N2_pluto_rein(1,j) + CH4_pluto_rein(im_plu+1,j) = CH4_pluto_rein(1,j) + CO_pluto_rein(im_plu+1,j) = CO_pluto_rein(1,j) + alb_rein(im_plu+1,j) = alb_rein(1,j) + end do + + call interp_horiz(N2_pluto_rein,N2_pluto_gcm, + & im_plu,jm_plu,iim,jjm,1,lonu_plu,latv_plu,rlonu,rlatv) + call interp_horiz(CH4_pluto_rein,CH4_pluto_gcm, + & im_plu,jm_plu,iim,jjm,1,lonu_plu,latv_plu,rlonu,rlatv) + call interp_horiz(CO_pluto_rein,CO_pluto_gcm, + & im_plu,jm_plu,iim,jjm,1,lonu_plu,latv_plu,rlonu,rlatv) + call interp_horiz(alb_rein,alb_gcm, + & im_plu,jm_plu,iim,jjm,1,lonu_plu,latv_plu,rlonu,rlatv) + + ! Switch dyn grid to fi grid + qsurf_tmp(:,:) =0. + CALL gr_dyn_fi(1,iip1,jjp1,ngridmx, + & N2_pluto_gcm,qsurf_tmp(1,igcm_n2)) + if (igcm_ch4_ice.ne.0) then + CALL gr_dyn_fi(1,iip1,jjp1,ngridmx, + & CH4_pluto_gcm,qsurf_tmp(1,igcm_ch4_ice)) + endif + if (igcm_co_ice.ne.0) then + CALL gr_dyn_fi(1,iip1,jjp1,ngridmx, + & CO_pluto_gcm,qsurf_tmp(1,igcm_co_ice)) + endif + alb=alb_gcm + CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,alb_gcm,albfi) + !print*, 'N2dat=',N2_pluto_gcm + !print*, 'COdat=',CO_pluto_gcm + !print*, 'CH4dat=',CH4_pluto_gcm + print*, 'N2dat=',qsurf_tmp(:,igcm_n2) + print*, 'COdat=',qsurf_tmp(:,igcm_co_ice) + print*, 'CH4dat=',qsurf_tmp(:,igcm_ch4_ice) + + ! Fill DO ig=1,ngridmx - DO j=1,iref - ithfi(ig,j)=ith_bb - ENDDO - ENDDO - - CALL gr_fi_dyn(nsoilmx,ngridmx,iip1,jjp1,ithfi,ith) - - - - - - - - - - - - - - - - - - - - - - - + qsurf(ig,igcm_n2)=qsurf(ig,igcm_n2)+qsurf_tmp(ig,igcm_n2) + if (igcm_ch4_ice.ne.0) then + qsurf(ig,igcm_ch4_ice)=qsurf(ig,igcm_ch4_ice)+ + & qsurf_tmp(ig,igcm_ch4_ice) + endif + if (igcm_co_ice.ne.0) then + qsurf(ig,igcm_co_ice)=qsurf(ig,igcm_co_ice)+ + & qsurf_tmp(ig,igcm_co_ice) + endif + ENDDO endif @@ -961,33 +2226,7 @@ 999 continue - -c======================================================================= -c Initialisation for cloud fraction and oceanic ice (added by BC 2010) -c======================================================================= -! DO ig=1, ngridmx -! totalfrac(ig)=0.5 -! DO l=1,llm -! cloudfrac(ig,l)=0.5 -! ENDDO -! Ehouarn, march 2012: also add some initialisation for hice -! hice(ig)=0.0 -! ENDDO - -c======================================================== - -! DO ig=1,ngridmx -! IF(tsurf(ig) .LT. 273.16-1.8) THEN -! hice(ig)=(273.16-1.8-tsurf(ig))/(273.16-1.8-240)*1. -! hice(ig)=min(hice(ig),1.0) -! ENDIF -! ENDDO - - - - c======================================================================= c Correct pressure on the new grid (menu 0) c======================================================================= - if ((choix_1.eq.0).and.(.not.flagps0)) then @@ -1008,12 +2247,7 @@ end if - -c======================================================================= -c======================================================================= - c======================================================================= c Initialisation de la physique / ecriture de newstartfi : c======================================================================= - CALL inifilr @@ -1111,2 +2345,20 @@ 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 +