| 1 | SUBROUTINE surfini(ngrid,piceco2,qsurf,psolaralb) |
|---|
| 2 | ! to use 'getin' |
|---|
| 3 | USE ioipsl_getincom |
|---|
| 4 | IMPLICIT NONE |
|---|
| 5 | c======================================================================= |
|---|
| 6 | c |
|---|
| 7 | c creation des calottes pour l'etat initial |
|---|
| 8 | c |
|---|
| 9 | c======================================================================= |
|---|
| 10 | c----------------------------------------------------------------------- |
|---|
| 11 | c Declarations: |
|---|
| 12 | c ------------- |
|---|
| 13 | #include "dimensions.h" |
|---|
| 14 | #include "dimphys.h" |
|---|
| 15 | #include "surfdat.h" |
|---|
| 16 | #include "callkeys.h" |
|---|
| 17 | #include "tracer.h" |
|---|
| 18 | #include "comgeomfi.h" |
|---|
| 19 | #include "comcstfi.h" |
|---|
| 20 | c |
|---|
| 21 | INTEGER ngrid,ig,icap,iq,alternate |
|---|
| 22 | REAL piceco2(ngrid),psolaralb(ngrid,2) |
|---|
| 23 | REAL qsurf(ngrid,nqmx) !tracer on surface (kg/m2) |
|---|
| 24 | REAL icedryness ! ice dryness |
|---|
| 25 | |
|---|
| 26 | ! longwatercaptag is watercaptag. Trick for some compilers |
|---|
| 27 | LOGICAL, DIMENSION(100000) :: longwatercaptag |
|---|
| 28 | |
|---|
| 29 | EXTERNAL ISMIN,ISMAX |
|---|
| 30 | INTEGER ISMIN,ISMAX |
|---|
| 31 | |
|---|
| 32 | ! Choose false to have a somewhat non resolution dependant water ice caps distribution, |
|---|
| 33 | ! i.e based only on lat & lon values of each physical point. |
|---|
| 34 | ! Choose true to get a carefully choosen distribution for GCM resolutions 32x24 or 64x48 |
|---|
| 35 | ! For vizualisation : /u/tnalmd/bin/watercaps gcm_txt_output |
|---|
| 36 | LOGICAL,SAVE :: improvedicelocation = .true. |
|---|
| 37 | c |
|---|
| 38 | c======================================================================= |
|---|
| 39 | |
|---|
| 40 | c water ice outliers |
|---|
| 41 | c ------------------------------------------ |
|---|
| 42 | |
|---|
| 43 | IF ((water) .and. (caps)) THEN |
|---|
| 44 | |
|---|
| 45 | c Perennial H20 north cap defined by watercaptag=true (allows surface to be |
|---|
| 46 | c hollowed by sublimation in vdifc). |
|---|
| 47 | |
|---|
| 48 | c We might not want albedodat to be modified because it is used to write |
|---|
| 49 | c restart files. Instead, albedo is directly modified when needed (i.e. |
|---|
| 50 | c if we have watercaptag and no co2 ice), below and in albedocaps.F90 |
|---|
| 51 | |
|---|
| 52 | c "Dryness coefficient" controlling the evaporation and |
|---|
| 53 | c sublimation from the ground water ice (close to 1) |
|---|
| 54 | c HERE, the goal is to correct for the fact |
|---|
| 55 | c that the simulated permanent water ice polar caps |
|---|
| 56 | c is larger than the actual cap and the atmospheric |
|---|
| 57 | c opacity not always realistic. |
|---|
| 58 | |
|---|
| 59 | alternate = 0 |
|---|
| 60 | |
|---|
| 61 | if (ngridmx .ne. 1) then |
|---|
| 62 | watercaptag(:) = .false. |
|---|
| 63 | longwatercaptag(:) = .false. |
|---|
| 64 | endif |
|---|
| 65 | |
|---|
| 66 | write(*,*) "Ice dryness ?" |
|---|
| 67 | icedryness=1. ! default value |
|---|
| 68 | call getin("icedryness",icedryness) |
|---|
| 69 | write(*,*) " icedryness = ",icedryness |
|---|
| 70 | |
|---|
| 71 | |
|---|
| 72 | #ifdef MESOSCALE |
|---|
| 73 | |
|---|
| 74 | do ig=1,ngridmx |
|---|
| 75 | |
|---|
| 76 | !write(*,*) "all qsurf to zero. dirty." |
|---|
| 77 | do iq=1,nqmx |
|---|
| 78 | qsurf(ig,iq)=0. !! on jette les inputs GCM |
|---|
| 79 | !! on regle juste watercaptag |
|---|
| 80 | !! il faudrait garder les inputs GCM |
|---|
| 81 | !! si elles sont consequentes |
|---|
| 82 | enddo |
|---|
| 83 | if ( ( lati(ig)*180./pi .gt. 70. ) .and. |
|---|
| 84 | . ( albedodat(ig) .ge. 0.26 ) ) then |
|---|
| 85 | write(*,*)"outlier ",ig |
|---|
| 86 | watercaptag(ig) = .true. |
|---|
| 87 | dryness(ig) = 1. |
|---|
| 88 | albedodat(ig) = albedo_h2o_ice !! pour output |
|---|
| 89 | else |
|---|
| 90 | watercaptag(ig) = .false. |
|---|
| 91 | dryness(ig) = 1. |
|---|
| 92 | endif |
|---|
| 93 | |
|---|
| 94 | enddo |
|---|
| 95 | |
|---|
| 96 | #else |
|---|
| 97 | |
|---|
| 98 | if (improvedicelocation) then |
|---|
| 99 | |
|---|
| 100 | if (ngridmx .eq. 738) then ! hopefully 32x24 |
|---|
| 101 | |
|---|
| 102 | print*,'water ice caps distribution for 32x24 resolution:' |
|---|
| 103 | longwatercaptag(1:9) = .true. ! central cap - core |
|---|
| 104 | longwatercaptag(26:33) = .true. ! central cap |
|---|
| 105 | |
|---|
| 106 | else if (ngridmx .eq. 3010) then ! hopefully 64x48 |
|---|
| 107 | |
|---|
| 108 | ! For latitudes: |
|---|
| 109 | ! [ 67.5 71.25 75. 78.75 82.5 86.25] |
|---|
| 110 | ! The water ice caps should be (according to TES temperatures): |
|---|
| 111 | ! [ 8.63e-03 1.02e+00 5.99e+00 2.66e+01 5.65e+01] |
|---|
| 112 | |
|---|
| 113 | print*,'water ice caps distribution for 64x48 resolution:' |
|---|
| 114 | longwatercaptag(1:65) = .true. ! central cap - core |
|---|
| 115 | longwatercaptag(75:85) = .true. ! central cap |
|---|
| 116 | longwatercaptag(93:114) = .true. ! central cap |
|---|
| 117 | !--------------------- OUTLIERS ---------------------------- |
|---|
| 118 | longwatercaptag(136) = .true. ! outlier, lat = 78.75 |
|---|
| 119 | longwatercaptag(138) = .true. ! outlier, lat = 78.75 |
|---|
| 120 | longwatercaptag(140) = .true. ! outlier, lat = 78.75 |
|---|
| 121 | longwatercaptag(142) = .true. ! outlier, lat = 78.75 |
|---|
| 122 | longwatercaptag(161) = .true. ! outlier, lat = 78.75 |
|---|
| 123 | longwatercaptag(163) = .true. ! outlier, lat = 78.75 |
|---|
| 124 | longwatercaptag(165) = .true. ! outlier, lat = 78.75 |
|---|
| 125 | longwatercaptag(183) = .true. ! outlier, lat = 78.75 |
|---|
| 126 | longwatercaptag(185) = .true. ! outlier, lat = 78.75 |
|---|
| 127 | longwatercaptag(187) = .true. ! outlier, lat = 78.75 |
|---|
| 128 | longwatercaptag(189) = .true. ! outlier, lat = 78.75 |
|---|
| 129 | longwatercaptag(191) = .true. ! outlier, lat = 78.75 |
|---|
| 130 | longwatercaptag(193) = .true. ! outlier, lat = 78.75 |
|---|
| 131 | longwatercaptag(194) = .true. ! outlier, lat = 75 |
|---|
| 132 | longwatercaptag(203) = .true. ! outlier, lat = 75 |
|---|
| 133 | longwatercaptag(207) = .true. ! outlier, lat = 75 |
|---|
| 134 | longwatercaptag(242) = .true. ! outlier, lat = 75 |
|---|
| 135 | longwatercaptag(244) = .true. ! outlier, lat = 75 |
|---|
| 136 | longwatercaptag(246) = .true. ! outlier, lat = 75 |
|---|
| 137 | longwatercaptag(248) = .true. ! outlier, lat = 75 |
|---|
| 138 | longwatercaptag(250) = .true. ! outlier, lat = 75 |
|---|
| 139 | longwatercaptag(252) = .true. ! outlier, lat = 75 |
|---|
| 140 | longwatercaptag(254) = .true. ! outlier, lat = 75 |
|---|
| 141 | longwatercaptag(256:257)= .true. ! outlier, lat = 75 |
|---|
| 142 | !!------------------------- OLD ---------------------------- |
|---|
| 143 | !! longwatercaptag(83:85) = .true. |
|---|
| 144 | !! longwatercaptag(135:142) = .true. |
|---|
| 145 | !! longwatercaptag(181:193) = .true. |
|---|
| 146 | !! longwatercaptag(242:257) = .true. |
|---|
| 147 | |
|---|
| 148 | |
|---|
| 149 | else if (ngridmx .ne. 1) then |
|---|
| 150 | print*,'No improved ice location for this resolution!' |
|---|
| 151 | print*,'Please set improvedicelocation to false in surfini.' |
|---|
| 152 | print*,'And check lat and lon conditions for ice caps in code.' |
|---|
| 153 | call exit(1) |
|---|
| 154 | |
|---|
| 155 | endif |
|---|
| 156 | |
|---|
| 157 | ! print caps locations |
|---|
| 158 | print*,'latitude | longitude | ig' |
|---|
| 159 | do ig=1,ngridmx |
|---|
| 160 | dryness (ig) = icedryness |
|---|
| 161 | |
|---|
| 162 | if (longwatercaptag(ig)) then |
|---|
| 163 | watercaptag(ig) = .True. |
|---|
| 164 | print*,'ice water cap', lati(ig)*180./pi, |
|---|
| 165 | . long(ig)*180./pi, ig |
|---|
| 166 | endif |
|---|
| 167 | enddo |
|---|
| 168 | |
|---|
| 169 | |
|---|
| 170 | else |
|---|
| 171 | |
|---|
| 172 | do ig=1,ngridmx |
|---|
| 173 | |
|---|
| 174 | dryness (ig) = icedryness |
|---|
| 175 | |
|---|
| 176 | !!c Towards olympia planitia water caps ('relatively' low latitude ones) |
|---|
| 177 | !!c---------------- proposition par AS -------------------- |
|---|
| 178 | !!c-------------------------------------------------------- |
|---|
| 179 | !!c if ( ( lati(ig)*180./pi .ge. 75 ) .and. |
|---|
| 180 | !!c . ( lati(ig)*180./pi .le. 77 ) .and. |
|---|
| 181 | !!c . ( ( ( long(ig)*180./pi .ge. 000. ) .and. |
|---|
| 182 | !!c . ( long(ig)*180./pi .le. 120. ) ) |
|---|
| 183 | !!c . .or. |
|---|
| 184 | !!c . ( ( long(ig)*180./pi .ge. -130. ) .and. |
|---|
| 185 | !!c . ( long(ig)*180./pi .le. -115. ) ) ) ) then |
|---|
| 186 | !!c---------------- proposition par TN -------------------- |
|---|
| 187 | !!c---------------- HIGHLY EXPERIMENTAL ------------------- |
|---|
| 188 | !!c-------------------------------------------------------- |
|---|
| 189 | !! if ( ( ( lati(ig)*180./pi .ge. 73. ) .and. ! cap #1 |
|---|
| 190 | !! . ( lati(ig)*180./pi .le. 75.1 ) .and. |
|---|
| 191 | !! . ( long(ig)*180./pi .ge. 95. ) .and. |
|---|
| 192 | !! . ( long(ig)*180./pi .le. 110. ) ) |
|---|
| 193 | !! . .or. |
|---|
| 194 | !! . ( ( lati(ig)*180./pi .ge. 77. ) .and. ! cap #2 |
|---|
| 195 | !! . ( lati(ig)*180./pi .le. 80. ) .and. |
|---|
| 196 | !! . ( long(ig)*180./pi .ge. 110. ) .and. |
|---|
| 197 | !! . ( long(ig)*180./pi .le. 140. ) ) |
|---|
| 198 | !! . .or. |
|---|
| 199 | !! . ( ( lati(ig)*180./pi .ge. 74.9 ) .and. ! cap #3 |
|---|
| 200 | !! . ( lati(ig)*180./pi .le. 78. ) .and. |
|---|
| 201 | !! . ( long(ig)*180./pi .ge. 155. ) .and. |
|---|
| 202 | !! . ( long(ig)*180./pi .le. 180. ) ) |
|---|
| 203 | !! . .or. |
|---|
| 204 | !! . ( ( lati(ig)*180./pi .ge. 71. ) .and. ! cap #4 (Korolev crater) |
|---|
| 205 | !! . ( lati(ig)*180./pi .le. 72. ) .and. |
|---|
| 206 | !! . ( long(ig)*180./pi .ge. 163. ) .and. |
|---|
| 207 | !! . ( long(ig)*180./pi .le. 168. ) ) |
|---|
| 208 | !! . .or. |
|---|
| 209 | !! . ( ( lati(ig)*180./pi .ge. 74.9 ) .and. ! cap #5 |
|---|
| 210 | !! . ( lati(ig)*180./pi .le. 78. ) .and. |
|---|
| 211 | !! . ( long(ig)*180./pi .ge. -160.) .and. |
|---|
| 212 | !! . ( long(ig)*180./pi .le. -120.) ) ) |
|---|
| 213 | !! . then |
|---|
| 214 | if ( ( ( lati(ig)*180./pi .ge. 77. ) .and. ! cap #2 |
|---|
| 215 | . ( lati(ig)*180./pi .le. 80. ) .and. |
|---|
| 216 | . ( long(ig)*180./pi .ge. 110. ) .and. |
|---|
| 217 | . ( long(ig)*180./pi .le. 181. ) ) |
|---|
| 218 | . .or. |
|---|
| 219 | |
|---|
| 220 | . ( ( lati(ig)*180./pi .ge. 75. ) .and. ! cap #4 (Korolev crater) |
|---|
| 221 | . ( lati(ig)*180./pi .le. 76. ) .and. |
|---|
| 222 | . ( long(ig)*180./pi .ge. 150. ) .and. |
|---|
| 223 | . ( long(ig)*180./pi .le. 168. ) ) |
|---|
| 224 | . .or. |
|---|
| 225 | . ( ( lati(ig)*180./pi .ge. 77 ) .and. ! cap #5 |
|---|
| 226 | . ( lati(ig)*180./pi .le. 80. ) .and. |
|---|
| 227 | . ( long(ig)*180./pi .ge. -150.) .and. |
|---|
| 228 | . ( long(ig)*180./pi .le. -110.) ) ) |
|---|
| 229 | . then |
|---|
| 230 | |
|---|
| 231 | |
|---|
| 232 | if ((alternate .eq. 0)) then ! 1/2 en 64x48 sinon trop large en lat |
|---|
| 233 | if (ngridmx.ne.1) watercaptag(ig)=.true. |
|---|
| 234 | print*,'ice water cap', lati(ig)*180./pi, |
|---|
| 235 | . long(ig)*180./pi, ig |
|---|
| 236 | dryness(ig) = 1. |
|---|
| 237 | alternate = 1 |
|---|
| 238 | else |
|---|
| 239 | alternate = 0 |
|---|
| 240 | endif !end if alternate = 0 |
|---|
| 241 | |
|---|
| 242 | |
|---|
| 243 | endif |
|---|
| 244 | |
|---|
| 245 | |
|---|
| 246 | c Opposite olympia planitia water cap |
|---|
| 247 | c---------------- proposition par AS -------------------- |
|---|
| 248 | c-------------------------------------------------------- |
|---|
| 249 | c if ( ( lati(ig)*180./pi .ge. 82 ) .and. |
|---|
| 250 | c . ( lati(ig)*180./pi .le. 84 ) .and. |
|---|
| 251 | c . ( ( long(ig)*180./pi .gt. -030. ) .and. |
|---|
| 252 | c . ( long(ig)*180./pi .lt. 090. ) ) ) then |
|---|
| 253 | c---------------- proposition par TN -------------------- |
|---|
| 254 | c-------------------------------------------------------- |
|---|
| 255 | if ( ( ( lati(ig)*180./pi .ge. 80 ) .and. |
|---|
| 256 | . ( lati(ig)*180./pi .le. 84 ) ) |
|---|
| 257 | . .and. |
|---|
| 258 | . ( ( long(ig)*180./pi .lt. -95. ) .or. !!! 32x24 |
|---|
| 259 | . ( long(ig)*180./pi .gt. 85. ) ) ) then !!! 32x24 |
|---|
| 260 | ! . ( ( ( long(ig)*180./pi .ge. -29. ) .and. !!! 64x48 |
|---|
| 261 | ! . ( long(ig)*180./pi .le. 90. ) ) .or. !!! 64x48 |
|---|
| 262 | ! . ( ( long(ig)*180./pi .ge. -77. ) .and. !!! 64x48 |
|---|
| 263 | ! . ( long(ig)*180./pi .le. -70. ) ) ) ) then !!! 64x48 |
|---|
| 264 | if (ngridmx.ne.1) then |
|---|
| 265 | watercaptag(ig)=.true. |
|---|
| 266 | print*,'ice water cap', lati(ig)*180./pi, |
|---|
| 267 | . long(ig)*180./pi, ig |
|---|
| 268 | endif |
|---|
| 269 | dryness(ig) = 1. |
|---|
| 270 | endif |
|---|
| 271 | |
|---|
| 272 | c Central cap |
|---|
| 273 | c---------------- anciens reglages FF ------------------- |
|---|
| 274 | c-------------------------------------------------------- |
|---|
| 275 | |
|---|
| 276 | if (abs(lati(ig)*180./pi).gt.80) then |
|---|
| 277 | print*,'ice water cap', lati(ig)*180./pi, |
|---|
| 278 | . long(ig)*180./pi, ig |
|---|
| 279 | if (ngridmx.ne.1) watercaptag(ig)=.true. |
|---|
| 280 | !dryness(ig) = 1. |
|---|
| 281 | c Use the following cap definition for high spatial resolution (latitudinal bin <= 5 deg) |
|---|
| 282 | c if (lati(ig)*180./pi.lt.85.and.long(ig).ge.0) then |
|---|
| 283 | c if (ngridmx.ne.1) watercaptag(ig)=.true. |
|---|
| 284 | c dryness(ig) = 1. |
|---|
| 285 | c endif |
|---|
| 286 | c if (lati(ig)*180./pi.ge.85) then |
|---|
| 287 | c if (ngridmx.ne.1) watercaptag(ig)=.true. |
|---|
| 288 | c dryness(ig) = 1. |
|---|
| 289 | c endif |
|---|
| 290 | endif ! (lati>80 deg) |
|---|
| 291 | |
|---|
| 292 | end do ! (ngridmx) |
|---|
| 293 | |
|---|
| 294 | endif ! of if (improvedicelocation) |
|---|
| 295 | #endif |
|---|
| 296 | |
|---|
| 297 | ENDIF ! (caps & water) |
|---|
| 298 | |
|---|
| 299 | c =============================================================== |
|---|
| 300 | c INITIAL ALBEDO |
|---|
| 301 | c =============================================================== |
|---|
| 302 | |
|---|
| 303 | write(*,*)"surfini: water frost thickness", |
|---|
| 304 | s frost_albedo_threshold |
|---|
| 305 | write(*,*)"surfini: water ice albedo:", albedo_h2o_ice |
|---|
| 306 | write(*,*)"surfini: water ice TI:", inert_h2o_ice |
|---|
| 307 | |
|---|
| 308 | c To start with : Initial albedo = observed dataset |
|---|
| 309 | c ------------------------------------------------- |
|---|
| 310 | DO ig=1,ngrid |
|---|
| 311 | psolaralb(ig,1)=albedodat(ig) |
|---|
| 312 | psolaralb(ig,2)=albedodat(ig) |
|---|
| 313 | END DO |
|---|
| 314 | PRINT*,'minimum albedo sans water caps', |
|---|
| 315 | s albedodat(ISMIN(ngrid,albedodat,1)) |
|---|
| 316 | PRINT*,'maximum albedo sans water caps', |
|---|
| 317 | s albedodat(ISMAX(ngrid,albedodat,1)) |
|---|
| 318 | |
|---|
| 319 | c initial albedo if permanent H2O ice is present |
|---|
| 320 | c ------------------------------------------------ |
|---|
| 321 | IF ((water) .and. (caps)) THEN |
|---|
| 322 | DO ig=1,ngrid |
|---|
| 323 | IF (watercaptag(ig)) THEN |
|---|
| 324 | psolaralb(ig,1) = albedo_h2o_ice |
|---|
| 325 | psolaralb(ig,2) = albedo_h2o_ice |
|---|
| 326 | ENDIF |
|---|
| 327 | END DO |
|---|
| 328 | PRINT*,'minimum albedo avec water caps', |
|---|
| 329 | s psolaralb(ISMIN(ngrid,psolaralb,1),1) |
|---|
| 330 | PRINT*,'maximum albedo avec water caps', |
|---|
| 331 | s psolaralb(ISMAX(ngrid,psolaralb,1),1) |
|---|
| 332 | |
|---|
| 333 | ENDIF |
|---|
| 334 | |
|---|
| 335 | c changing initial albedo if CO2 ice is present |
|---|
| 336 | c ------------------------------------------- |
|---|
| 337 | |
|---|
| 338 | DO ig=1,ngrid |
|---|
| 339 | IF (piceco2(ig) .GT. 0.) THEN |
|---|
| 340 | IF(ig.GT.ngrid/2+1) THEN |
|---|
| 341 | icap=2 |
|---|
| 342 | ELSE |
|---|
| 343 | icap=1 |
|---|
| 344 | ENDIF |
|---|
| 345 | psolaralb(ig,1) = albedice(icap) |
|---|
| 346 | psolaralb(ig,2) = albedice(icap) |
|---|
| 347 | END IF |
|---|
| 348 | END DO |
|---|
| 349 | |
|---|
| 350 | c changing initial albedo if water ice frost is present |
|---|
| 351 | c ------------------------------------------- |
|---|
| 352 | IF (water) THEN |
|---|
| 353 | do iq=1,nqmx |
|---|
| 354 | c if there is frost and surface albedo is set to albedo_h2o_ice |
|---|
| 355 | if(noms(iq).eq."h2o_ice") then |
|---|
| 356 | do ig=1,ngrid |
|---|
| 357 | if ((piceco2(ig) .eq. 0.).and. |
|---|
| 358 | & (qsurf(ig,iq).gt.frost_albedo_threshold)) then |
|---|
| 359 | psolaralb(ig,1) = albedo_h2o_ice |
|---|
| 360 | psolaralb(ig,2) = albedo_h2o_ice |
|---|
| 361 | c PRINT*,'surfini.F frost', |
|---|
| 362 | c & lati(ig)*180./pi, long(ig)*180./pi |
|---|
| 363 | endif |
|---|
| 364 | enddo |
|---|
| 365 | endif |
|---|
| 366 | end do |
|---|
| 367 | PRINT*,'minimum albedo avec givre et co2', |
|---|
| 368 | s psolaralb(ISMIN(ngrid,psolaralb,1),1) |
|---|
| 369 | PRINT*,'maximum albedo avec givre et co2', |
|---|
| 370 | s psolaralb(ISMAX(ngrid,psolaralb,1),1) |
|---|
| 371 | END IF |
|---|
| 372 | |
|---|
| 373 | |
|---|
| 374 | RETURN |
|---|
| 375 | END |
|---|
