Changeset 3371 for trunk/LMDZ.PLUTO/libf
- Timestamp:
- Jun 12, 2024, 8:43:20 PM (6 months ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/LMDZ.PLUTO/libf/dynphy_lonlat/phypluto/newstart.F
r3228 r3371 6 6 c Auteur: Christophe Hourdin/Francois Forget/Yann Wanherdrick 7 7 c ------ 8 c Derniere modif : 12/039 c 8 c Adapted to Pluto: Tanguy Bertrand 9 c Derniere modif : 06/2024 10 10 c 11 11 c Objet: Create or modify the initial state for the LMD Mars GCM 12 12 c ----- (fichiers NetCDF start et startfi) 13 13 c 14 c15 14 c======================================================================= 16 15 … … 19 18 & is_master 20 19 use infotrac, only: infotrac_init, nqtot, tname 21 USE tracer_h, ONLY: igcm_n2 20 USE tracer_h, ONLY: igcm_n2,igcm_ch4_ice,igcm_co_ice,igcm_haze, 21 & igcm_prec_haze,igcm_co_gas,igcm_ch4_gas,noms 22 22 USE comsoil_h, ONLY: nsoilmx, layer, mlayer, inertiedat 23 23 USE surfdat_h, ONLY: phisfi, albedodat, 24 24 & zmea, zstd, zsig, zgam, zthe 25 ! TB24 USE nonoro_gwd_ran_mod, ONLY: du_nonoro_gwd, dv_nonoro_gwd,26 ! & east_gwstress, west_gwstress27 25 use datafile_mod, only: datadir, surfdir 28 26 use ioipsl_getin_p_mod, only: getin_p … … 30 28 use phyredem, only: physdem0, physdem1 31 29 use iostart, only: open_startphy 32 ! use slab_ice_h, only:noceanmx33 ! USE ocean_slab_mod, ONLY: nslay34 30 use filtreg_mod, only: inifilr 35 31 USE mod_const_mpi, ONLY: COMM_LMDZ … … 66 62 CHARACTER relief*3 67 63 68 69 64 c Variables pour les lectures NetCDF des fichiers "start_archive" 70 65 c-------------------------------------------------- … … 90 85 REAL w(iip1,jjp1,llm+1) 91 86 REAL pbaru(ip1jmp1,llm),pbarv(ip1jm,llm) 92 ! REAL dv(ip1jm,llm),du(ip1jmp1,llm)93 ! REAL dh(ip1jmp1,llm),dp(ip1jmp1)94 87 REAL phi(iip1,jjp1,llm) 95 88 … … 106 99 REAL tsurf(ngridmx) ! surface temperature 107 100 REAL,ALLOCATABLE :: tsoil(:,:) ! soil temperature 108 ! REAL n2ice(ngridmx) ! N2 ice layer109 101 REAL emis(ngridmx) ! surface emissivity 110 102 real emisread ! added by RW 111 103 REAL,ALLOCATABLE :: qsurf(:,:) 104 REAL,ALLOCATABLE :: qsurf_tmp(:,:) 112 105 REAL q2(ngridmx,llm+1) 113 ! REAL rnaturfi(ngridmx)114 106 real alb(iip1,jjp1),albfi(ngridmx) ! albedos 115 107 real,ALLOCATABLE :: ith(:,:,:),ithfi(:,:) ! thermal inertia (3D) 116 108 real surfith(iip1,jjp1),surfithfi(ngridmx) ! surface thermal inertia (2D) 117 ! REALlatfi(ngridmx),lonfi(ngridmx),airefi(ngridmx)118 119 INTEGER i,j,l,isoil,ig,idum 109 REAL :: latfi(ngridmx),lonfi(ngridmx),airefi(ngridmx) 110 111 INTEGER i,j,l,isoil,ig,idum,k 120 112 real mugaz ! molar mass of the atmosphere 121 113 122 integer ierr ,iref114 integer ierr 123 115 124 116 c Variables on the new grid along scalar points 125 117 c------------------------------------------------------ 126 ! REAL p(iip1,jjp1)127 118 REAL t(iip1,jjp1,llm) 128 119 REAL tset(iip1,jjp1,llm) 129 120 real phisold_newgrid(iip1,jjp1) 121 real phisold(iip1,jjp1) 130 122 REAL :: teta(iip1, jjp1, llm) 131 123 REAL :: pk(iip1,jjp1,llm) … … 135 127 REAL :: xpn,xps,xppn(iim),xpps(iim) 136 128 REAL :: p3d(iip1, jjp1, llm+1) 137 ! REAL dteta(ip1jmp1,llm)138 129 139 130 c Variable de l'ancienne grille … … 145 136 c variables diverses 146 137 c------------------- 147 real choix_1,pp 138 real choix_1,pp,choice 148 139 character*80 fichnom 149 140 character*250 filestring … … 153 144 real tsud,albsud,alb_bb,ith_bb,Tiso,Tabove 154 145 real ptoto,pcap,patm,airetot,ptotn,patmn,psea 155 ! real ssum156 146 character*1 yes 157 147 logical :: flagtset=.false. , flagps0=.false. 158 real val, val2, val3, val4 ! to store temporary variables 148 real val, val1, val2, val3, val4, dist, ref ! to store temporary variables 149 real val5, val6, val7, val8, val9, val10,val11, val12 ! to store temporary variables 159 150 real :: iceith=2000 ! thermal inertia of subterranean ice 151 integer :: iref,jref,iref1,jref1,iref2,jref2 152 integer :: igref,igref1,igref2,igref3 160 153 161 154 INTEGER :: itau … … 174 167 real fact2 175 168 169 c Special Pluto Map from Lellouch & Grundy 170 c------------------------------------------ 171 integer,parameter :: im_plu=360 ! from topo 172 integer,parameter :: jm_plu=180 173 real latv_plu(jm_plu+1),lonu_plu(im_plu+1) 174 real map_pluto_dat(im_plu,jm_plu+1) 175 176 real N2_pluto_dat(im_plu,jm_plu+1) 177 real CH4_pluto_dat(im_plu,jm_plu+1) 178 real CO_pluto_dat(im_plu,jm_plu+1) 179 real alb_dat(im_plu,jm_plu+1) 180 181 real N2_pluto_rein(im_plu+1,jm_plu+1) 182 real CH4_pluto_rein(im_plu+1,jm_plu+1) 183 real CO_pluto_rein(im_plu+1,jm_plu+1) 184 real alb_rein(im_plu+1,jm_plu+1) 185 real N2_pluto_gcm(iip1,jjp1) 186 real CH4_pluto_gcm(iip1,jjp1) 187 real CO_pluto_gcm(iip1,jjp1) 188 real alb_gcm(iip1,jjp1) 189 190 c Special Topo Map mountain 191 real lat0, lon0 192 integer ig0 193 real dist_pluto 194 real radm,height, phisprev, temp 195 real preffnew,panew 196 c Special copy of terrain 197 real actualmin,angle 198 integer array_ind(ngridmx) 199 real array_dist(ngridmx) 200 real array_angle(ngridmx) 176 201 177 202 c sortie visu pour les champs dynamiques 178 c---------------------------------------179 ! INTEGER :: visuid180 ! real :: time_step,t_ops,t_wrt181 ! CHARACTER*80 :: visu_file182 !TB24183 203 CALL conf_gcm( 99, .TRUE. ) 184 204 185 186 205 cpp = 0. 187 preff = 0.188 pa = 0. ! to ensure disaster rather than minor error if we don`t206 preff = 2. 207 pa = 0.2 ! to ensure disaster rather than minor error if we don`t 189 208 ! make deliberate choice of these values elsewhere. 190 209 … … 204 223 allocate(q(iip1,jjp1,llm,nqtot)) 205 224 allocate(qsurf(ngridmx,nqtot)) 225 allocate(qsurf_tmp(ngridmx,nqtot)) 206 226 207 227 ! get value of nsoilmx and allocate corresponding arrays … … 272 292 273 293 endif 274 275 294 276 295 c======================================================================= … … 354 373 . nqtot,day_ini,time, 355 374 . tsurf,tsoil,emis,q2,qsurf) !) ! temporary modif by RDW 356 ! . cloudfrac,totalfrac,hice,rnat,pctsrf_sic,tslab,tsea_ice,357 ! . sea_ice)358 375 359 376 ! copy albedo and soil thermal inertia on (local) physics grid … … 432 449 endif 433 450 451 c Initialisation coordonnees /aires 452 c ------------------------------- 453 ! Note: rlatu(:) and rlonv(:) are commons defined in "comgeom.h" 454 ! rlatu() and rlonv() are given in radians 455 latfi(1)=rlatu(1) 456 lonfi(1)=0. 457 DO j=2,jjm 458 DO i=1,iim 459 latfi((j-2)*iim+1+i)=rlatu(j) 460 lonfi((j-2)*iim+1+i)=rlonv(i) 461 ENDDO 462 ENDDO 463 latfi(ngridmx)=rlatu(jjp1) 464 lonfi(ngridmx)=0. 465 466 CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,aire,airefi) 467 434 468 rad = p_rad 435 469 omeg = p_omeg … … 438 472 mugaz = p_mugaz 439 473 daysec = p_daysec 440 441 if (p_omeg .eq. -9999.) then442 p_omeg = 8.*atan(1.)/p_daysec443 print*,"new value of omega",p_omeg444 endif445 446 474 kappa = 8.314*1000./(p_mugaz*p_cpp) ! added by RDW 447 475 … … 536 564 endif ! of if (choix_1.eq.0) 537 565 538 539 566 c======================================================================= 540 567 c Lecture des fichiers (start ou start_archive) … … 549 576 & surfith,nid) 550 577 551 ! & rnat,pctsrf_sic,tslab,tsea_ice,sea_ice)552 ! TB24 & du_nonoro_gwd,dv_nonoro_gwd,east_gwstress,west_gwstress)553 578 write(*,*) "OK, read start_archive file" 554 579 ! copy soil thermal inertia … … 572 597 do ! infinite loop on list of changes 573 598 574 write(*,*) 575 write(*,*) 576 write(*,*) 'List of possible changes :' 577 write(*,*) '~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~' 578 write(*,*) 579 write(*,*) 'flat : no topography ("aquaplanet")' 580 write(*,*) 'set_ps_to_preff : used if changing preff with topo' 581 write(*,*) 'qname : change tracer name' 582 write(*,*) 't=profile : read temperature profile in profile.in' 583 write(*,*) 'q=0 : ALL tracer =zero' 584 write(*,*) 'q=x : give a specific uniform value to one tracer' 585 write(*,*) 'qs=x : give a uniform value to a surface tracer' 586 write(*,*) 'q=profile : specify a profile for a tracer' 587 write(*,*) 'subsoil_all : set seasonal subsurface thermal inertia' 588 write(*,*) 'diurnal_TI : set diurnal subsurface thermal inertia' 599 write(*,*) 600 write(*,*) 601 write(*,*) 'List of possible changes :' 602 write(*,*) '~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~' 603 write(*,*) 604 write(*,*) 'flat : no topography ("aquaplanet")' 605 write(*,*) 'set_ps_to_preff : used if changing preff with topo' 606 write(*,*) 'qname : change tracer name' 607 write(*,*) 't=profile : read temperature profile in profile.in' 608 write(*,*) 'q=0 : ALL tracer =zero' 609 write(*,*) 'q=x : give a specific uniform value to one tracer' 610 write(*,*) 'qs=x : give a uniform value to a surface tracer' 611 write(*,*) 'q=profile : specify a profile for a tracer' 612 write(*,*) 'qnogcm: initial tracer for GCM from nogcm (CH4,CO)' 613 write(*,*) 'isotherm : Isothermal Temperatures' 614 write(*,*) 'tempprof : specific Temperature profile' 615 write(*,*) 'initsoil : initialize soil Temperatures' 616 write(*,*) 'ptot : change total pressure' 617 write(*,*) 'therm_ini_s: set soil TI to reference surf. values' 618 write(*,*) 'inert3d: give uniform value of thermal inertia' 619 write(*,*) 'subsoilice_n: put deep underground ice in N. hemis' 620 write(*,*) 'subsoilice_s: put deep underground ice in S. hemis' 621 write(*,*) 'reservoir: increase/decrease reservoir of ice' 622 write(*,*) 'reservoir_SP: increase/decrease reservoir in SP' 623 write(*,*) 'plutomap: initialize surface like pluto from ...' 624 write(*,*) 'subsoil_n2: diu-sea thermal inertia for N2 only' 625 write(*,*) 'subsoil_ch4: diu-sea thermal inertia for CH4 only' 626 write(*,*) 'subsoil_all: diu-sea thermal inertia for all terr' 627 write(*,*) 'subsoil_alb: diu-sea thermal inertia from albedo' 628 write(*,*) 'diurnal_TI: diurnal thermal inertia for all terr' 629 write(*,*) 'initsurf: initial surface temperature (global)' 630 write(*,*) 'modtsurf: initial surface temperature (global)' 631 write(*,*) 'copylat: copy tsurf and tsoil latitude' 632 write(*,*) 'copylatlon: copy tsurf and tsoil lat/lon' 633 write(*,*) 'copylon: copy tsurf tsoil lat, n2surf and phisfi' 634 write(*,*) 'tsurfice: temperature depending on N2 ice' 635 write(*,*) 'icarus: option to change surface/soil temperature' 636 write(*,*) 'icarus_ch4: special option to add ch4' 637 write(*,*) 'delfrostch4: delete ch4 frost' 638 write(*,*) 'modch4: remove/modify amount ch4 frost' 639 write(*,*) 'modch4n2: modify amount ch4 frost according N2' 640 write(*,*) 'modco: remove/modify amount co frost' 641 write(*,*) 'modn2: remove/modify amount n2 ice' 642 write(*,*) 'modcoch4: remove/modify co ch4 where no n2 ' 643 write(*,*) 'checktsoil: change tsoil where no n2 ' 644 write(*,*) 'non2noco: if no n2, no co ' 645 write(*,*) 'modn2_topo: modify n2 ice topo according to topo' 646 write(*,*) 'modwheren2: modify n2 where already n2' 647 write(*,*) 'modn2HD: modify n2 for HD runs' 648 write(*,*) 'modn2HD_SP: modify n2 for HD runs in SP' 649 write(*,*) 'globch4co: add/remove global amount ch4co frost' 650 write(*,*) 'source_ch4: add source ch4' 651 write(*,*) 'nomountain: remove pluto mountains (!)' 652 write(*,*) 'relief: add pluto crater or mountain' 653 write(*,*) 'seticeNH: set ice for initial start with NH topo' 654 write(*,*) 'topo_sp: change topography of SP' 655 write(*,*) 'fill_sp: fill sp with N2 ice and adjust phisfi' 656 write(*,*) 'fill_sp_inv: fill inverted sp and adjust' 657 write(*,*) 'subsoil_SP: diu-sea thermal inertia for SP ' 658 write(*,*) 'bladed: add ch4 on bladed terrains' 659 write(*,*) 'cpbladed: add extension bladed terrains' 660 write(*,*) 'slope: add slope over all long. (for triton)' 661 write(*,*) 'digsp: specific to dig SP' 662 write(*,*) 'bugn2: correct bug of warm n2 due to HighRes' 663 write(*,*) 'bugch4: correct bug of warm ch4 due to HighRes' 664 write(*,*) 'flatnosp : topo flat except SP (topo controled)' 665 write(*,*) 'flatregion: topo flat for specific region' 666 write(*,*) 'changetopo: change topo' 667 write(*,*) 'asymtopo: N-S asym topo in tanh' 668 write(*,*) 'corrtopo: correction topo bug' 669 write(*,*) 'adjustphi: adjust phisfi according to N2 mass' 670 write(*,*) 'correctphi: adjust phisfi' 671 write(*,*) 'correctps: test to correct ps' 672 write(*,*) 'toponoise: no flat topo' 673 write(*,*) 'asymtriton: asymetry in longitude for triton' 674 write(*,*) 'copytsoil: copy 2D tsoil field from nc file' 675 write(*,*) 'albedomap: read in an albedomap albedo.nc' 676 write(*,*) 'mod_distrib_ice: modify ice distrib. from albedo' 677 678 write(*,*) 679 write(*,*) 'Please note that emis and albedo are set in physiq' 680 write(*,*) 589 681 590 682 write(*,*) … … 825 917 endif 826 918 827 828 829 c subsoil_all : initialize subsurface thermal inertia 830 c -------------------------------------------------- 831 else if (modif(1:len_trim(modif)) .eq. 'subsoil_all') then 832 833 write(*,*) 'New value for subsoil thermal inertia ?' 834 104 read(*,*,iostat=ierr) ith_bb 835 if(ierr.ne.0) goto 104 836 write(*,*) 'thermal inertia (new value):',ith_bb 837 838 write(*,*)'At which depth (in m.) does the ice layer begin?' 839 write(*,*)'(here, the deepest soil layer extends down to:' 840 & ,layer(1),' - ',layer(nsoilmx),')' 841 write(*,*)'write 0 for uniform value for all subsurf levels?' 842 ierr=1 843 do while (ierr.ne.0) 844 read(*,*,iostat=ierr) val2 845 write(*,*)'val2 in subsoil_all:',val2,'ierr=',ierr 846 if (ierr.eq.0) then ! got a value, but do a sanity check 847 if(val2.gt.layer(nsoilmx)) then 848 write(*,*)'Depth should be less than ',layer(nsoilmx) 919 c qnogcm : initialise tracer from nogcm (CH4, CO) 920 c -------------------------------- 921 else if (modif(1:len_trim(modif)).eq.'qnogcm') then 922 write(*,*) 'Which tracer do you want to modify ?' 923 write(*,*) 'Should be ch4_gas and co_gas' 924 do iq=1,nqtot 925 write(*,*)iq,' : ',trim(noms(iq)),' : ',q(1,1,1,iq) 926 enddo 927 write(*,*) '(choose between 1 and ',nqtot,')' 928 read(*,*) iq 929 DO l=1,llm 930 DO j=1,jjp1 931 DO i=1,iip1 932 q(i,j,l,iq)=q(i,j,1,iq) 933 ENDDO 934 ENDDO 935 ENDDO 936 937 c isothermal temperatures 938 c ------------------------------------------------ 939 else if (modif(1:len_trim(modif)) .eq. 'isotherm') then 940 941 write(*,*)'Isothermal temperature of the atmosphere' 942 write(*,*) 'Value of THIS temperature ? :' 943 203 read(*,*,iostat=ierr) Tset(1,1,1) 944 if(ierr.ne.0) goto 203 945 flagtset=.true. 946 DO l=1,llm 947 DO j=1,jjp1 948 DO i=1,iip1 949 Tset(i,j,l)=Tset(1,1,1) 950 ENDDO 951 ENDDO 952 ENDDO 953 write(*,*) 'atmospheric temp =', Tset(2,2,2) 954 955 c specific temperature profile 956 c ------------------------------------------------ 957 else if (modif(1:len_trim(modif)) .eq. 'tempprof') then 958 959 write(*,*) 'phi=' 960 write(*,*) phi(1,1,:) 961 write(*,*) 'temperature profile in the atmosphere' 962 write(*,*) 'Value of THIS temperature ? :' 963 204 read(*,*,iostat=ierr) Tset(1,1,1) 964 if(ierr.ne.0) goto 204 965 flagtset=.true. 966 DO l=1,llm 967 DO j=1,jjp1 968 DO i=1,iip1 969 Tset(i,j,l)=Tset(1,1,1) 970 ENDDO 971 ENDDO 972 ENDDO 973 write(*,*) 'atmospheric temp =', Tset(2,2,2) 974 975 c initsoil: subsurface temperature 976 c --------------------------------- 977 else if (modif(1:len_trim(modif)) .eq. 'initsoil') then 978 979 write(*,*)'Isothermal temperature of the subsurface' 980 write(*,*) 'Value of this temperature ? :' 981 303 read(*,*,iostat=ierr) Tiso 982 if(ierr.ne.0) goto 303 983 984 do l=1,nsoilmx 985 do ig=1, ngridmx 986 tsoil(ig,l) = Tiso 987 end do 988 end do 989 990 c icarus : changing tsoil tsurf on latitudinal bands 991 c ----------------------------------------------------- 992 else if (modif(1:len_trim(modif)) .eq. 'icarus') then 993 994 write(*,*) 'At which lat lon do you want to extract the 995 & reference tsurf/tsoil profile ?' 996 407 read(*,*,iostat=ierr) val,val2 997 if(ierr.ne.0) goto 407 998 write(*,*) 'You want lat =',val 999 write(*,*) 'You want lon =',val2 1000 dist=0 1001 ref=1000 1002 do j=1,ngridmx-1 1003 dist=sqrt((lonfi(j)*180./pi-val2)**2+ 1004 & (latfi(j)*180./pi-val)**2) 1005 IF (dist.lt.ref) then 1006 ref=dist 1007 jref=j 1008 ENDIF 1009 ENDDO 1010 1011 write(*,*)'Will use nearest grid point which is:', 1012 & latfi(jref)*180./pi,lonfi(jref)*180./pi 1013 1014 ! Extraction of the profile 1015 write(*,*) 'Profile is =',tsoil(jref,:) 1016 write(*,*) 'tsurf is =',tsurf(jref) 1017 write(*,*) 'Choice lat for latitudinal band with this profile' 1018 408 read(*,*,iostat=ierr) val3 1019 if(ierr.ne.0) goto 408 1020 write(*,*) 'Choice delta lat for this latitudinal band' 1021 409 read(*,*,iostat=ierr) val4 1022 if(ierr.ne.0) goto 409 1023 write(*,*) 'Choice delta temp (K) for this latitudinal band' 1024 410 read(*,*,iostat=ierr) val5 1025 if(ierr.ne.0) goto 410 1026 write(*,*) 'How much N2 ice should I put on this band (kg/m2)' 1027 415 read(*,*,iostat=ierr) choice 1028 if(ierr.ne.0) goto 415 1029 write(*,*) 'Choice lat2' 1030 411 read(*,*,iostat=ierr) val6 1031 if(ierr.ne.0) goto 411 1032 write(*,*) 'Choice delta lat for this latitudinal band' 1033 412 read(*,*,iostat=ierr) val7 1034 if(ierr.ne.0) goto 412 1035 write(*,*) 'Choice delta temp (K) for this latitudinal band' 1036 413 read(*,*,iostat=ierr) val8 1037 if(ierr.ne.0) goto 413 1038 1039 DO ig=1,ngridmx 1040 IF ( (latfi(ig)*180./pi.ge.(val3-val4/2.)) .and. 1041 & (latfi(ig)*180./pi.le.(val3+val4/2.)) .and. 1042 & (qsurf(ig,igcm_n2).lt.0.001) ) then 1043 tsurf(ig)=tsurf(jref)+val5 1044 DO l=1,nsoilmx 1045 tsoil(ig,l)=tsoil(jref,l)+val5 1046 ENDDO 1047 qsurf(ig,igcm_n2)=choice 1048 ENDIF 1049 IF ( (latfi(ig)*180./pi.ge.(val6-val7/2.)) .and. 1050 & (latfi(ig)*180./pi.le.(val6+val7/2.)) .and. 1051 & (qsurf(ig,igcm_n2).lt.0.001) ) then 1052 tsurf(ig)=tsurf(jref)+val8 1053 DO l=1,nsoilmx 1054 tsoil(ig,l)=tsoil(jref,l)+val8 1055 ENDDO 1056 ENDIF 1057 ENDDO 1058 1059 c icarus_ch4 : changing tsoil tsurf and adding ch4 1060 c ----------------------------------------------------- 1061 else if (modif(1:len_trim(modif)) .eq. 'icarus_ch4') then 1062 1063 write(*,*) 'At which lat lon do you want to extract the 1064 & reference tsurf/tsoil profile ?' 1065 416 read(*,*,iostat=ierr) val,val2 1066 if(ierr.ne.0) goto 416 1067 write(*,*) 'You want lat =',val 1068 write(*,*) 'You want lon =',val2 1069 dist=0 1070 ref=1000 1071 do j=1,ngridmx-1 1072 dist=sqrt((lonfi(j)*180./pi-val2)**2+ 1073 & (latfi(j)*180./pi-val)**2) 1074 IF (dist.lt.ref) then 1075 ref=dist 1076 jref=j 1077 ENDIF 1078 ENDDO 1079 1080 write(*,*)'Will use nearest grid point which is:', 1081 & latfi(jref)*180./pi,lonfi(jref)*180./pi 1082 1083 ! Extraction of the profile 1084 write(*,*) 'Profile is =',tsoil(jref,:) 1085 write(*,*) 'tsurf is =',tsurf(jref) 1086 write(*,*) 'Choice band : lat1 and lat2 ?' 1087 417 read(*,*,iostat=ierr) val3,val4 1088 if(ierr.ne.0) goto 417 1089 write(*,*) 'Choice temp coefficient for this latitudinal band' 1090 418 read(*,*,iostat=ierr) val5 1091 if(ierr.ne.0) goto 418 1092 write(*,*) 'How much CH4 ice do I put on this band (kg/m2)' 1093 419 read(*,*,iostat=ierr) choice 1094 if(ierr.ne.0) goto 419 1095 1096 DO ig=1,ngridmx 1097 IF ( (latfi(ig)*180./pi.ge.val3) .and. 1098 & (latfi(ig)*180./pi.le.val4) .and. 1099 & (qsurf(ig,igcm_ch4_ice).lt.0.001) ) then 1100 tsurf(ig)=tsurf(jref)*val5 1101 DO l=1,nsoilmx 1102 tsoil(ig,l)=tsoil(jref,l)*val5 1103 ENDDO 1104 qsurf(ig,igcm_ch4_ice)=choice 1105 ENDIF 1106 ENDDO 1107 1108 c globch4co : adding/remove global ch4 co 1109 c ---------------------------------------------------------- 1110 else if (modif(1:len_trim(modif)) .eq. 'globch4co') then 1111 1112 write(*,*) 'Adding or removing ch4 co ' 1113 write(*,*) 'Choice amount add/less ch4 ice (0 = rm all)' 1114 431 read(*,*,iostat=ierr) val3 1115 if(ierr.ne.0) goto 431 1116 write(*,*) 'Choice amount add/less co ice (0 = rm all)' 1117 432 read(*,*,iostat=ierr) val4 1118 if(ierr.ne.0) goto 432 1119 IF (val3.eq.0.) then 1120 DO ig=1,ngridmx 1121 qsurf(ig,igcm_ch4_ice)=0. 1122 ENDDO 1123 ELSE 1124 DO ig=1,ngridmx 1125 qsurf(ig,igcm_ch4_ice)=qsurf(ig,igcm_ch4_ice)+val3 1126 ENDDO 1127 ENDIF 1128 IF (val4.eq.0.) then 1129 DO ig=1,ngridmx 1130 qsurf(ig,igcm_co_ice)=0. 1131 ENDDO 1132 ELSE 1133 DO ig=1,ngridmx 1134 qsurf(ig,igcm_co_ice)=qsurf(ig,igcm_co_ice)+val4 1135 ENDDO 1136 ENDIF 1137 1138 c bladed : adding/remove ch4 on bladed terrains 1139 c ---------------------------------------------------------- 1140 else if (modif(1:len_trim(modif)).eq.'bladed') then 1141 1142 write(*,*) 'Adding or removing ch4 on the bladed terrains' 1143 write(*,*) 'Choice band : lat1 and lat2 ?' 1144 450 read(*,*,iostat=ierr) val,val2 1145 if(ierr.ne.0) goto 450 1146 write(*,*) 'Choice band : lon1 and lon2 ?' 1147 451 read(*,*,iostat=ierr) val3,val4 1148 if(ierr.ne.0) goto 451 1149 write(*,*) 'Choice phisfi minimum ?' 1150 454 read(*,*,iostat=ierr) choice 1151 if(ierr.ne.0) goto 454 1152 write(*,*) 'Choice multiplicative factor' 1153 452 read(*,*,iostat=ierr) val5 1154 if(ierr.ne.0) goto 452 1155 write(*,*) 'Choice additional ch4' 1156 453 read(*,*,iostat=ierr) val6 1157 if(ierr.ne.0) goto 453 1158 write(*,*) 'Choice index for tsurf tsoil' 1159 449 read(*,*,iostat=ierr) iref 1160 if(ierr.ne.0) goto 449 1161 1162 ! shift 1163 DO ig=1,ngridmx 1164 IF ( (latfi(ig)*180./pi.ge.val) .and. 1165 & (latfi(ig)*180./pi.le.val2) .and. 1166 & (lonfi(ig)*180./pi.ge.val3) .and. 1167 & (lonfi(ig)*180./pi.le.val4) .and. 1168 & (phisfi(ig).gt.choice) ) then 1169 qsurf(ig,igcm_ch4_ice)=qsurf(ig,igcm_ch4_ice)*val5 1170 qsurf(ig,igcm_ch4_ice)=qsurf(ig,igcm_ch4_ice)+val6 1171 tsurf(ig)=tsurf(iref) 1172 DO l=1,nsoilmx 1173 tsoil(ig,l)=tsoil(iref,l) 1174 ENDDO 1175 ENDIF 1176 ENDDO 1177 1178 c cpbladed : add extension bladed terrains 1179 c ---------------------------------------------------------- 1180 else if (modif(1:len_trim(modif)).eq.'cpbladed') then 1181 1182 write(*,*) 'Choice lat,lon, centre of band to be copied ?' 1183 560 read(*,*,iostat=ierr) val,val2 1184 if(ierr.ne.0) goto 560 1185 write(*,*) 'Choice distance (km) from there defining the area' 1186 561 read(*,*,iostat=ierr) val3 1187 if(ierr.ne.0) goto 561 1188 write(*,*) 'Nb of copies ?' 1189 562 read(*,*,iostat=ierr) l 1190 if(ierr.ne.0) goto 562 1191 1192 ! Get index correponding to central points 1193 ! 1/ Reference point 1194 igref=1 1195 actualmin=1000. 1196 do ig=1,ngridmx 1197 val6=dist_pluto(latfi(ig),lonfi(ig), 1198 & val*pi/180.,val2*pi/180.) 1199 if (val6.lt.actualmin) then 1200 actualmin=val6 1201 igref=ig 1202 endif 1203 enddo 1204 1205 do k=1,l 1206 1207 write(*,*) 'Choice lat,lon where terrain is copied' 1208 563 read(*,*,iostat=ierr) val4,val5 1209 if(ierr.ne.0) goto 563 1210 1211 if (val5.gt.180.) then 1212 val5=val5-360. 1213 endif 1214 1215 ! 2/ Target point 1216 igref2=1 1217 actualmin=1000. 1218 do ig=1,ngridmx 1219 val6=dist_pluto(latfi(ig),lonfi(ig), 1220 & val4*pi/180.,val5*pi/180.) 1221 if (val6.lt.actualmin) then 1222 actualmin=val6 1223 igref2=ig 1224 endif 1225 enddo 1226 1227 ! for each point within A1, get distance and angle 1228 ! save in a table 1229 i=1 1230 array_ind(:)=0 1231 array_dist(:)=1000. 1232 array_angle(:)=0. 1233 do ig=1,ngridmx 1234 val6=dist_pluto(latfi(ig),lonfi(ig), 1235 & val*pi/180.,val2*pi/180.) 1236 if (val6.lt.val3) then 1237 array_ind(i)=ig 1238 array_dist(i)=val6 1239 array_angle(i)=atan2(val-latfi(ig)*180./pi, 1240 & val2-lonfi(ig)*180./pi) 1241 i=i+1 1242 endif 1243 enddo 1244 1245 ! for each point within A2, get distance and angle 1246 ! and look where fit with previous table, and set value 1247 1248 do ig=1,ngridmx 1249 val6=dist_pluto(latfi(ig),lonfi(ig), 1250 & val4*pi/180.,val5*pi/180.) 1251 if (val6.lt.val3) then 1252 ! dist = val6 1253 ! get angle: 1254 angle=atan2(val4-latfi(ig)*180./pi, 1255 & val5-lonfi(ig)*180./pi) 1256 ! Loop where min 1257 actualmin=1000. 1258 do j=1,i 1259 if ( (array_angle(j).lt.angle+0.52).and. 1260 & (array_angle(j).gt.angle-0.52).and. 1261 & (array_dist(j)-val6.lt.actualmin) ) then 1262 actualmin=array_dist(j)-val6 1263 igref3=j 1264 endif 1265 enddo 1266 phisfi(ig)=phisfi(array_ind(igref3)) 1267 qsurf(ig,igcm_ch4_ice)= 1268 & qsurf(array_ind(igref3),igcm_ch4_ice) 1269 tsurf(ig)=tsurf(array_ind(igref3)) 1270 endif 1271 enddo 1272 enddo 1273 CALL gr_fi_dyn(1,ngridmx,iip1,jjp1,phisfi,phis) 1274 1275 c delfrostch4/ delete ch4 frost on a latitudinal band 1276 c ---------------------------------------------------------- 1277 else if (modif(1:len_trim(modif)) .eq. 'delfrostch4') then 1278 1279 write(*,*) 'removing ch4 latitudinally' 1280 write(*,*) 'Choice band : lat1 and lat2 ?' 1281 read(*,*,iostat=ierr) val,val2 1282 write(*,*) 'Choice band : lon1 and lon2 ?' 1283 522 read(*,*,iostat=ierr) val4,val5 1284 if(ierr.ne.0) goto 522 1285 write(*,*) 'Choice amount max below whcih ch4 is removed' 1286 521 read(*,*,iostat=ierr) val3 1287 if(ierr.ne.0) goto 521 1288 1289 DO ig=1,ngridmx 1290 IF ( (latfi(ig)*180./pi.ge.val) .and. 1291 & (latfi(ig)*180./pi.le.val2) .and. 1292 & (lonfi(ig)*180./pi.ge.val4) .and. 1293 & (lonfi(ig)*180./pi.lt.val5) .and. 1294 & (qsurf(ig,igcm_ch4_ice).lt.val3) ) then 1295 qsurf(ig,igcm_ch4_ice)=0. 1296 ENDIF 1297 ENDDO 1298 1299 c modch4 : adding/remove ch4 frost on a latitudinal band 1300 c ---------------------------------------------------------- 1301 else if (modif(1:len_trim(modif)) .eq. 'modch4') then 1302 1303 write(*,*) 'Adding or removing ch4 latitudinally' 1304 write(*,*) 'Choice band : lat1 and lat2 ?' 1305 read(*,*,iostat=ierr) val,val2 1306 write(*,*) 'Choice band : lon1 and lon2 ?' 1307 422 read(*,*,iostat=ierr) val4,val5 1308 if(ierr.ne.0) goto 422 1309 write(*,*) 'Choice multiplicative factor' 1310 421 read(*,*,iostat=ierr) val3 1311 if(ierr.ne.0) goto 421 1312 write(*,*) 'Choice additional ch4' 1313 423 read(*,*,iostat=ierr) val6 1314 if(ierr.ne.0) goto 423 1315 1316 DO ig=1,ngridmx 1317 IF ( (latfi(ig)*180./pi.ge.val) .and. 1318 & (latfi(ig)*180./pi.le.val2) .and. 1319 & (lonfi(ig)*180./pi.ge.val4) .and. 1320 & (lonfi(ig)*180./pi.lt.val5)) then 1321 ! & (qsurf(ig,igcm_n2).gt.50)) then 1322 ! & (qsurf(ig,igcm_ch4_ice).lt.10) ) then 1323 qsurf(ig,igcm_ch4_ice)=qsurf(ig,igcm_ch4_ice)*val3 1324 qsurf(ig,igcm_ch4_ice)=qsurf(ig,igcm_ch4_ice)+val6 1325 ENDIF 1326 ENDDO 1327 1328 c modch4n2 : adding/remove ch4 frost on a latitudinal band 1329 c ---------------------------------------------------------- 1330 else if (modif(1:len_trim(modif)) .eq. 'modch4n2') then 1331 1332 write(*,*) 'Adding or removing ch4 latitudinally' 1333 write(*,*) 'Choice band : lat1 and lat2 ?' 1334 read(*,*,iostat=ierr) val,val2 1335 write(*,*) 'Choice band : lon1 and lon2 ?' 1336 read(*,*,iostat=ierr) val4,val5 1337 write(*,*) 'Choice range n2 for modif' 1338 read(*,*,iostat=ierr) val7,val8 1339 write(*,*) 'Choice multiplicative factor ch4' 1340 read(*,*,iostat=ierr) val3 1341 write(*,*) 'Choice additional ch4' 1342 read(*,*,iostat=ierr) val6 1343 1344 DO ig=1,ngridmx 1345 IF ( (latfi(ig)*180./pi.ge.val) .and. 1346 & (latfi(ig)*180./pi.le.val2) .and. 1347 & (lonfi(ig)*180./pi.ge.val4) .and. 1348 & (lonfi(ig)*180./pi.lt.val5) .and. 1349 & (qsurf(ig,igcm_n2).gt.val7) .and. 1350 & (qsurf(ig,igcm_n2).lt.val8) ) then 1351 qsurf(ig,igcm_ch4_ice)=qsurf(ig,igcm_ch4_ice)*val3 1352 qsurf(ig,igcm_ch4_ice)=qsurf(ig,igcm_ch4_ice)+val6 1353 ENDIF 1354 ENDDO 1355 1356 c non2noco : if no n2 no co 1357 c ---------------------------------------------------------- 1358 else if (modif(1:len_trim(modif)) .eq. 'non2noco') then 1359 DO ig=1,ngridmx 1360 IF ( (qsurf(ig,igcm_n2).lt.0.001)) then 1361 qsurf(ig,igcm_co_ice)=0. 1362 ENDIF 1363 ENDDO 1364 1365 c modco : adding/remove co frost on a latitudinal band 1366 c ---------------------------------------------------------- 1367 else if (modif(1:len_trim(modif)) .eq. 'modco') then 1368 1369 write(*,*) 'Adding or removing co where N2 is present ' 1370 write(*,*) 'Choice limit mini n2 pour mettre co?' 1371 460 read(*,*,iostat=ierr) val7 1372 if(ierr.ne.0) goto 460 1373 write(*,*) 'Choice band : lat1 and lat2 ?' 1374 461 read(*,*,iostat=ierr) val,val2 1375 if(ierr.ne.0) goto 461 1376 write(*,*) 'Choice band : lon1 and lon2 ?' 1377 462 read(*,*,iostat=ierr) val4,val5 1378 if(ierr.ne.0) goto 462 1379 write(*,*) 'Choice multiplicative factor' 1380 463 read(*,*,iostat=ierr) val3 1381 if(ierr.ne.0) goto 463 1382 write(*,*) 'Choice additional co' 1383 464 read(*,*,iostat=ierr) val6 1384 if(ierr.ne.0) goto 464 1385 1386 DO ig=1,ngridmx 1387 IF ( (latfi(ig)*180./pi.ge.val) .and. 1388 & (latfi(ig)*180./pi.le.val2) .and. 1389 & (lonfi(ig)*180./pi.ge.val4) .and. 1390 & (lonfi(ig)*180./pi.lt.val5) .and. 1391 & (qsurf(ig,igcm_n2).lt.val7)) then 1392 qsurf(ig,igcm_co_ice)=qsurf(ig,igcm_co_ice)*val3 1393 qsurf(ig,igcm_co_ice)=qsurf(ig,igcm_co_ice)+val6 1394 ENDIF 1395 ENDDO 1396 1397 c modn2 : adding/remove n2 ice 1398 c ---------------------------------------------------------- 1399 else if (modif(1:len_trim(modif)) .eq. 'modn2') then 1400 1401 write(*,*) 'Adding or removing n2 ' 1402 write(*,*) 'Choice band : lat1 and lat2 ?' 1403 425 read(*,*,iostat=ierr) val,val2 1404 if(ierr.ne.0) goto 425 1405 write(*,*) 'Choice band : lon1 and lon2 ?' 1406 426 read(*,*,iostat=ierr) val4,val5 1407 if(ierr.ne.0) goto 426 1408 write(*,*) 'Choice multiplicative factor' 1409 427 read(*,*,iostat=ierr) val3 1410 if(ierr.ne.0) goto 427 1411 write(*,*) 'Choice additional n2' 1412 428 read(*,*,iostat=ierr) val6 1413 if(ierr.ne.0) goto 428 1414 1415 DO ig=1,ngridmx 1416 IF ( (latfi(ig)*180./pi.ge.val) .and. 1417 & (latfi(ig)*180./pi.le.val2) .and. 1418 & (lonfi(ig)*180./pi.ge.val4) .and. 1419 & (lonfi(ig)*180./pi.lt.val5) ) then 1420 c & (qsurf(ig,igcm_n2).lt.50)) then 1421 qsurf(ig,igcm_n2)=qsurf(ig,igcm_n2)*val3 1422 qsurf(ig,igcm_n2)=qsurf(ig,igcm_n2)+val6 1423 ENDIF 1424 ! IF ((phisfi(ig).gt.-1000.)) then 1425 ! qsurf(ig,igcm_n2)=0. 1426 ! ENDIF 1427 ENDDO 1428 1429 c modcoch4 : adding/remove co and ch4 frost where co without n2 1430 c ---------------------------------------------------------- 1431 else if (modif(1:len_trim(modif)) .eq. 'modcoch4') then 1432 1433 write(*,*) 'Adding or removing co where N2 is not present ' 1434 write(*,*) 'Choice band : lat1 and lat2 ?' 1435 491 read(*,*,iostat=ierr) val,val2 1436 if(ierr.ne.0) goto 491 1437 write(*,*) 'Choice band : lon1 and lon2 ?' 1438 492 read(*,*,iostat=ierr) val4,val5 1439 if(ierr.ne.0) goto 492 1440 write(*,*) 'Choice multiplicative factor' 1441 493 read(*,*,iostat=ierr) val3 1442 if(ierr.ne.0) goto 493 1443 write(*,*) 'Choice additional co ch4' 1444 494 read(*,*,iostat=ierr) val6 1445 if(ierr.ne.0) goto 494 1446 1447 DO ig=1,ngridmx 1448 IF ( (latfi(ig)*180./pi.ge.val) .and. 1449 & (latfi(ig)*180./pi.le.val2) .and. 1450 & (lonfi(ig)*180./pi.ge.val4) .and. 1451 & (lonfi(ig)*180./pi.le.val5) .and. 1452 & (qsurf(ig,igcm_n2).lt.10.)) then 1453 qsurf(ig,igcm_co_ice)=qsurf(ig,igcm_co_ice)*val3 1454 qsurf(ig,igcm_co_ice)=qsurf(ig,igcm_co_ice)+val6 1455 qsurf(ig,igcm_ch4_ice)=qsurf(ig,igcm_ch4_ice)*val3 1456 qsurf(ig,igcm_ch4_ice)=qsurf(ig,igcm_ch4_ice)+val6 1457 ENDIF 1458 ENDDO 1459 1460 c modn2HD : adding/remove n2 ice for HD runs 1461 c ---------------------------------------------------------- 1462 else if (modif(1:len_trim(modif)) .eq. 'modn2HD') then 1463 1464 write(*,*) 'Adding or removing n2 ' 1465 write(*,*) 'Choice band : lat1 and lat2 ?' 1466 480 read(*,*,iostat=ierr) val,val1 1467 if(ierr.ne.0) goto 480 1468 write(*,*) 'Choice band : lon1 and lon2 ?' 1469 481 read(*,*,iostat=ierr) val2,val3 1470 if(ierr.ne.0) goto 481 1471 write(*,*) 'Choice amount N2 min max where to modify' 1472 482 read(*,*,iostat=ierr) val4,val5 1473 if(ierr.ne.0) goto 482 1474 write(*,*) 'Choice phisfi min max where change n2' 1475 483 read(*,*,iostat=ierr) val6,val11 1476 if(ierr.ne.0) goto 483 1477 write(*,*) 'Choice multiplicative factor' 1478 484 read(*,*,iostat=ierr) val7 1479 if(ierr.ne.0) goto 484 1480 write(*,*) 'Choice additional n2' 1481 485 read(*,*,iostat=ierr) val8 1482 if(ierr.ne.0) goto 485 1483 ! write(*,*) 'Choice ind lon equivalent for change tsurf tsoil' 1484 ! 486 read(*,*,iostat=ierr) val9 1485 ! if(ierr.ne.0) goto 486 1486 1487 DO ig=1,ngridmx 1488 IF ( (latfi(ig)*180./pi.ge.val) .and. 1489 & (latfi(ig)*180./pi.le.val1) .and. 1490 & (lonfi(ig)*180./pi.ge.val2) .and. 1491 & (lonfi(ig)*180./pi.lt.val3) .and. 1492 & (qsurf(ig,igcm_n2).ge.val4) .and. 1493 & (qsurf(ig,igcm_n2).lt.val5) .and. 1494 & (phisfi(ig).ge.val6) .and. 1495 & (phisfi(ig).le.val11) ) then 1496 1497 qsurf(ig,igcm_n2)=qsurf(ig,igcm_n2)*val7 1498 qsurf(ig,igcm_n2)=qsurf(ig,igcm_n2)+val8 1499 !qsurf(ig,igcm_ch4_ice)=0. 1500 qsurf(ig,igcm_co_ice)=0. 1501 ! index of ig 1502 !if (val9.eq.0.) then 1503 ! iref=2546 1504 !else 1505 ! val10=modulo((ig-1),96) 1506 ! choice=ig+int(96/2)-val10 1507 ! !choice=int(1+96*(val9-1)+val10) 1508 !endif 1509 !tsurf(ig)=tsurf(iref) 1510 !DO l=1,nsoilmx 1511 ! tsoil(ig,l)=tsoil(iref,l) 1512 !ENDDO 1513 tsurf(ig)=tsurf(ig)+4 1514 DO l=1,nsoilmx 1515 tsoil(ig,l)=tsoil(ig,l)+4 1516 ENDDO 1517 ENDIF 1518 ! IF ((phisfi(ig).gt.-1000.)) then 1519 ! qsurf(ig,igcm_n2)=0. 1520 ! ENDIF 1521 ENDDO 1522 1523 c modn2HD_SP : adding/remove n2 ice for HD runs 1524 c ---------------------------------------------------------- 1525 else if (modif(1:len_trim(modif)) .eq. 'modn2HD_SP') then 1526 1527 write(*,*) 'Adding or removing n2 ' 1528 write(*,*) 'Choice band : lat1 and lat2 ?' 1529 501 read(*,*,iostat=ierr) val,val1 1530 if(ierr.ne.0) goto 501 1531 write(*,*) 'Choice band : lon1 and lon2 ?' 1532 502 read(*,*,iostat=ierr) val2,val3 1533 if(ierr.ne.0) goto 502 1534 write(*,*) 'Choice amount N2 min max where to modify' 1535 503 read(*,*,iostat=ierr) val4,val5 1536 if(ierr.ne.0) goto 503 1537 write(*,*) 'Choice phisfi min max where change n2' 1538 504 read(*,*,iostat=ierr) val6,val11 1539 if(ierr.ne.0) goto 504 1540 write(*,*) 'Choice multiplicative factor' 1541 505 read(*,*,iostat=ierr) val7 1542 if(ierr.ne.0) goto 505 1543 write(*,*) 'Choice additional n2' 1544 506 read(*,*,iostat=ierr) val8 1545 if(ierr.ne.0) goto 506 1546 write(*,*) 'Choice ind lon equivalent for change tsurf tsoil' 1547 507 read(*,*,iostat=ierr) iref 1548 if(ierr.ne.0) goto 507 1549 1550 DO ig=1,ngridmx 1551 IF ( (latfi(ig)*180./pi.ge.val) .and. 1552 & (latfi(ig)*180./pi.le.val1) .and. 1553 & (lonfi(ig)*180./pi.ge.val2) .and. 1554 & (lonfi(ig)*180./pi.lt.val3) .and. 1555 & (qsurf(ig,igcm_n2).ge.val4) .and. 1556 & (qsurf(ig,igcm_n2).lt.val5) .and. 1557 & (phisfi(ig).ge.val6) .and. 1558 & (phisfi(ig).le.val11) ) then 1559 1560 qsurf(ig,igcm_n2)=qsurf(ig,igcm_n2)*val7 1561 qsurf(ig,igcm_n2)=qsurf(ig,igcm_n2)+val8 1562 qsurf(ig,igcm_ch4_ice)=0. 1563 qsurf(ig,igcm_co_ice)=0. 1564 ! index of ig 1565 !if (val9.eq.0.) then 1566 ! iref=2546 1567 !else 1568 !val10=modulo((ig-1),96) 1569 !choice=ig+int(96/2)-val10 1570 !choice=int(1+96*(val9-1)+val10) 1571 if (iref.ge.1) then 1572 tsurf(ig)=tsurf(iref) 1573 DO l=1,nsoilmx 1574 tsoil(ig,l)=tsoil(iref,l) 1575 ENDDO 1576 else if (iref.eq.0) then 1577 choice=int(ig/96.)*96+84 1578 print*, 'choice=',choice 1579 tsurf(ig)=tsurf(int(choice)) 1580 DO l=1,nsoilmx 1581 tsoil(ig,l)=tsoil(int(choice),l) 1582 ENDDO 1583 endif 1584 ENDIF 1585 ENDDO 1586 1587 1588 c modn2_topo : adding/remove n2 ice 1589 c ---------------------------------------------------------- 1590 else if (modif(1:len_trim(modif)) .eq. 'modn2_topo') then 1591 1592 write(*,*) 'Adding or removing n2 ' 1593 write(*,*) 'Choice band : lat1 and lat2 ?' 1594 441 read(*,*,iostat=ierr) val,val2 1595 if(ierr.ne.0) goto 441 1596 write(*,*) 'Choice band : lon1 and lon2 ?' 1597 442 read(*,*,iostat=ierr) val4,val5 1598 if(ierr.ne.0) goto 442 1599 write(*,*) 'Choice topo 1 et 2 (phi) where change is active' 1600 443 read(*,*,iostat=ierr) val7,val8 1601 if(ierr.ne.0) goto 443 1602 write(*,*) 'Choice multiplicative factor' 1603 444 read(*,*,iostat=ierr) val3 1604 if(ierr.ne.0) goto 444 1605 write(*,*) 'Choice additional n2' 1606 445 read(*,*,iostat=ierr) val6 1607 if(ierr.ne.0) goto 445 1608 1609 DO ig=1,ngridmx 1610 IF ( (latfi(ig)*180./pi.ge.val) .and. 1611 & (latfi(ig)*180./pi.le.val2) .and. 1612 & (lonfi(ig)*180./pi.ge.val4) .and. 1613 & (lonfi(ig)*180./pi.lt.val5) .and. 1614 & (phisfi(ig).le.val8) .and. 1615 & (phisfi(ig).ge.val7) ) then 1616 qsurf(ig,igcm_n2)=qsurf(ig,igcm_n2)*val3 1617 qsurf(ig,igcm_n2)=qsurf(ig,igcm_n2)+val6 1618 ENDIF 1619 ENDDO 1620 1621 c modwheren2 : adding/remove n2 ice where already n2 1622 c ---------------------------------------------------------- 1623 else if (modif(1:len_trim(modif)) .eq. 'modwheren2') then 1624 1625 write(*,*) 'Choice multiplicative factor' 1626 471 read(*,*,iostat=ierr) val3 1627 if(ierr.ne.0) goto 471 1628 write(*,*) 'Choice additional n2' 1629 472 read(*,*,iostat=ierr) val6 1630 if(ierr.ne.0) goto 472 1631 1632 DO ig=1,ngridmx 1633 IF ( qsurf(ig,igcm_n2).gt.0.001 ) then 1634 qsurf(ig,igcm_n2)=qsurf(ig,igcm_n2)*val3 1635 qsurf(ig,igcm_n2)=qsurf(ig,igcm_n2)+val6 1636 ENDIF 1637 ENDDO 1638 1639 c checktsoil : changing tsoil if no n2 1640 c ---------------------------------------------------------- 1641 else if (modif(1:len_trim(modif)) .eq. 'checktsoil') then 1642 1643 write(*,*) 'selecting area where tsoil to be check ' 1644 write(*,*) 'Choice band : lat1 and lat2 ?' 1645 511 read(*,*,iostat=ierr) val,val1 1646 if(ierr.ne.0) goto 511 1647 write(*,*) 'Choice band : lon1 and lon2 ?' 1648 512 read(*,*,iostat=ierr) val2,val3 1649 if(ierr.ne.0) goto 512 1650 write(*,*) 'Choice temp min max in which change is made' 1651 513 read(*,*,iostat=ierr) val4,val5 1652 if(ierr.ne.0) goto 513 1653 write(*,*) 'Choice phisfi min max where change n2' 1654 514 read(*,*,iostat=ierr) val6,val11 1655 if(ierr.ne.0) goto 514 1656 1657 DO ig=1,ngridmx 1658 IF ( (latfi(ig)*180./pi.ge.val) .and. 1659 & (latfi(ig)*180./pi.le.val1) .and. 1660 & (lonfi(ig)*180./pi.ge.val2) .and. 1661 & (lonfi(ig)*180./pi.le.val3) .and. 1662 & (((tsurf(ig).ge.val4) .and. 1663 & (tsurf(ig).le.val5)) .or. 1664 & ((tsoil(ig,nsoilmx).ge.val4) .and. 1665 & (tsoil(ig,nsoilmx).le.val5))) .and. 1666 & (phisfi(ig).ge.val6) .and. 1667 & (phisfi(ig).le.val11) .and. 1668 & (qsurf(ig,igcm_n2).le.0.001) ) then 1669 1670 ! DO i=1,ngridmx 1671 ! IF ( (latfi(i).eq.latfi(ig)) .and. 1672 ! & (lonfi(i).eq.0.) ) then 1673 ! 1674 tsurf(ig)=50. 1675 !qsurf(ig,igcm_ch4_ice)=qsurf(i,igcm_ch4_ice) 1676 ! 1677 DO l=1,nsoilmx 1678 tsoil(ig,l)=50. !tsoil(i,l) 1679 ENDDO 1680 !ENDIF 1681 !ENDDO 1682 ENDIF 1683 ENDDO 1684 1685 c asymtriton : changing ice, tsurf and tsoil on triton 1686 c ---------------------------------------------------------- 1687 else if (modif(1:len_trim(modif)) .eq. 'asymtriton') then 1688 1689 write(*,*) 'selecting area where tsoil to be check ' 1690 write(*,*) 'Choice center latitude of sinuisoid distrib?' 1691 531 read(*,*,iostat=ierr) val1 1692 if(ierr.ne.0) goto 531 1693 write(*,*) 'Choice maginutde in latitude?' 1694 532 read(*,*,iostat=ierr) val2 1695 if(ierr.ne.0) goto 532 1696 write(*,*) 'Choice center longitude?' 1697 533 read(*,*,iostat=ierr) val3 1698 if(ierr.ne.0) goto 533 1699 ! write(*,*) 'iip1,jjp1=',iip1,jjp1,ngridmx 1700 1701 DO ig=1,ngridmx 1702 ! Latitude of the sinusoid: 1703 val11=val1+val2*cos(lonfi(ig)*1.57079633*2/pi-val3*pi/180.) 1704 ! If above line and ice: remove ice 1705 IF ( (latfi(ig)*180./pi.ge.val11) .and. 1706 & (latfi(ig)*180./pi.le.val1+val2) .and. 1707 & (qsurf(ig,igcm_n2).gt.0.) ) then 1708 ! Look for same longitude point where no ice 1709 val5=1. 1710 jref=ig 1711 ! 1 1712 ! ... iip1 ... x (jjp1-2) 32 x 23 1713 ! 1 1714 do while (val5.ge.1..and.jref.gt.iip1) 1715 ! northward point 1716 jref=jref-iip1+1 1717 ! ice at the point 1718 val5=qsurf(jref,igcm_n2) 1719 ! write(*,*) jref, 1720 ! & latfi(jref)*180./pi,lonfi(jref)*180/pi,val5 1721 enddo 1722 if (val5.ge.1) write(*,*) 'NO POINT FOUND WITH NO ICE' 1723 1724 ! Copy point in the selected area 1725 tsurf(ig)=tsurf(jref) 1726 qsurf(ig,igcm_n2)=qsurf(jref,igcm_n2) 1727 qsurf(ig,igcm_ch4_ice)=qsurf(jref,igcm_ch4_ice) 1728 qsurf(ig,igcm_co_ice)=qsurf(jref,igcm_co_ice) 1729 DO l=1,nsoilmx 1730 tsoil(ig,l)=tsoil(jref,l) 1731 ENDDO 1732 ENDIF 1733 ! If below line and no ice: add ice 1734 IF ( (latfi(ig)*180./pi.le.val11) .and. 1735 & (latfi(ig)*180./pi.le.val1+val2) .and. 1736 & (qsurf(ig,igcm_n2).eq.0.) ) then 1737 ! Look for same longitude point where ice 1738 val5=1. 1739 jref=ig 1740 do while (val5.le.1.and.jref.lt.ngridmx-iip1) 1741 ! southward point 1742 jref=jref+iip1-1 1743 ! ice at the point 1744 val5=qsurf(jref,igcm_n2) 1745 write(*,*) jref, 1746 & latfi(jref)*180./pi,lonfi(jref)*180/pi,val5 1747 enddo 1748 if (val5.le.1) write(*,*) 'NO POINT FOUND WITH ICE' 1749 1750 ! Copy point in the selected area 1751 tsurf(ig)=tsurf(jref) 1752 qsurf(ig,igcm_n2)=qsurf(jref,igcm_n2) 1753 qsurf(ig,igcm_ch4_ice)=qsurf(jref,igcm_ch4_ice) 1754 qsurf(ig,igcm_co_ice)=qsurf(jref,igcm_co_ice) 1755 DO l=1,nsoilmx 1756 tsoil(ig,l)=tsoil(jref,l) 1757 ENDDO 1758 ENDIF 1759 1760 ENDDO 1761 1762 c source_ch4 : adding source ch4 1763 c ---------------------------------------------------------- 1764 else if (modif(1:len_trim(modif)) .eq. 'source_ch4') then 1765 1766 write(*,*) 'Adding ch4 ice latitudinally ' 1767 write(*,*) 'Choice : lat1 and lat2 ?' 1768 433 read(*,*,iostat=ierr) val,val2 1769 if(ierr.ne.0) goto 433 1770 write(*,*) 'Choice : lon1 and lon2 ?' 1771 434 read(*,*,iostat=ierr) val3,val4 1772 if(ierr.ne.0) goto 434 1773 write(*,*) 'Choice amount (kg/m2)' 1774 435 read(*,*,iostat=ierr) val5 1775 if(ierr.ne.0) goto 435 1776 1777 DO ig=1,ngridmx 1778 IF ( (latfi(ig)*180./pi.ge.val) .and. 1779 & (latfi(ig)*180./pi.le.val2) .and. 1780 & (lonfi(ig)*180./pi.ge.val3) .and. 1781 & (lonfi(ig)*180./pi.lt.val4) ) then 1782 qsurf(ig,igcm_ch4_ice)=qsurf(ig,igcm_ch4_ice)+val5 1783 ENDIF 1784 ENDDO 1785 1786 c source_co : adding source co 1787 c ---------------------------------------------------------- 1788 else if (modif(1:len_trim(modif)) .eq. 'source_co') then 1789 1790 write(*,*) 'Adding co ice latitudinally ' 1791 write(*,*) 'Choice : lat1 and lat2 ?' 1792 436 read(*,*,iostat=ierr) val,val2 1793 if(ierr.ne.0) goto 436 1794 write(*,*) 'Choice : lon1 and lon2 ?' 1795 437 read(*,*,iostat=ierr) val3,val4 1796 if(ierr.ne.0) goto 437 1797 write(*,*) 'Choice amount (kg/m2)' 1798 438 read(*,*,iostat=ierr) val5 1799 if(ierr.ne.0) goto 438 1800 1801 DO ig=1,ngridmx 1802 IF ( (latfi(ig)*180./pi.ge.val) .and. 1803 & (latfi(ig)*180./pi.le.val2) .and. 1804 & (lonfi(ig)*180./pi.ge.val3) .and. 1805 & (lonfi(ig)*180./pi.lt.val4) ) then 1806 qsurf(ig,igcm_co_ice)=qsurf(ig,igcm_co_ice)+val5 1807 ENDIF 1808 ENDDO 1809 1810 ! therm_ini_s: (re)-set soil thermal inertia to reference surface values 1811 ! ---------------------------------------------------------------------- 1812 else if (modif(1:len_trim(modif)).eq.'therm_ini_s') then 1813 ! write(*,*)"surfithfi(1):",surfithfi(1) 1814 do isoil=1,nsoilmx 1815 inertiedat(1:ngridmx,isoil)=surfithfi(1:ngridmx) 1816 enddo 1817 write(*,*)'OK: Soil thermal inertia has been reset to referenc 1818 &e surface values' 1819 write(*,*)"inertiedat(1,1):",inertiedat(1,1) 1820 ithfi(:,:)=inertiedat(:,:) 1821 ! recast ithfi() onto ith() 1822 call gr_fi_dyn(nsoilmx,ngridmx,iip1,jjp1,ithfi,ith) 1823 1824 ! inert3d: set soil thermal inertia to specific uniform value 1825 ! ---------------------------------------------------------------------- 1826 else if (modif(1:len_trim(modif)).eq.'inert3d') then 1827 write(*,*) 'Actual value of surf Thermal inertia at ig=1: ', 1828 & inertiedat(1,1), ' SI' 1829 write(*,*) 'Value of Thermal inertia (SI) ? ' 1830 read(*,*) val 1831 do isoil=1,nsoilmx 1832 do ig=1,ngridmx 1833 inertiedat(ig,isoil)=val 1834 enddo 1835 enddo 1836 write(*,*)'OK: Soil TI set to ',val,' SI everywhere' 1837 ithfi(:,:)=inertiedat(:,:) 1838 ! recast ithfi() onto ith() 1839 call gr_fi_dyn(nsoilmx,ngridmx,iip1,jjp1,ithfi,ith) 1840 1841 ! subsoilice_n: Put deep ice layer in northern hemisphere soil 1842 ! ------------------------------------------------------------ 1843 else if (modif(1:len_trim(modif)).eq.'subsoilice_n') then 1844 1845 write(*,*)'From which latitude (in deg.), up to the north pole, 1846 &should we put subterranean ice?' 1847 ierr=1 1848 do while (ierr.ne.0) 1849 read(*,*,iostat=ierr) val 1850 if (ierr.eq.0) then ! got a value 1851 ! do a sanity check 1852 if((val.lt.0.).or.(val.gt.90)) then 1853 write(*,*)'Latitude should be between 0 and 90 deg. !!!' 849 1854 ierr=1 850 endif 851 if(val2.lt.layer(1)) then 852 if(val2.eq.0) then 853 write(*,*)'Thermal inertia set for all subsurface layers' 854 ierr=0 855 else 856 write(*,*)'Depth should be more than ',layer(1) 857 ierr=1 858 endif 859 endif 860 endif 861 enddo ! of do while 862 863 ! find the reference index iref the depth corresponds to 864 if(val2.eq.0) then 865 iref=1 866 write(*,*)'Level selected is first level: ',layer(iref),' m' 867 else 868 do isoil=1,nsoilmx-1 869 if ((val2.gt.layer(isoil)).and.(val2.lt.layer(isoil+1))) 870 & then 871 iref=isoil+1 872 write(*,*)'Level selected : ',layer(isoil+1),' m' 873 exit 874 endif 875 enddo 876 endif 877 878 DO ig=1,ngridmx 879 DO j=iref,nsoilmx 880 ithfi(ig,j)=ith_bb 881 ENDDO 882 ENDDO 883 884 CALL gr_fi_dyn(nsoilmx,ngridmx,iip1,jjp1,ithfi,ith) 885 886 c diurnal_TI : choice of thermal inertia values (global) 887 c ---------------------------------------------------------------- 888 else if (modif(1:len_trim(modif)) .eq. 'diurnal_TI') then 889 890 write(*,*) 'New value for diurnal thermal inertia ?' 891 106 read(*,*,iostat=ierr) ith_bb 892 if(ierr.ne.0) goto 106 893 write(*,*) 'Diurnal thermal inertia (new value):',ith_bb 894 895 write(*,*)'At which depth (in m.) does the ice layer ends?' 896 write(*,*)'(currently, the soil layer 1 and nsoil are:' 897 & ,layer(1),' - ',layer(nsoilmx),')' 1855 else ! find corresponding jref (nearest latitude) 1856 ! note: rlatu(:) contains decreasing values of latitude 1857 ! starting from PI/2 to -PI/2 1858 do j=1,jjp1 1859 if ((rlatu(j)*180./pi.ge.val).and. 1860 & (rlatu(j+1)*180./pi.le.val)) then 1861 ! find which grid point is nearest to val: 1862 if (abs(rlatu(j)*180./pi-val).le. 1863 & abs((rlatu(j+1)*180./pi-val))) then 1864 jref=j 1865 else 1866 jref=j+1 1867 endif 1868 1869 write(*,*)'Will use nearest grid latitude which is:', 1870 & rlatu(jref)*180./pi 1871 endif 1872 enddo ! of do j=1,jjp1 1873 endif ! of if((val.lt.0.).or.(val.gt.90)) 1874 endif !of if (ierr.eq.0) 1875 enddo ! of do while 1876 1877 ! Build layers() (as in soil_settings.F) 1878 val2=sqrt(mlayer(0)*mlayer(1)) 1879 val3=mlayer(1)/mlayer(0) 1880 do isoil=1,nsoilmx 1881 layer(isoil)=val2*(val3**(isoil-1)) 1882 enddo 1883 1884 write(*,*)'At which depth (in m.) does the ice layer begin?' 1885 write(*,*)'(currently, the deepest soil layer extends down to:' 1886 & ,layer(nsoilmx),')' 898 1887 ierr=1 899 1888 do while (ierr.ne.0) 900 1889 read(*,*,iostat=ierr) val2 901 write(*,*)'val2 in diurnal_TI:',val2,'ierr=',ierr1890 ! write(*,*)'val2:',val2,'ierr=',ierr 902 1891 if (ierr.eq.0) then ! got a value, but do a sanity check 903 1892 if(val2.gt.layer(nsoilmx)) then … … 911 1900 endif 912 1901 enddo ! of do while 913 914 ! find the reference index iref the depth corresponds to 1902 1903 ! find the reference index iref the depth corresponds to 1904 ! if (val2.lt.layer(1)) then 1905 ! iref=1 1906 ! else 915 1907 do isoil=1,nsoilmx-1 916 !write(*,*)'isoil, ',isoil,val2 917 !write(*,*)'lay(i),lay(i+1):',layer(isoil),layer(isoil+1),' m' 918 if ((val2.gt.layer(isoil)).and.(val2.lt.layer(isoil+1))) 1908 if((val2.gt.layer(isoil)).and.(val2.lt.layer(isoil+1))) 919 1909 & then 920 iref=isoil+1 921 write(*,*)'Level selected : ',layer(isoil+1),' m' 1910 iref=isoil 922 1911 exit 1912 endif 1913 enddo 1914 1915 ! thermal inertia of the ice: 1916 ierr=1 1917 do while (ierr.ne.0) 1918 write(*,*)'What is the value of subterranean ice thermal inert 1919 &ia? (e.g.: 2000)' 1920 read(*,*,iostat=ierr)iceith 1921 enddo ! of do while 1922 1923 ! recast ithfi() onto ith() 1924 call gr_fi_dyn(nsoilmx,ngridmx,iip1,jjp1,ithfi,ith) 1925 1926 do j=1,jref 1927 ! write(*,*)'j:',j,'rlatu(j)*180./pi:',rlatu(j)*180./pi 1928 do i=1,iip1 ! loop on longitudes 1929 ! Build "equivalent" thermal inertia for the mixed layer 1930 ith(i,j,iref+1)=sqrt((layer(iref+1)-layer(iref))/ 1931 & (((val2-layer(iref))/(ith(i,j,iref)**2))+ 1932 & ((layer(iref+1)-val2)/(iceith)**2))) 1933 ! Set thermal inertia of lower layers 1934 do isoil=iref+2,nsoilmx 1935 ith(i,j,isoil)=iceith ! ice 1936 enddo 1937 enddo ! of do i=1,iip1 1938 enddo ! of do j=1,jjp1 1939 1940 1941 CALL gr_dyn_fi(nsoilmx,iip1,jjp1,ngridmx,ith,ithfi) 1942 1943 ! subsoilice_s: Put deep ice layer in southern hemisphere soil 1944 ! ------------------------------------------------------------ 1945 else if (modif(1:len_trim(modif)).eq.'subsoilice_s') then 1946 1947 write(*,*)'From which latitude (in deg.), down to the south pol 1948 &e, should we put subterranean ice?' 1949 ierr=1 1950 do while (ierr.ne.0) 1951 read(*,*,iostat=ierr) val 1952 if (ierr.eq.0) then ! got a value 1953 ! do a sanity check 1954 if((val.gt.0.).or.(val.lt.-90)) then 1955 write(*,*)'Latitude should be between 0 and -90 deg. !!!' 1956 ierr=1 1957 else ! find corresponding jref (nearest latitude) 1958 ! note: rlatu(:) contains decreasing values of latitude 1959 ! starting from PI/2 to -PI/2 1960 do j=1,jjp1 1961 if ((rlatu(j)*180./pi.ge.val).and. 1962 & (rlatu(j+1)*180./pi.le.val)) then 1963 ! find which grid point is nearest to val: 1964 if (abs(rlatu(j)*180./pi-val).le. 1965 & abs((rlatu(j+1)*180./pi-val))) then 1966 jref=j 1967 else 1968 jref=j+1 1969 endif 1970 1971 write(*,*)'Will use nearest grid latitude which is:', 1972 & rlatu(jref)*180./pi 1973 endif 1974 enddo ! of do j=1,jjp1 1975 endif ! of if((val.lt.0.).or.(val.gt.90)) 1976 endif !of if (ierr.eq.0) 1977 enddo ! of do while 1978 1979 ! Build layers() (as in soil_settings.F) 1980 val2=sqrt(mlayer(0)*mlayer(1)) 1981 val3=mlayer(1)/mlayer(0) 1982 do isoil=1,nsoilmx 1983 layer(isoil)=val2*(val3**(isoil-1)) 1984 enddo 1985 1986 write(*,*)'At which depth (in m.) does the ice layer begin?' 1987 write(*,*)'(currently, the deepest soil layer extends down to:' 1988 & ,layer(nsoilmx),')' 1989 ierr=1 1990 do while (ierr.ne.0) 1991 read(*,*,iostat=ierr) val2 1992 ! write(*,*)'val2:',val2,'ierr=',ierr 1993 if (ierr.eq.0) then ! got a value, but do a sanity check 1994 if(val2.gt.layer(nsoilmx)) then 1995 write(*,*)'Depth should be less than ',layer(nsoilmx) 1996 ierr=1 923 1997 endif 1998 if(val2.lt.layer(1)) then 1999 write(*,*)'Depth should be more than ',layer(1) 2000 ierr=1 2001 endif 2002 endif 2003 enddo ! of do while 2004 2005 ! find the reference index iref the depth corresponds to 2006 do isoil=1,nsoilmx-1 2007 if((val2.gt.layer(isoil)).and.(val2.lt.layer(isoil+1))) 2008 & then 2009 iref=isoil 2010 exit 2011 endif 2012 enddo 2013 2014 ! thermal inertia of the ice: 2015 ierr=1 2016 do while (ierr.ne.0) 2017 write(*,*)'What is the value of subterranean ice thermal inert 2018 &ia? (e.g.: 2000)' 2019 read(*,*,iostat=ierr)iceith 2020 enddo ! of do while 2021 2022 ! recast ithfi() onto ith() 2023 call gr_fi_dyn(nsoilmx,ngridmx,iip1,jjp1,ithfi,ith) 2024 2025 do j=jref,jjp1 2026 ! write(*,*)'j:',j,'rlatu(j)*180./pi:',rlatu(j)*180./pi 2027 do i=1,iip1 ! loop on longitudes 2028 ! Build "equivalent" thermal inertia for the mixed layer 2029 ith(i,j,iref+1)=sqrt((layer(iref+1)-layer(iref))/ 2030 & (((val2-layer(iref))/(ith(i,j,iref)**2))+ 2031 & ((layer(iref+1)-val2)/(iceith)**2))) 2032 ! Set thermal inertia of lower layers 2033 do isoil=iref+2,nsoilmx 2034 ith(i,j,isoil)=iceith ! ice 2035 enddo 2036 enddo ! of do i=1,iip1 2037 enddo ! of do j=jref,jjp1 2038 2039 CALL gr_dyn_fi(nsoilmx,iip1,jjp1,ngridmx,ith,ithfi) 2040 2041 c ---------------------------------------------------------- 2042 c 'reservoir_SP' : increase or decrase ices reservoir in SP by factor 2043 c ---------------------------------------------------------- 2044 else if (modif(1:len_trim(modif)).eq.'reservoir_SP') then 2045 write(*,*) "Increase/Decrease N2 reservoir by factor :" 2046 read(*,*) val7 2047 write(*,*) "Increase/Decrease CH4 reservoir by factor :" 2048 read(*,*) val8 2049 write(*,*) "Increase/Decrease CO reservoir by factor :" 2050 read(*,*) val9 2051 2052 ! Definition SP: 2053 val3=-33. !lat1 2054 val4=50. !lat2 2055 val5=140. ! lon1 2056 val6=-155. ! lon2 2057 choice=-50. ! min phisfi 2058 write(*,*) 'def SP :' 2059 write(*,*) 'lat :',val3,val4 2060 write(*,*) 'lon :',val5,'180 / -180',val6 2061 write(*,*) 'choice phisfi min :',choice 2062 2063 DO ig=1,ngridmx 2064 2065 IF ( (latfi(ig)*180./pi.ge.val3) .and. 2066 & (latfi(ig)*180./pi.le.val4) .and. 2067 & ( ((lonfi(ig)*180./pi.ge.-180.) .and. 2068 & (lonfi(ig)*180./pi.le.val6)) .or. 2069 & ((lonfi(ig)*180./pi.ge.val5) .and. 2070 & (lonfi(ig)*180./pi.le.180.))) ) then 2071 2072 IF ((phisfi(ig).le.choice)) then !.and. 2073 c & (qsurf(ig,igcm_n2).ge.50)) then 2074 2075 qsurf(ig,igcm_n2)=qsurf(ig,igcm_n2)*val7 2076 qsurf(ig,igcm_ch4_ice)=qsurf(ig,igcm_ch4_ice)*val8 2077 qsurf(ig,igcm_co_ice)=qsurf(ig,igcm_co_ice)*val9 2078 ENDIF 2079 ENDIF 2080 ENDDO 2081 2082 c ---------------------------------------------------------- 2083 c 'reservoir' : increase or decrase ices reservoir by factor 2084 c ---------------------------------------------------------- 2085 else if (modif(1:len_trim(modif)).eq.'reservoir') then 2086 write(*,*) "Increase/Decrease N2 reservoir by factor :" 2087 read(*,*) val 2088 write(*,*) "Increase/Decrease CH4 reservoir by factor :" 2089 read(*,*) val2 2090 write(*,*) "Increase/Decrease CO reservoir by factor :" 2091 read(*,*) val3 2092 2093 DO ig=1,ngridmx 2094 qsurf(ig,igcm_n2)=qsurf(ig,igcm_n2)*val 2095 qsurf(ig,igcm_ch4_ice)=qsurf(ig,igcm_ch4_ice)*val2 2096 qsurf(ig,igcm_co_ice)=qsurf(ig,igcm_co_ice)*val3 2097 ENDDO 2098 2099 c -------------------------------------------------------- 2100 c 'plutomap' : initialize pluto ices distribution 2101 c -------------------------------------------------------- 2102 else if (modif(1:len_trim(modif)).eq.'plutomap') then 2103 2104 write(*,*) 'pluto_map.dat has to be in your simulation file !!' 2105 open(27,file='pluto_map.dat') 2106 N2_pluto_dat(:,:) =0. 2107 CH4_pluto_dat(:,:) =0. 2108 CO_pluto_dat(:,:) =0. 2109 2110 ! Dimension file pluto_map 2111 IF (jm_plu.ne.180) then 2112 write(*,*) 'Newstart : you should set jm_plu to 180' 2113 stop 2114 ENDIF 2115 ! Read values 2116 do j=1,jm_plu+1 2117 read(27,*) (map_pluto_dat(i,j),i=1,im_plu) 2118 do i=1,im_plu 2119 2120 !!!!!! BTD_v2 2121 if (map_pluto_dat(i,j).eq.3) then 2122 N2_pluto_dat(i,j)=100000. 2123 elseif (map_pluto_dat(i,j).eq.4) then 2124 CH4_pluto_dat(i,j)=100000. 2125 elseif (map_pluto_dat(i,j).eq.0) then 2126 alb_dat(i,j)=0.3 2127 elseif (map_pluto_dat(i,j).eq.6) then 2128 alb_dat(i,j)=0.6 2129 elseif (map_pluto_dat(i,j).eq.7) then 2130 alb_dat(i,j)=0.1 2131 endif 2132 2133 !!!!!! BTD 2134 !if (map_pluto_dat(i,j).eq.3) then 2135 ! CH4_pluto_dat(i,j)=100000. 2136 !endif 2137 end do 2138 end do 2139 close (27) 2140 ! Interpolate 2141 2142 !! latitude and longitude in REindexed pluto map : 2143 latv_plu(1)=90. *pi/180. 2144 do j=2,jm_plu 2145 latv_plu(j)= (90-(j-1 -0.5)*180./(jm_plu-1))*pi/180. 2146 end do 2147 latv_plu(jm_plu+1)= -90. *pi/180. 2148 2149 do i=1,im_plu+1 2150 lonu_plu(i)=(-180+ (i-1)*360./(im_plu) ) *pi/180. 924 2151 enddo 925 2152 2153 !Reindexation to shift the longitude grid like a LMD GCM grid... 2154 do j=1,jm_plu+1 2155 N2_pluto_rein(1,j)=(N2_pluto_dat(1,j)+ 2156 & N2_pluto_dat(im_plu,j))/2 2157 CH4_pluto_rein(1,j)=(CH4_pluto_dat(1,j)+ 2158 & CH4_pluto_dat(im_plu,j))/2 2159 CO_pluto_rein(1,j)=(CO_pluto_dat(1,j)+ 2160 & CO_pluto_dat(im_plu,j))/2 2161 alb_rein(1,j)=(alb_dat(1,j)+ 2162 & alb_dat(im_plu,j))/2 2163 do i=2,im_plu 2164 N2_pluto_rein(i,j)=(N2_pluto_dat(i-1,j)+ 2165 & N2_pluto_dat(i,j))/2 2166 CH4_pluto_rein(i,j)=(CH4_pluto_dat(i-1,j)+ 2167 & CH4_pluto_dat(i,j))/2 2168 CO_pluto_rein(i,j)=(CO_pluto_dat(i-1,j)+ 2169 & CO_pluto_dat(i,j))/2 2170 alb_rein(i,j)=(alb_dat(i-1,j)+ 2171 & alb_dat(i,j))/2 2172 end do 2173 N2_pluto_rein(im_plu+1,j) = N2_pluto_rein(1,j) 2174 CH4_pluto_rein(im_plu+1,j) = CH4_pluto_rein(1,j) 2175 CO_pluto_rein(im_plu+1,j) = CO_pluto_rein(1,j) 2176 alb_rein(im_plu+1,j) = alb_rein(1,j) 2177 end do 2178 2179 call interp_horiz(N2_pluto_rein,N2_pluto_gcm, 2180 & im_plu,jm_plu,iim,jjm,1,lonu_plu,latv_plu,rlonu,rlatv) 2181 call interp_horiz(CH4_pluto_rein,CH4_pluto_gcm, 2182 & im_plu,jm_plu,iim,jjm,1,lonu_plu,latv_plu,rlonu,rlatv) 2183 call interp_horiz(CO_pluto_rein,CO_pluto_gcm, 2184 & im_plu,jm_plu,iim,jjm,1,lonu_plu,latv_plu,rlonu,rlatv) 2185 call interp_horiz(alb_rein,alb_gcm, 2186 & im_plu,jm_plu,iim,jjm,1,lonu_plu,latv_plu,rlonu,rlatv) 2187 2188 ! Switch dyn grid to fi grid 2189 qsurf_tmp(:,:) =0. 2190 CALL gr_dyn_fi(1,iip1,jjp1,ngridmx, 2191 & N2_pluto_gcm,qsurf_tmp(1,igcm_n2)) 2192 if (igcm_ch4_ice.ne.0) then 2193 CALL gr_dyn_fi(1,iip1,jjp1,ngridmx, 2194 & CH4_pluto_gcm,qsurf_tmp(1,igcm_ch4_ice)) 2195 endif 2196 if (igcm_co_ice.ne.0) then 2197 CALL gr_dyn_fi(1,iip1,jjp1,ngridmx, 2198 & CO_pluto_gcm,qsurf_tmp(1,igcm_co_ice)) 2199 endif 2200 alb=alb_gcm 2201 CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,alb_gcm,albfi) 2202 !print*, 'N2dat=',N2_pluto_gcm 2203 !print*, 'COdat=',CO_pluto_gcm 2204 !print*, 'CH4dat=',CH4_pluto_gcm 2205 print*, 'N2dat=',qsurf_tmp(:,igcm_n2) 2206 print*, 'COdat=',qsurf_tmp(:,igcm_co_ice) 2207 print*, 'CH4dat=',qsurf_tmp(:,igcm_ch4_ice) 2208 2209 ! Fill 926 2210 DO ig=1,ngridmx 927 DO j=1,iref 928 ithfi(ig,j)=ith_bb 929 ENDDO 930 ENDDO 931 932 CALL gr_fi_dyn(nsoilmx,ngridmx,iip1,jjp1,ithfi,ith) 933 934 935 936 937 938 939 940 941 942 943 944 945 946 947 948 949 950 951 952 953 954 955 2211 qsurf(ig,igcm_n2)=qsurf(ig,igcm_n2)+qsurf_tmp(ig,igcm_n2) 2212 if (igcm_ch4_ice.ne.0) then 2213 qsurf(ig,igcm_ch4_ice)=qsurf(ig,igcm_ch4_ice)+ 2214 & qsurf_tmp(ig,igcm_ch4_ice) 2215 endif 2216 if (igcm_co_ice.ne.0) then 2217 qsurf(ig,igcm_co_ice)=qsurf(ig,igcm_co_ice)+ 2218 & qsurf_tmp(ig,igcm_co_ice) 2219 endif 2220 ENDDO 956 2221 957 2222 endif … … 961 2226 999 continue 962 2227 963 964 c=======================================================================965 c Initialisation for cloud fraction and oceanic ice (added by BC 2010)966 c=======================================================================967 ! DO ig=1, ngridmx968 ! totalfrac(ig)=0.5969 ! DO l=1,llm970 ! cloudfrac(ig,l)=0.5971 ! ENDDO972 ! Ehouarn, march 2012: also add some initialisation for hice973 ! hice(ig)=0.0974 ! ENDDO975 976 c========================================================977 978 ! DO ig=1,ngridmx979 ! IF(tsurf(ig) .LT. 273.16-1.8) THEN980 ! hice(ig)=(273.16-1.8-tsurf(ig))/(273.16-1.8-240)*1.981 ! hice(ig)=min(hice(ig),1.0)982 ! ENDIF983 ! ENDDO984 985 986 987 988 2228 c======================================================================= 989 2229 c Correct pressure on the new grid (menu 0) 990 2230 c======================================================================= 991 992 2231 993 2232 if ((choix_1.eq.0).and.(.not.flagps0)) then … … 1008 2247 end if 1009 2248 1010 1011 c=======================================================================1012 c=======================================================================1013 1014 2249 c======================================================================= 1015 2250 c Initialisation de la physique / ecriture de newstartfi : 1016 2251 c======================================================================= 1017 1018 2252 1019 2253 CALL inifilr … … 1111 2345 end 1112 2346 2347 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 2348 function dist_pluto(lat1,lon1,lat2,lon2) 2349 implicit none 2350 real dist_pluto 2351 real lat1,lon1,lat2,lon2 2352 real dlon, dlat 2353 real a,c 2354 real rad 2355 rad=1190. ! km 2356 2357 dlon = lon2 - lon1 2358 dlat = lat2 - lat1 2359 a = (sin(dlat/2))**2 + cos(lat1) * cos(lat2) * (sin(dlon/2))**2 2360 c = 2 * atan2( sqrt(a), sqrt(1-a) ) 2361 dist_pluto = rad * c 2362 return 2363 end 2364
Note: See TracChangeset
for help on using the changeset viewer.
