| 1 | PROGRAM SUPERFIT |
|---|
| 2 | |
|---|
| 3 | c ------------------------------------------------------------------ |
|---|
| 4 | c #### Special version with ice depth and ice inertia ##### |
|---|
| 5 | c Program used to compute best-fit cap albedo, ICE DEPTH and total CO2 |
|---|
| 6 | c inventory |
|---|
| 7 | c Make use of 2 annual GCM run performed with 2 different cap albedos. |
|---|
| 8 | c Based on Hourdin et al. (JGR, 1995) |
|---|
| 9 | c Modified to account for non-linear sensitivity of polar cap mass |
|---|
| 10 | c To cap albedo (FF,2000) |
|---|
| 11 | c Modified/cleaned-up to use 4 input files and minimise cap albedos, |
|---|
| 12 | c MONS-derived ice depth and total pressure (EM,2008) |
|---|
| 13 | c Modified to minimise wrt ice thermal inertia coefficients (EM, 2009) |
|---|
| 14 | c ------------------------------------------------------------------ |
|---|
| 15 | |
|---|
| 16 | implicit none |
|---|
| 17 | |
|---|
| 18 | integer ngcm,nvl1,nsol,time_unit |
|---|
| 19 | |
|---|
| 20 | c ***** THE FOLLOWING LINES MUST BE ADAPTED FOR EACH CASES : ************** |
|---|
| 21 | ! parameter(ngcm=10704) ! size of the 1 year raw GCM data |
|---|
| 22 | ! parameter(ngcm=2672) ! "*4 & *6" files |
|---|
| 23 | parameter(ngcm=8028) ! "*3 & *5" files |
|---|
| 24 | parameter(nvl1=669) ! size of the Vl1 Ps Observation |
|---|
| 25 | c ************************************************************* |
|---|
| 26 | |
|---|
| 27 | |
|---|
| 28 | c size the smoothed data (1/sol) which are used for the simulation |
|---|
| 29 | parameter(nsol=669) ! size the smoothed data (1/sol) |
|---|
| 30 | c parameter(nsol=445) ! size the smoothed data (1/sol) |
|---|
| 31 | |
|---|
| 32 | c OBSERVED VL1 pressure |
|---|
| 33 | c ~~~~~~~~~~~~~~~~~~~~~ |
|---|
| 34 | real pvl_obs(nvl1), solvl(nvl1), lsvl(nvl1) |
|---|
| 35 | |
|---|
| 36 | c Raw Data from 4 simulations with 2 different albedo/iced : |
|---|
| 37 | c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
|---|
| 38 | c runs are numbered 1-4 such than runs 1&2 are at same ice depth |
|---|
| 39 | c and runs 3&4 are at same ice depth |
|---|
| 40 | c and runs 1&3 are at same ice inertia |
|---|
| 41 | c and runs 2&4 are at same ice inertia |
|---|
| 42 | real iceigcmN(4) , iceigcmS(4) |
|---|
| 43 | real icedgcm(4) |
|---|
| 44 | integer refrun ! run number of the 'reference' run |
|---|
| 45 | ! reference values for functions (e.g. values of run #1) |
|---|
| 46 | real iceigcm_refN, icedgcm_ref |
|---|
| 47 | real iceigcm_refS |
|---|
| 48 | |
|---|
| 49 | real ptot(4) |
|---|
| 50 | real patm(ngcm,4), pcapn(ngcm,4),pcaps(ngcm,4) |
|---|
| 51 | real pvl_gcm(ngcm,4) ! Simulated VL1 pressure |
|---|
| 52 | real solgcm(ngcm) ! time in runs in sol |
|---|
| 53 | real lsgcm(ngcm) ! time in runs in ls |
|---|
| 54 | integer year_xvik ! year of xvik files to use for fit |
|---|
| 55 | integer plot_test ! =1 if output file is to plot COST(DN,DS), =2 for COST(IN,IS)s |
|---|
| 56 | |
|---|
| 57 | c Smoothed Data from the 2 simulations with 2 different albedo : |
|---|
| 58 | c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
|---|
| 59 | real box |
|---|
| 60 | integer n ,i |
|---|
| 61 | real sol(nsol) ! time in runs |
|---|
| 62 | real ls(nsol) |
|---|
| 63 | real solconv |
|---|
| 64 | real pvl_sm(nsol,4) ! Simulated smooth VL1 pressure |
|---|
| 65 | real patm_sm(nsol,4), pcapn_sm(nsol,4),pcaps_sm(nsol,4) |
|---|
| 66 | real alphavl1(nsol) ! pvl/patm ratio |
|---|
| 67 | c Mass cap sensitivity to thermal inertia (derivative) |
|---|
| 68 | real dpdicein(nsol),dpdiceis(nsol) |
|---|
| 69 | real dpdicein_fltr(nsol),dpdiceis_fltr(nsol) |
|---|
| 70 | c Mass cap sensitivity to ice inertia (derivative) |
|---|
| 71 | real dpden(nsol),dpdes(nsol) |
|---|
| 72 | real dpden_fltr(nsol),dpdes_fltr(nsol) |
|---|
| 73 | |
|---|
| 74 | real iceis, icein, ptry ,pvl1 |
|---|
| 75 | real iceds, icedn |
|---|
| 76 | ! ice inertia range to explore, and maximum allowed difference between N&S ice inertias: |
|---|
| 77 | real iceimin,iceimax,maxiceidiff |
|---|
| 78 | real deltaicei ! step in ice thermal inertia |
|---|
| 79 | ! ice depth range to explore, and maximum allowed difference between N&S ice depths: |
|---|
| 80 | real icedmin,icedmax,maxiceddiff |
|---|
| 81 | real deltaiced ! step in ice depth coefficient |
|---|
| 82 | ! total pressure range to explore, and pressure step: |
|---|
| 83 | real pmin,pmax,deltap |
|---|
| 84 | real cost,cost4plot |
|---|
| 85 | real iceis4plot,ps4plot,icein4plot !,iceds4plot |
|---|
| 86 | real iceds4plot,icedn4plot |
|---|
| 87 | real costmin, pfit, iceinfit, iceisfit |
|---|
| 88 | real pcapn_new, pcaps_new |
|---|
| 89 | real icednfit, icedsfit |
|---|
| 90 | logical fit |
|---|
| 91 | |
|---|
| 92 | real fonc, fonc2n, fonc2s |
|---|
| 93 | |
|---|
| 94 | c Inputs files (from xvik) variables : |
|---|
| 95 | c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
|---|
| 96 | character(len=33) filename1, filename2 |
|---|
| 97 | character(len=99) dset_DminImin,dset_DmaxImin |
|---|
| 98 | character(len=99) dset_DminImax,dset_DmaxImax |
|---|
| 99 | integer ie |
|---|
| 100 | |
|---|
| 101 | |
|---|
| 102 | c---------------------------------------------------------------------- |
|---|
| 103 | c Initialisation : |
|---|
| 104 | c ~~~~~~~~~~~~~~ |
|---|
| 105 | c Implicit function : behavior of Pcap = fonc(alb) |
|---|
| 106 | c (e.g. fonc(alb) = alb if linear) |
|---|
| 107 | !! fonc(albn) = exp(3.46*albn) ! for albedo, both north & south |
|---|
| 108 | fonc(icein) = icein |
|---|
| 109 | c fonc2n(icedn) = log(icedn +10) |
|---|
| 110 | fonc2n(icedn) = icedn ! for northern ice depth coefficient |
|---|
| 111 | c fonc2s(iceds) = log(iceds +1) |
|---|
| 112 | fonc2s(iceds) = iceds ! for southern ice depth coefficient |
|---|
| 113 | |
|---|
| 114 | c ***** THE FOLLOWING LINES MUST BE ADAPTED FOR EACH CASES : ************** |
|---|
| 115 | |
|---|
| 116 | write(*,*) 'Program written for xvik outputs files ''xpsol'' ', |
|---|
| 117 | &'and ''xprestot'' from xvik.F' |
|---|
| 118 | write(*,*) 'Enter extreme parameters values used in xvik.F' |
|---|
| 119 | write(*,*) 'Minimum ice depth coefficient' |
|---|
| 120 | read(*,*) icedmin |
|---|
| 121 | write(*,*) icedmin |
|---|
| 122 | write(*,*) 'Maximum ice depth coefficient' |
|---|
| 123 | read(*,*) icedmax |
|---|
| 124 | write(*,*) icedmax |
|---|
| 125 | write(*,*) 'Minimum ice thermal inertia ' |
|---|
| 126 | read(*,*) iceimin |
|---|
| 127 | write(*,*) iceimin |
|---|
| 128 | write(*,*) 'Maximum ice thermal inertia ' |
|---|
| 129 | read(*,*) iceimax |
|---|
| 130 | write(*,*) iceimax |
|---|
| 131 | |
|---|
| 132 | |
|---|
| 133 | maxiceidiff = 3000 ! maximum allowed difference between N and S best fit ice thermal inertia |
|---|
| 134 | deltaicei= 20 ! step in ice thermal inertia |
|---|
| 135 | maxiceddiff = 30.e-4 ! maximum allowed difference between N and S best fit ice depth |
|---|
| 136 | deltaiced= 0.5e-4 ! step in ice depth coefficient |
|---|
| 137 | pmin = 690 ! minimum allowed best fit CO2 total pressure (Pa) |
|---|
| 138 | pmax = 720 ! maximum allowed best fit CO2 total pressure (Pa) |
|---|
| 139 | deltap=1.0 ! step in pressure |
|---|
| 140 | |
|---|
| 141 | c File to be read : |
|---|
| 142 | c ~~~~~~~~~~~~~~~ |
|---|
| 143 | c runs are numbered 1-4 such than runs 1&2 are at same ice depth |
|---|
| 144 | c and runs 3&4 are at same ice depth |
|---|
| 145 | c and runs 1&3 are at same ice inertia |
|---|
| 146 | c and runs 2&4 are at same ice inertia |
|---|
| 147 | ! First run: |
|---|
| 148 | iceigcmN(1)= iceimin ! Northern Cap ice thermal inertia |
|---|
| 149 | iceigcmS(1)= iceimin ! Southern Cap ice thermal inertia |
|---|
| 150 | icedgcm(1) = icedmin ! Ice depth coefficient |
|---|
| 151 | ! Second run: |
|---|
| 152 | iceigcmN(2) = iceimax ! Northern Cap ice thermal inertia |
|---|
| 153 | iceigcmS(2) = iceimax ! Southern Cap ice thermal inertia |
|---|
| 154 | icedgcm(2) = icedmin ! Ice depth coefficient |
|---|
| 155 | ! Third run: |
|---|
| 156 | iceigcmN(3) = iceimin ! Northern Cap ice thermal inertia |
|---|
| 157 | iceigcmS(3) = iceimin ! Southern Cap ice thermal inertia |
|---|
| 158 | icedgcm(3) = icedmax ! Ice depth coefficient |
|---|
| 159 | ! Fourth run: |
|---|
| 160 | iceigcmN(4) = iceimax ! Northern Cap ice thermal inertia |
|---|
| 161 | iceigcmS(4) = iceimax ! Southern Cap ice thermal inertia |
|---|
| 162 | icedgcm(4) = icedmax ! Ice depth coefficient |
|---|
| 163 | |
|---|
| 164 | ! set reference values: those of one of the files (here file 1) |
|---|
| 165 | refrun=1 |
|---|
| 166 | ! ice inertia of the reference run used to simulate function |
|---|
| 167 | iceigcm_refN=iceigcmN(refrun) |
|---|
| 168 | ! ice inertia of the reference run used to simulate function |
|---|
| 169 | iceigcm_refS=iceigcmS(refrun) |
|---|
| 170 | ! icedepth coefficient of the reference run used to simulate function |
|---|
| 171 | icedgcm_ref=icedgcm(refrun) |
|---|
| 172 | |
|---|
| 173 | write(*,*) 'Path to xvik outputs with minimum ice depth ', |
|---|
| 174 | &'coefficient and minimum ice thermal inertia ?' |
|---|
| 175 | read (*,'(a)') dset_DminImin |
|---|
| 176 | write(*,*) dset_DminImin |
|---|
| 177 | write(*,*) 'Path to xvik outputs with minimum ice depth ', |
|---|
| 178 | &'coefficient and maximum ice thermal inertia ?' |
|---|
| 179 | read (*,'(a)') dset_DminImax |
|---|
| 180 | write(*,*) dset_DminImax |
|---|
| 181 | write(*,*) 'Path to xvik outputs with maximum ice depth ', |
|---|
| 182 | &'coefficient and minimum ice thermal inertia ?' |
|---|
| 183 | read (*,'(a)') dset_DmaxImin |
|---|
| 184 | write(*,*) dset_DmaxImin |
|---|
| 185 | write(*,*) 'Path to xvik outputs with maximum ice depth ', |
|---|
| 186 | &'coefficient and maximum ice thermal inertia ?' |
|---|
| 187 | read (*,'(a)') dset_DmaxImax |
|---|
| 188 | write(*,*) dset_DmaxImax |
|---|
| 189 | |
|---|
| 190 | write(*,*) 'Which year of xvik outputs do you want to use ?' |
|---|
| 191 | read (*,*) year_xvik |
|---|
| 192 | write(*,*) year_xvik |
|---|
| 193 | |
|---|
| 194 | write(*,*) 'Xvik files in sol (1), ls (2) or both (3)' |
|---|
| 195 | read(*,*) time_unit |
|---|
| 196 | write(*,*) time_unit |
|---|
| 197 | |
|---|
| 198 | |
|---|
| 199 | write(*,*) 'COST being the model/obs difference' |
|---|
| 200 | write(*,*) 'IN ans IS being the ice thermal inertia of ', |
|---|
| 201 | &'Northern and Southern Cap' |
|---|
| 202 | write(*,*) 'DN ans DS being the ice depth coefficient of ', |
|---|
| 203 | &'of Northern and Southern Cap' |
|---|
| 204 | write(*,*) 'Do you want to use output file ', |
|---|
| 205 | &'''minimization.txt'' to plot COST(DN,DS) (1) or COST(IN,IS) (2)' |
|---|
| 206 | read(*,*) plot_test |
|---|
| 207 | write(*,*) plot_test |
|---|
| 208 | |
|---|
| 209 | c For four simulation, files : "time (VL1 sol) , Ps Vl1 (Pa)" |
|---|
| 210 | c "time (VL1 sol), Patm,PcapN,PcapS (Pa)" |
|---|
| 211 | |
|---|
| 212 | |
|---|
| 213 | |
|---|
| 214 | write(filename1,fmt='(a6,i1)') 'xpsol1',year_xvik |
|---|
| 215 | write(filename2,fmt='(a8,i1)') 'xprestot',year_xvik |
|---|
| 216 | |
|---|
| 217 | |
|---|
| 218 | ! Dmin Imin files |
|---|
| 219 | |
|---|
| 220 | write(*,*) 'Opening ', trim(dset_DminImin)//'/'//trim(filename1) |
|---|
| 221 | open(21,file=trim(dset_DminImin)//'/'//trim(filename1),iostat=ie) |
|---|
| 222 | |
|---|
| 223 | if (ie.ne.0) then |
|---|
| 224 | write(*,*) 'Error opening file ',trim(dset_DminImin)//'/' |
|---|
| 225 | &//trim(filename1) |
|---|
| 226 | stop |
|---|
| 227 | endif |
|---|
| 228 | |
|---|
| 229 | write(*,*) 'Opening ', trim(dset_DminImin)//'/'//trim(filename2) |
|---|
| 230 | open(11,file=trim(dset_DminImin)//'/'//trim(filename2),iostat=ie) |
|---|
| 231 | |
|---|
| 232 | if (ie.ne.0) then |
|---|
| 233 | write(*,*) 'Error opening file ',trim(dset_DminImin)//'/' |
|---|
| 234 | &//trim(filename2) |
|---|
| 235 | stop |
|---|
| 236 | endif |
|---|
| 237 | |
|---|
| 238 | ! Dmin Imax files |
|---|
| 239 | |
|---|
| 240 | write(*,*) 'Opening ', trim(dset_DminImax)//'/'//trim(filename1) |
|---|
| 241 | open(22,file=trim(dset_DminImax)//'/'//trim(filename1),iostat=ie) |
|---|
| 242 | |
|---|
| 243 | if (ie.ne.0) then |
|---|
| 244 | write(*,*) 'Error opening file ',trim(dset_DminImax)//'/' |
|---|
| 245 | &//trim(filename1) |
|---|
| 246 | stop |
|---|
| 247 | endif |
|---|
| 248 | |
|---|
| 249 | write(*,*) 'Opening ', trim(dset_DminImax)//'/'//trim(filename2) |
|---|
| 250 | open(12,file=trim(dset_DminImax)//'/'//trim(filename2),iostat=ie) |
|---|
| 251 | |
|---|
| 252 | if (ie.ne.0) then |
|---|
| 253 | write(*,*) 'Error opening file ',trim(dset_DminImax)//'/' |
|---|
| 254 | &//trim(filename2) |
|---|
| 255 | stop |
|---|
| 256 | endif |
|---|
| 257 | |
|---|
| 258 | ! Dmax Imin files |
|---|
| 259 | |
|---|
| 260 | write(*,*) 'Opening ', trim(dset_DmaxImin)//'/'//trim(filename1) |
|---|
| 261 | open(23,file=trim(dset_DmaxImin)//'/'//trim(filename1),iostat=ie) |
|---|
| 262 | |
|---|
| 263 | if (ie.ne.0) then |
|---|
| 264 | write(*,*) 'Error opening file ',trim(dset_DmaxImin)//'/' |
|---|
| 265 | &//trim(filename1) |
|---|
| 266 | stop |
|---|
| 267 | endif |
|---|
| 268 | |
|---|
| 269 | write(*,*) 'Opening ', trim(dset_DmaxImin)//'/'//trim(filename2) |
|---|
| 270 | open(13,file=trim(dset_DmaxImin)//'/'//trim(filename2),iostat=ie) |
|---|
| 271 | |
|---|
| 272 | if (ie.ne.0) then |
|---|
| 273 | write(*,*) 'Error opening file ',trim(dset_DmaxImin)//'/' |
|---|
| 274 | &//trim(filename2) |
|---|
| 275 | stop |
|---|
| 276 | endif |
|---|
| 277 | |
|---|
| 278 | |
|---|
| 279 | ! Dmax Imax files |
|---|
| 280 | |
|---|
| 281 | write(*,*) 'Opening ', trim(dset_DmaxImax)//'/'//trim(filename1) |
|---|
| 282 | open(24,file=trim(dset_DmaxImax)//'/'//trim(filename1),iostat=ie) |
|---|
| 283 | |
|---|
| 284 | if (ie.ne.0) then |
|---|
| 285 | write(*,*) 'Error opening file ',trim(dset_DmaxImax)//'/' |
|---|
| 286 | &//trim(filename1) |
|---|
| 287 | stop |
|---|
| 288 | endif |
|---|
| 289 | |
|---|
| 290 | write(*,*) 'Opening ', trim(dset_DmaxImax)//'/'//trim(filename2) |
|---|
| 291 | open(14,file=trim(dset_DmaxImax)//'/'//trim(filename2),iostat=ie) |
|---|
| 292 | |
|---|
| 293 | if (ie.ne.0) then |
|---|
| 294 | write(*,*) 'Error opening file ',trim(dset_DmaxImax)//'/' |
|---|
| 295 | &//trim(filename2) |
|---|
| 296 | stop |
|---|
| 297 | endif |
|---|
| 298 | |
|---|
| 299 | |
|---|
| 300 | |
|---|
| 301 | c ************************************************************* |
|---|
| 302 | c Observed Smooth Viking Curves |
|---|
| 303 | c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
|---|
| 304 | open(30,file='VL1') |
|---|
| 305 | do n=1,nsol |
|---|
| 306 | read(30,*) solvl(n), lsvl(n), pvl_obs(n) |
|---|
| 307 | end do |
|---|
| 308 | close(30) |
|---|
| 309 | |
|---|
| 310 | c Opening output file |
|---|
| 311 | |
|---|
| 312 | !open(33, file = 'minimization.txt') |
|---|
| 313 | |
|---|
| 314 | |
|---|
| 315 | c reading Simulation results |
|---|
| 316 | c ~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
|---|
| 317 | |
|---|
| 318 | |
|---|
| 319 | do n=1,ngcm |
|---|
| 320 | |
|---|
| 321 | c Reading Viking Lander 1 simulated pressure for 4 runs |
|---|
| 322 | |
|---|
| 323 | if (time_unit == 1) then |
|---|
| 324 | read(21,*) solgcm(n), pvl_gcm(n,1) |
|---|
| 325 | read(22,*) solgcm(n), pvl_gcm(n,2) |
|---|
| 326 | read(23,*) solgcm(n), pvl_gcm(n,3) |
|---|
| 327 | read(24,*) solgcm(n), pvl_gcm(n,4) |
|---|
| 328 | |
|---|
| 329 | c Reading atmospheric and "cap" total pressure for all runs |
|---|
| 330 | read(11,*) solgcm(n), patm(n,1), pcapn(n,1),pcaps(n,1) |
|---|
| 331 | read(12,*) solgcm(n), patm(n,2), pcapn(n,2),pcaps(n,2) |
|---|
| 332 | read(13,*) solgcm(n), patm(n,3), pcapn(n,3),pcaps(n,3) |
|---|
| 333 | read(14,*) solgcm(n), patm(n,4), pcapn(n,4),pcaps(n,4) |
|---|
| 334 | |
|---|
| 335 | elseif (time_unit == 2) then |
|---|
| 336 | read(21,*) lsgcm(n), pvl_gcm(n,1) |
|---|
| 337 | read(22,*) lsgcm(n), pvl_gcm(n,2) |
|---|
| 338 | read(23,*) lsgcm(n), pvl_gcm(n,3) |
|---|
| 339 | read(24,*) lsgcm(n), pvl_gcm(n,4) |
|---|
| 340 | |
|---|
| 341 | c Reading atmospheric and "cap" total pressure for all runs |
|---|
| 342 | read(11,*) lsgcm(n), patm(n,1), pcapn(n,1),pcaps(n,1) |
|---|
| 343 | read(12,*) lsgcm(n), patm(n,2), pcapn(n,2),pcaps(n,2) |
|---|
| 344 | read(13,*) lsgcm(n), patm(n,3), pcapn(n,3),pcaps(n,3) |
|---|
| 345 | read(14,*) lsgcm(n), patm(n,4), pcapn(n,4),pcaps(n,4) |
|---|
| 346 | call ls2sol(lsgcm(n),solconv) |
|---|
| 347 | solgcm(n) = solconv |
|---|
| 348 | |
|---|
| 349 | elseif (time_unit == 3) then |
|---|
| 350 | read(21,*) solgcm(n), lsgcm(n), pvl_gcm(n,1) |
|---|
| 351 | read(22,*) solgcm(n), lsgcm(n), pvl_gcm(n,2) |
|---|
| 352 | read(23,*) solgcm(n), lsgcm(n), pvl_gcm(n,3) |
|---|
| 353 | read(24,*) solgcm(n), lsgcm(n), pvl_gcm(n,4) |
|---|
| 354 | |
|---|
| 355 | c Reading atmospheric and "cap" total pressure for all runs |
|---|
| 356 | read(11,*) solgcm(n), lsgcm(n), patm(n,1), pcapn(n,1),pcaps(n,1) |
|---|
| 357 | read(12,*) solgcm(n), lsgcm(n), patm(n,2), pcapn(n,2),pcaps(n,2) |
|---|
| 358 | read(13,*) solgcm(n), lsgcm(n), patm(n,3), pcapn(n,3),pcaps(n,3) |
|---|
| 359 | read(14,*) solgcm(n), lsgcm(n), patm(n,4), pcapn(n,4),pcaps(n,4) |
|---|
| 360 | |
|---|
| 361 | else |
|---|
| 362 | write(*,*) 'Wrong integer for xvik files format :', |
|---|
| 363 | &' must be 1, 2 or 3' |
|---|
| 364 | stop |
|---|
| 365 | |
|---|
| 366 | endif |
|---|
| 367 | |
|---|
| 368 | c Checking total CO2 inventory for all runs : |
|---|
| 369 | do i=1,4 |
|---|
| 370 | if(n.eq.1) then |
|---|
| 371 | ptot(i) = patm(n,i)+pcapn(n,i)+pcaps(n,i) |
|---|
| 372 | write(*,*) 'For run = ',i,' Ptot= ',ptot(i) |
|---|
| 373 | else |
|---|
| 374 | if(abs(patm(n,i)+pcapn(n,i)+pcaps(n,i)-ptot(i)).gt.3)then |
|---|
| 375 | write(*,*)'total pressure not constant for run i= ',i |
|---|
| 376 | write(*,*) 'n=',1,' ptot=',ptot(i) |
|---|
| 377 | write(*,*)'n=',n,' ptot=',patm(n,i)+pcapn(n,i)+pcaps(n,i) |
|---|
| 378 | end if |
|---|
| 379 | end if |
|---|
| 380 | end do |
|---|
| 381 | end do |
|---|
| 382 | |
|---|
| 383 | close(11) |
|---|
| 384 | close(12) |
|---|
| 385 | close(13) |
|---|
| 386 | close(14) |
|---|
| 387 | close(21) |
|---|
| 388 | close(22) |
|---|
| 389 | close(23) |
|---|
| 390 | close(24) |
|---|
| 391 | |
|---|
| 392 | |
|---|
| 393 | c Smoothing simulated GCM Viking 1 pressure curves |
|---|
| 394 | c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
|---|
| 395 | do n=1,nsol |
|---|
| 396 | sol(n)=float(n) |
|---|
| 397 | end do |
|---|
| 398 | |
|---|
| 399 | |
|---|
| 400 | c Running average for both runs (with "boxsize" box in days) |
|---|
| 401 | box=20. |
|---|
| 402 | do i=1,4 |
|---|
| 403 | call runave(solgcm,pvl_gcm(1,i),ngcm,669.,box, |
|---|
| 404 | & sol,pvl_sm(1,i),nsol) |
|---|
| 405 | call runave(solgcm,patm(1,i),ngcm,669.,box, |
|---|
| 406 | & sol,patm_sm(1,i),nsol) |
|---|
| 407 | call runave(solgcm,pcapn(1,i),ngcm,669.,box, |
|---|
| 408 | & sol,pcapn_sm(1,i),nsol) |
|---|
| 409 | call runave(solgcm,pcaps(1,i),ngcm,669.,box, |
|---|
| 410 | & sol,pcaps_sm(1,i),nsol) |
|---|
| 411 | end do |
|---|
| 412 | |
|---|
| 413 | |
|---|
| 414 | c Computing mass cap sensitivity to albedo (derivative) and alphaVl1 |
|---|
| 415 | c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
|---|
| 416 | c alphaVL1 = Pvl1/patm at a given time (see Hourdin et al., JGR, 1995) |
|---|
| 417 | |
|---|
| 418 | do n=1,nsol |
|---|
| 419 | ! Evaluate derivative of pressure wrt Northern inertia, as the |
|---|
| 420 | ! mean of the 2 values obtained using runs 1&2 and 3&4 |
|---|
| 421 | dpdicein(n)=(((pcapn_sm(n,1)-pcapn_sm(n,2))/ |
|---|
| 422 | & (fonc(iceigcmN(1))-fonc(iceigcmN(2))))+ |
|---|
| 423 | & ((pcapn_sm(n,3)-pcapn_sm(n,4))/ |
|---|
| 424 | & (fonc(iceigcmN(3))-fonc(iceigcmN(4)))))*(1./2.) |
|---|
| 425 | ! Evaluate derivative of pressure wrt Southern inertia, as the |
|---|
| 426 | ! mean of the 2 values obtained using runs 1&2 and 3&4 |
|---|
| 427 | dpdiceis(n)=(((pcaps_sm(n,1)-pcaps_sm(n,2))/ |
|---|
| 428 | & (fonc(iceigcmS(1))-fonc(iceigcmS(2))))+ |
|---|
| 429 | & ((pcaps_sm(n,3)-pcaps_sm(n,4))/ |
|---|
| 430 | & (fonc(iceigcmS(3))-fonc(iceigcmS(4)))))*(1./2.) |
|---|
| 431 | ! Evaluate derivative of pressure wrt Northern ice depth coefficient, |
|---|
| 432 | ! as the mean of the 2 values obtained using runs 1&3 and 2&4 |
|---|
| 433 | dpden(n)=(((pcapn_sm(n,1)-pcapn_sm(n,3))/ |
|---|
| 434 | & (fonc2n(icedgcm(1))-fonc2n(icedgcm(3))))+ |
|---|
| 435 | & ((pcapn_sm(n,2)-pcapn_sm(n,4))/ |
|---|
| 436 | & (fonc2n(icedgcm(2))-fonc2n(icedgcm(4)))))*(1./2.) |
|---|
| 437 | ! Evaluate derivative of pressure wrt Southern ice depth coefficient, |
|---|
| 438 | ! as the mean of the 2 values obtained using runs 1&3 and 2&4 |
|---|
| 439 | dpdes(n)=(((pcaps_sm(n,1)-pcaps_sm(n,3))/ |
|---|
| 440 | & (fonc2s(icedgcm(1))-fonc2s(icedgcm(3))))+ |
|---|
| 441 | & ((pcaps_sm(n,2)-pcaps_sm(n,4))/ |
|---|
| 442 | & (fonc2s(icedgcm(2))-fonc2s(icedgcm(4)))))*(1./2.) |
|---|
| 443 | ! Evaluate alphaVL1 coefficient, as the mean of alphaVL1 coefficients |
|---|
| 444 | ! of all 4 runs: |
|---|
| 445 | alphavl1(n) = (1./4.)*(pvl_sm(n,1)/patm_sm(n,1) |
|---|
| 446 | & +pvl_sm(n,2)/patm_sm(n,2)+pvl_sm(n,3)/patm_sm(n,3) |
|---|
| 447 | & +pvl_sm(n,4)/patm_sm(n,4)) |
|---|
| 448 | !write(91,*) sol(n), dpden(n), dpdes(n) |
|---|
| 449 | c write(*,*) 'pcapn_sm(n,1),pcapn_sm(n,2)' |
|---|
| 450 | c & ,pcapn_sm(n,1),pcapn_sm(n,2) |
|---|
| 451 | c write(*,*)'fonc(albgcmN(1)),fonc(albgcmS(2))', |
|---|
| 452 | c & fonc(albgcmN(1)),fonc(albgcmS(2)) |
|---|
| 453 | c write(*,*)'albgcmN(1)),albgcmS(2)', |
|---|
| 454 | c & albgcmN(1),albgcmS(2) |
|---|
| 455 | !write(90,*)pvl_sm(n,1)/patm_sm(n,1),pvl_sm(n,2)/patm_sm(n,2), |
|---|
| 456 | c & pvl_sm(n,3)/patm_sm(n,3),pvl_sm(n,4)/patm_sm(n,4) |
|---|
| 457 | end do ! of do n=1,nsol |
|---|
| 458 | |
|---|
| 459 | c ------------------------------------------------- |
|---|
| 460 | c Smooth the derivative like Frederic did ? : |
|---|
| 461 | call runave(sol,dpdicein,nsol,669.,60., |
|---|
| 462 | & sol,dpdicein_fltr,nsol) |
|---|
| 463 | call runave(sol,dpdiceis,nsol,669.,60., |
|---|
| 464 | & sol,dpdiceis_fltr,nsol) |
|---|
| 465 | call runave(sol,dpden,nsol,669.,60., |
|---|
| 466 | & sol,dpden_fltr,nsol) |
|---|
| 467 | call runave(sol,dpdes,nsol,669.,60., |
|---|
| 468 | & sol,dpdes_fltr,nsol) |
|---|
| 469 | do n=1,nsol |
|---|
| 470 | dpdiceis(n)=dpdiceis_fltr(n) |
|---|
| 471 | dpdicein(n)=dpdicein_fltr(n) |
|---|
| 472 | dpdes(n)=dpdes_fltr(n) |
|---|
| 473 | dpden(n)=dpden_fltr(n) |
|---|
| 474 | !write(92,*) sol(n), dpden_fltr(n), dpdes_fltr(n) |
|---|
| 475 | end do |
|---|
| 476 | c ------------------------------------------------- |
|---|
| 477 | |
|---|
| 478 | |
|---|
| 479 | c Compute best fit parameters by minimizing Cost function |
|---|
| 480 | c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
|---|
| 481 | c STUPID ultra-robust minimization method: waisting CPU to save human time... |
|---|
| 482 | |
|---|
| 483 | open(33, file = 'minimization.txt') |
|---|
| 484 | c 'minimzation.txt' will contain data to plot COST(DN,DS) or COST(IN,IS) depending on plot_test |
|---|
| 485 | |
|---|
| 486 | |
|---|
| 487 | c to plot COST(DN,DS) (plot_test=1), loops must be in this order : DN->DS->IN->IS |
|---|
| 488 | |
|---|
| 489 | costmin = 1.e30 ! initialization |
|---|
| 490 | |
|---|
| 491 | |
|---|
| 492 | if(plot_test.eq.1) then |
|---|
| 493 | c to plot COST(DN,DS) (plot_test=1), loops must be in this order : DN->DS->IN->IS |
|---|
| 494 | write(*,*)' DN DS cost IN IS Ps' |
|---|
| 495 | icedn=icedmin ! initialization |
|---|
| 496 | do while (icedn.le.icedmax) ! loop on northern ice depth coefficient |
|---|
| 497 | ! albn=albmin ! initialization |
|---|
| 498 | ! do while (albn.le.albmax) ! loop on northern cap albedo |
|---|
| 499 | iceds=max(icedn-maxiceddiff,icedmin) ! initialization |
|---|
| 500 | do while (iceds.le.min(icedn+maxiceddiff,icedmax)) |
|---|
| 501 | cost4plot = 1.e30 ! initializationqsub |
|---|
| 502 | ! iceds=max(icedn-maxiceddiff,icedmin) ! initialization |
|---|
| 503 | ! do while (iceds.le.min(icedn+maxiceddiff,icedmax)) |
|---|
| 504 | icein=iceimin ! initialization |
|---|
| 505 | do while (icein.le.iceimax) ! loop on northern ice inertia |
|---|
| 506 | iceis=max(icein-maxiceidiff,iceimin) ! initialization |
|---|
| 507 | do while (iceis.le.min(icein+maxiceidiff,iceimax)) |
|---|
| 508 | ptry=pmin ! initialization |
|---|
| 509 | do while (ptry.le.pmax) ! loop on total pressure |
|---|
| 510 | cost =0. !initialization |
|---|
| 511 | do n=1,nsol |
|---|
| 512 | !write(*,*) n |
|---|
| 513 | ! Pressure corresponding to Northern cap |
|---|
| 514 | pcapn_new=pcapn_sm(n,refrun)+ |
|---|
| 515 | & (fonc(icein) -fonc(iceigcm_refN)) * dpdicein(n) |
|---|
| 516 | & + (fonc2n(icedn) -fonc2n(icedgcm_ref)) * dpden(n) |
|---|
| 517 | pcapn_new= max(pcapn_new,0.) |
|---|
| 518 | ! Pressure corresponding to Southern Cap |
|---|
| 519 | pcaps_new=pcaps_sm(n,refrun)+ |
|---|
| 520 | & (fonc(iceis) -fonc(iceigcm_refS)) * dpdiceis(n) |
|---|
| 521 | & + (fonc2s(iceds) -fonc2s(icedgcm_ref)) * dpdes(n) |
|---|
| 522 | pcaps_new= max(pcaps_new,0.) |
|---|
| 523 | ! Pressure at VL1 site |
|---|
| 524 | pvl1 = alphavl1(n)* |
|---|
| 525 | & (ptry - pcapn_new - pcaps_new) |
|---|
| 526 | ! cumulative squared differences between predicted |
|---|
| 527 | ! and observed pressures at VL1 site |
|---|
| 528 | cost= cost+ ( pvl_obs(n) - pvl1)**2 |
|---|
| 529 | |
|---|
| 530 | end do ! of do n=1,nsol |
|---|
| 531 | ! store parameters which lead to minimum cost |
|---|
| 532 | ! (i.e. best fit so far) |
|---|
| 533 | if(cost.lt.costmin) then |
|---|
| 534 | costmin = cost |
|---|
| 535 | pfit=ptry |
|---|
| 536 | iceinfit=icein |
|---|
| 537 | iceisfit=iceis |
|---|
| 538 | icednfit=icedn |
|---|
| 539 | icedsfit=iceds |
|---|
| 540 | end if |
|---|
| 541 | |
|---|
| 542 | if(cost.lt.cost4plot) then |
|---|
| 543 | c RMS, best pressure and best albedos for these icedsivities value |
|---|
| 544 | cost4plot = cost |
|---|
| 545 | iceis4plot=iceis |
|---|
| 546 | ! iceds4plot=iceds |
|---|
| 547 | icein4plot=icein |
|---|
| 548 | ps4plot=ptry |
|---|
| 549 | end if |
|---|
| 550 | ptry=ptry+deltap ! increment ptry |
|---|
| 551 | end do ! of do while (ptry.le.pmax) |
|---|
| 552 | iceis=iceis+deltaicei ! increment iceis |
|---|
| 553 | end do ! of do while (iceis.le.min(icein+maxiceidiff,iceimax)) |
|---|
| 554 | icein=icein+deltaicei ! increment icein |
|---|
| 555 | enddo ! of do while (icein.le.iceimax) |
|---|
| 556 | ! iceds=iceds+deltaiced ! increment iceds |
|---|
| 557 | ! end do ! of do while (iceds.le.min(icedn+maxiceddiff,icedmax)) |
|---|
| 558 | ! write(*,fmt='(1pe9.2,f5.2,f9.3,1pe9.2, |
|---|
| 559 | ! & f5.2,f7.2)') |
|---|
| 560 | write(*,fmt='(6(1pe10.3))') ! output to screen |
|---|
| 561 | & icedn,iceds,sqrt(cost4plot/float(nsol)), |
|---|
| 562 | & icein4plot,iceis4plot,ps4plot |
|---|
| 563 | write(33,fmt='(6(1pe10.3))') ! output to file |
|---|
| 564 | & icedn,iceds,sqrt(cost4plot/float(nsol)), |
|---|
| 565 | & icein4plot,iceis4plot,ps4plot |
|---|
| 566 | iceds=iceds+deltaiced ! increment iceds |
|---|
| 567 | enddo ! of do while (iceds.le.min(icedn+maxiceddiff,icedmax)) |
|---|
| 568 | write(33,*)' ' ! blank line in output file |
|---|
| 569 | ! albn=albn+deltaalb ! increment albn |
|---|
| 570 | ! end do ! of do while (albn.le.albmax) |
|---|
| 571 | icedn=icedn+deltaiced |
|---|
| 572 | end do ! do while (icedn.le.icedmax) |
|---|
| 573 | |
|---|
| 574 | elseif(plot_test.eq.2) then |
|---|
| 575 | c to plot COST(IN,IS) (plot_test=2), loops must be in this order : IN->IS->DN->DS |
|---|
| 576 | write(*,*)' DN DS cost IN IS Ps' |
|---|
| 577 | write(33,*)' DN DS cost IN IS Ps' |
|---|
| 578 | icein=iceimin ! initialization |
|---|
| 579 | do while (icein.le.iceimax) ! loop on northern ice depth coefficient |
|---|
| 580 | ! albn=albmin ! initialization |
|---|
| 581 | ! do while (albn.le.albmax) ! loop on northern cap albedo |
|---|
| 582 | iceis=max(icein-maxiceidiff,iceimin)! initialization |
|---|
| 583 | do while (iceis.le.min(icein+maxiceidiff,iceimax)) |
|---|
| 584 | cost4plot = 1.e30 ! initializationqsub |
|---|
| 585 | ! iceds=max(icedn-maxiceddiff,icedmin) ! initialization |
|---|
| 586 | ! do while (iceds.le.min(icedn+maxiceddiff,icedmax)) |
|---|
| 587 | icedn=icedmin ! initialization |
|---|
| 588 | do while (icedn.le.icedmax) ! loop on northern ice inertia |
|---|
| 589 | iceds=max(icedn-maxiceddiff,icedmin) ! initialization |
|---|
| 590 | do while (iceds.le.min(icedn+maxiceddiff,icedmax)) |
|---|
| 591 | ptry=pmin ! initialization |
|---|
| 592 | do while (ptry.le.pmax) ! loop on total pressure |
|---|
| 593 | cost =0. !initialization |
|---|
| 594 | do n=1,nsol |
|---|
| 595 | !write(*,*) n |
|---|
| 596 | ! Pressure corresponding to Northern cap |
|---|
| 597 | pcapn_new=pcapn_sm(n,refrun)+ |
|---|
| 598 | & (fonc(icein) -fonc(iceigcm_refN)) * dpdicein(n) |
|---|
| 599 | & + (fonc2n(icedn) -fonc2n(icedgcm_ref)) * dpden(n) |
|---|
| 600 | pcapn_new= max(pcapn_new,0.) |
|---|
| 601 | ! Pressure corresponding to Southern Cap |
|---|
| 602 | pcaps_new=pcaps_sm(n,refrun)+ |
|---|
| 603 | & (fonc(iceis) -fonc(iceigcm_refS)) * dpdiceis(n) |
|---|
| 604 | & + (fonc2s(iceds) -fonc2s(icedgcm_ref)) * dpdes(n) |
|---|
| 605 | pcaps_new= max(pcaps_new,0.) |
|---|
| 606 | ! Pressure at VL1 site |
|---|
| 607 | pvl1 = alphavl1(n)* |
|---|
| 608 | & (ptry - pcapn_new - pcaps_new) |
|---|
| 609 | ! cumulative squared differences between predicted |
|---|
| 610 | ! and observed pressures at VL1 site |
|---|
| 611 | cost= cost+ ( pvl_obs(n) - pvl1)**2 |
|---|
| 612 | |
|---|
| 613 | end do ! of do n=1,nsol |
|---|
| 614 | ! store parameters which lead to minimum cost |
|---|
| 615 | ! (i.e. best fit so far) |
|---|
| 616 | if(cost.lt.costmin) then |
|---|
| 617 | costmin = cost |
|---|
| 618 | pfit=ptry |
|---|
| 619 | iceinfit=icein |
|---|
| 620 | iceisfit=iceis |
|---|
| 621 | icednfit=icedn |
|---|
| 622 | icedsfit=iceds |
|---|
| 623 | end if |
|---|
| 624 | |
|---|
| 625 | if(cost.lt.cost4plot) then |
|---|
| 626 | c RMS, best pressure and best albedos for these icedsivities value |
|---|
| 627 | cost4plot = cost |
|---|
| 628 | iceds4plot=iceds |
|---|
| 629 | ! iceds4plot=iceds |
|---|
| 630 | icedn4plot=icedn |
|---|
| 631 | ps4plot=ptry |
|---|
| 632 | end if |
|---|
| 633 | ptry=ptry+deltap ! increment ptry |
|---|
| 634 | end do ! of do while (ptry.le.pmax) |
|---|
| 635 | iceds=iceds+deltaiced ! increment iceis |
|---|
| 636 | end do ! of do while (iceis.le.min(icein+maxiceidiff,iceimax)) |
|---|
| 637 | icedn=icedn+deltaiced ! increment icein |
|---|
| 638 | enddo ! of do while (icein.le.iceimax) |
|---|
| 639 | ! iceds=iceds+deltaiced ! increment iceds |
|---|
| 640 | ! end do ! of do while (iceds.le.min(icedn+maxiceddiff,icedmax)) |
|---|
| 641 | ! write(*,fmt='(1pe9.2,f5.2,f9.3,1pe9.2, |
|---|
| 642 | ! & f5.2,f7.2)') |
|---|
| 643 | write(*,fmt='(6(1pe10.3))') ! output to screen |
|---|
| 644 | & icedn4plot,iceds4plot,sqrt(cost4plot/float(nsol)), |
|---|
| 645 | & icein,iceis,ps4plot |
|---|
| 646 | write(33,fmt='(6(1pe10.3))') ! output to file |
|---|
| 647 | & icedn4plot,iceds4plot,sqrt(cost4plot/float(nsol)), |
|---|
| 648 | & icein,iceis,ps4plot |
|---|
| 649 | iceis=iceis+deltaicei ! increment iceds |
|---|
| 650 | enddo ! of do while (iceds.le.min(icedn+maxiceddiff,icedmax)) |
|---|
| 651 | write(33,*)' ' ! blank line in output file |
|---|
| 652 | ! albn=albn+deltaalb ! increment albn |
|---|
| 653 | ! end do ! of do while (albn.le.albmax) |
|---|
| 654 | icein=icein+deltaicei |
|---|
| 655 | end do ! do while (icedn.le.icedmax) |
|---|
| 656 | else |
|---|
| 657 | write(*,*) 'Wrong integer for plot (must be 1 or 2)' |
|---|
| 658 | stop |
|---|
| 659 | |
|---|
| 660 | endif |
|---|
| 661 | |
|---|
| 662 | close(33) |
|---|
| 663 | |
|---|
| 664 | write(*,*) 'Best fit Ptot=', pfit |
|---|
| 665 | write(*,*) 'Best fit In=', iceinfit |
|---|
| 666 | write(*,*) 'Best fit IS=', iceisfit |
|---|
| 667 | write(*,*) 'Best fit Dn=', icednfit |
|---|
| 668 | write(*,*) 'Best fit DS=', icedsfit |
|---|
| 669 | write(*,*) 'RMS difference model/obs=',sqrt(costmin/float(nsol)) |
|---|
| 670 | |
|---|
| 671 | c Synthethic VL1 pressure curves for information |
|---|
| 672 | c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
|---|
| 673 | open(41,file='pvlfit') |
|---|
| 674 | open(42,file='xprestotfit') |
|---|
| 675 | do n=1,nsol |
|---|
| 676 | pcapn_new=pcapn_sm(n,refrun)+ |
|---|
| 677 | & (fonc(iceinfit) -fonc(iceigcm_refN)) * dpdicein(n) |
|---|
| 678 | & + (fonc2n(icednfit) -fonc2n(icedgcm_ref)) * dpden(n) |
|---|
| 679 | pcapn_new= max(pcapn_new,0.) |
|---|
| 680 | pcaps_new=pcaps_sm(n,refrun)+ |
|---|
| 681 | & (fonc(iceisfit) -fonc(iceigcm_refS)) * dpdiceis(n) |
|---|
| 682 | & + (fonc2s(icedsfit) -fonc2s(icedgcm_ref)) * dpdes(n) |
|---|
| 683 | pcaps_new= max(pcaps_new,0.) |
|---|
| 684 | |
|---|
| 685 | call sol2ls(sol(n),solconv) |
|---|
| 686 | |
|---|
| 687 | write(41,*) sol(n), solconv, alphavl1(n)* |
|---|
| 688 | & ( pfit - pcapn_new - pcaps_new) |
|---|
| 689 | write(42,*) sol(n),solconv, pfit - pcapn_new - pcaps_new, |
|---|
| 690 | & pcapn_new, pcaps_new, pfit |
|---|
| 691 | end do |
|---|
| 692 | close(41) |
|---|
| 693 | close(42) |
|---|
| 694 | |
|---|
| 695 | end |
|---|
| 696 | |
|---|
| 697 | c ***************************************************************** |
|---|
| 698 | SUBROUTINE RUNAVE (x,y,nmax,xperiod,xave,xsmooth,ysmooth,nsm) |
|---|
| 699 | |
|---|
| 700 | |
|---|
| 701 | IMPLICIT NONE |
|---|
| 702 | |
|---|
| 703 | c ----------------------------------------------------------- |
|---|
| 704 | c Computing runnig average ysmooth on xsmooth coordinate for a periodic |
|---|
| 705 | c field y in coordinate x PERIODIC between 0 and xperiod |
|---|
| 706 | c Averaging box size : xave |
|---|
| 707 | c F.Forget 1999 |
|---|
| 708 | c ----------------------------------------------------------- |
|---|
| 709 | |
|---|
| 710 | integer nmax,nsm |
|---|
| 711 | real y(nmax), ysmooth(nsm) |
|---|
| 712 | real x(nmax), xsmooth(nsm) |
|---|
| 713 | real xperiod,xave |
|---|
| 714 | integer n,i, nave, imin |
|---|
| 715 | |
|---|
| 716 | integer nbig |
|---|
| 717 | parameter (nbig=99999999) |
|---|
| 718 | real xx(nbig) |
|---|
| 719 | |
|---|
| 720 | c ----------------------------------------------------------- |
|---|
| 721 | |
|---|
| 722 | |
|---|
| 723 | if (nbig.lt.3*nmax) then |
|---|
| 724 | write(*,*) 'Must increase nbig in runave' |
|---|
| 725 | stop |
|---|
| 726 | endif |
|---|
| 727 | |
|---|
| 728 | c Reindexation des donnees |
|---|
| 729 | do n=1,nmax |
|---|
| 730 | xx(n) = x(n) - xperiod |
|---|
| 731 | xx(n+nmax) = x(n) |
|---|
| 732 | xx(n+2*nmax) = x(n) +xperiod |
|---|
| 733 | end do |
|---|
| 734 | |
|---|
| 735 | c Moyenne glissante |
|---|
| 736 | imin=1 |
|---|
| 737 | do n=1,nsm |
|---|
| 738 | ysmooth(n) =0. |
|---|
| 739 | nave =0 |
|---|
| 740 | do i=imin,3*nmax |
|---|
| 741 | if (xx(i).ge.(xsmooth(n)-0.5*xave)) then |
|---|
| 742 | if (xx(i).gt.(xsmooth(n)+0.5*xave)) goto 999 |
|---|
| 743 | ysmooth(n) = ysmooth(n) + y(mod(i-1,nmax)+1) |
|---|
| 744 | nave = nave +1 |
|---|
| 745 | end if |
|---|
| 746 | end do |
|---|
| 747 | 999 continue |
|---|
| 748 | imin = i -1 -nave |
|---|
| 749 | ysmooth(n) = ysmooth(n)/ float (nave) |
|---|
| 750 | end do |
|---|
| 751 | end |
|---|
| 752 | |
|---|
| 753 | subroutine ls2sol(ls,sol) |
|---|
| 754 | |
|---|
| 755 | implicit none |
|---|
| 756 | !================================================================ |
|---|
| 757 | ! Arguments: |
|---|
| 758 | !================================================================ |
|---|
| 759 | real,intent(in) :: ls |
|---|
| 760 | real,intent(out) :: sol |
|---|
| 761 | |
|---|
| 762 | !================================================================ |
|---|
| 763 | ! Local: |
|---|
| 764 | !================================================================ |
|---|
| 765 | double precision xref,zx0,zteta,zz |
|---|
| 766 | !xref: mean anomaly, zteta: true anomaly, zx0: eccentric anomaly |
|---|
| 767 | double precision year_day |
|---|
| 768 | double precision peri_day,timeperi,e_elips |
|---|
| 769 | double precision pi,degrad |
|---|
| 770 | parameter (year_day=668.6d0) ! number of sols in a martian year |
|---|
| 771 | parameter (peri_day=485.35d0) ! date (in sols) of perihelion |
|---|
| 772 | !timeperi: 2*pi*( 1 - Ls(perihelion)/ 360 ); Ls(perihelion)=250.99 |
|---|
| 773 | parameter (timeperi=1.90258341759902d0) |
|---|
| 774 | parameter (e_elips=0.0934d0) ! eccentricity of orbit |
|---|
| 775 | parameter (pi=3.14159265358979d0) |
|---|
| 776 | parameter (degrad=57.2957795130823d0) |
|---|
| 777 | |
|---|
| 778 | if (abs(ls).lt.1.0e-5) then |
|---|
| 779 | if (ls.ge.0.0) then |
|---|
| 780 | sol = 0.0 |
|---|
| 781 | else |
|---|
| 782 | sol = real(year_day) |
|---|
| 783 | end if |
|---|
| 784 | return |
|---|
| 785 | end if |
|---|
| 786 | |
|---|
| 787 | zteta = ls/degrad + timeperi |
|---|
| 788 | zx0 = 2.0*datan(dtan(0.5*zteta)/dsqrt((1.+e_elips)/(1.-e_elips))) |
|---|
| 789 | xref = zx0-e_elips*dsin(zx0) |
|---|
| 790 | zz = xref/(2.*pi) |
|---|
| 791 | sol = real(zz*year_day + peri_day) |
|---|
| 792 | if (sol.lt.0.0) sol = sol + real(year_day) |
|---|
| 793 | if (sol.ge.year_day) sol = sol - real(year_day) |
|---|
| 794 | |
|---|
| 795 | |
|---|
| 796 | end subroutine ls2sol |
|---|
| 797 | |
|---|
| 798 | subroutine sol2ls(sol,Ls) |
|---|
| 799 | !============================================================================== |
|---|
| 800 | ! Purpose: |
|---|
| 801 | ! Convert a date/time, given in sol (martian day), |
|---|
| 802 | ! into solar longitude date/time, in Ls (in degrees), |
|---|
| 803 | ! where sol=0 is (by definition) the northern hemisphere |
|---|
| 804 | ! spring equinox (where Ls=0). |
|---|
| 805 | !============================================================================== |
|---|
| 806 | ! Notes: |
|---|
| 807 | ! Even though "Ls" is cyclic, if "sol" is greater than N (martian) year, |
|---|
| 808 | ! "Ls" will be increased by N*360 |
|---|
| 809 | ! Won't work as expected if sol is negative (then again, |
|---|
| 810 | ! why would that ever happen?) |
|---|
| 811 | !============================================================================== |
|---|
| 812 | |
|---|
| 813 | implicit none |
|---|
| 814 | |
|---|
| 815 | !============================================================================== |
|---|
| 816 | ! Arguments: |
|---|
| 817 | !============================================================================== |
|---|
| 818 | real,intent(in) :: sol |
|---|
| 819 | real,intent(out) :: Ls |
|---|
| 820 | |
|---|
| 821 | !============================================================================== |
|---|
| 822 | ! Local variables: |
|---|
| 823 | !============================================================================== |
|---|
| 824 | real year_day,peri_day,timeperi,e_elips,twopi,degrad |
|---|
| 825 | data year_day /669./ ! # of sols in a martian year |
|---|
| 826 | data peri_day /485.0/ |
|---|
| 827 | data timeperi /1.9082314/ |
|---|
| 828 | data e_elips /0.093358/ |
|---|
| 829 | data twopi /6.2831853/ ! 2.*pi |
|---|
| 830 | data degrad /57.2957795/ ! pi/180 |
|---|
| 831 | |
|---|
| 832 | real zanom,xref,zx0,zdx,zteta,zz |
|---|
| 833 | |
|---|
| 834 | integer count_years |
|---|
| 835 | integer iter |
|---|
| 836 | |
|---|
| 837 | !============================================================================== |
|---|
| 838 | ! 1. Compute Ls |
|---|
| 839 | !============================================================================== |
|---|
| 840 | |
|---|
| 841 | zz=(sol-peri_day)/year_day |
|---|
| 842 | zanom=twopi*(zz-nint(zz)) |
|---|
| 843 | xref=abs(zanom) |
|---|
| 844 | |
|---|
| 845 | ! The equation zx0 - e * sin (zx0) = xref, solved by Newton |
|---|
| 846 | zx0=xref+e_elips*sin(xref) |
|---|
| 847 | do iter=1,20 ! typically, 2 or 3 iterations are enough |
|---|
| 848 | zdx=-(zx0-e_elips*sin(zx0)-xref)/(1.-e_elips*cos(zx0)) |
|---|
| 849 | zx0=zx0+zdx |
|---|
| 850 | if(abs(zdx).le.(1.e-7)) then |
|---|
| 851 | ! write(*,*)'iter:',iter,' |zdx|:',abs(zdx) |
|---|
| 852 | exit |
|---|
| 853 | endif |
|---|
| 854 | enddo |
|---|
| 855 | |
|---|
| 856 | if(zanom.lt.0.) zx0=-zx0 |
|---|
| 857 | |
|---|
| 858 | zteta=2.*atan(sqrt((1.+e_elips)/(1.-e_elips))*tan(zx0/2.)) |
|---|
| 859 | Ls=zteta-timeperi |
|---|
| 860 | |
|---|
| 861 | if(Ls.lt.0.) then |
|---|
| 862 | Ls=Ls+twopi |
|---|
| 863 | else |
|---|
| 864 | if(Ls.gt.twopi) then |
|---|
| 865 | Ls=Ls-twopi |
|---|
| 866 | endif |
|---|
| 867 | endif |
|---|
| 868 | |
|---|
| 869 | Ls=degrad*Ls |
|---|
| 870 | ! Ls is now in degrees |
|---|
| 871 | |
|---|
| 872 | !============================================================================== |
|---|
| 873 | ! 1. Account for (eventual) years included in input date/time sol |
|---|
| 874 | !============================================================================== |
|---|
| 875 | |
|---|
| 876 | count_years=0 ! initialize |
|---|
| 877 | zz=sol ! use "zz" to store (and work on) the value of sol |
|---|
| 878 | do while (zz.ge.year_day) |
|---|
| 879 | count_years=count_years+1 |
|---|
| 880 | zz=zz-year_day |
|---|
| 881 | enddo |
|---|
| 882 | |
|---|
| 883 | ! Add 360 degrees to Ls for every year |
|---|
| 884 | if (count_years.ne.0) then |
|---|
| 885 | Ls=Ls+360.*count_years |
|---|
| 886 | endif |
|---|
| 887 | |
|---|
| 888 | |
|---|
| 889 | end subroutine sol2ls |
|---|
| 890 | |
|---|
