SUBROUTINE surfini(ngrid,piceco2,qsurf,psolaralb) ! to use 'getin' USE ioipsl_getincom IMPLICIT NONE c======================================================================= c c creation des calottes pour l'etat initial c c======================================================================= c----------------------------------------------------------------------- c Declarations: c ------------- #include "dimensions.h" #include "dimphys.h" #include "surfdat.h" #include "callkeys.h" #include "tracer.h" #include "comgeomfi.h" #include "comcstfi.h" c INTEGER ngrid,ig,icap,iq,alternate REAL piceco2(ngrid),psolaralb(ngrid,2) REAL qsurf(ngrid,nqmx) !tracer on surface (kg/m2) REAL icedryness ! ice dryness ! longwatercaptag is watercaptag. Trick for some compilers LOGICAL, DIMENSION(100000) :: longwatercaptag EXTERNAL ISMIN,ISMAX INTEGER ISMIN,ISMAX ! Choose false to have a somewhat non resolution dependant water ice caps distribution, ! i.e based only on lat & lon values of each physical point. ! Choose true to get a carefully choosen distribution for GCM resolutions 32x24 or 64x48 ! For vizualisation : /u/tnalmd/bin/watercaps gcm_txt_output LOGICAL,SAVE :: improvedicelocation = .true. c c======================================================================= c water ice outliers c ------------------------------------------ IF ((water) .and. (caps)) THEN c Perennial H20 north cap defined by watercaptag=true (allows surface to be c hollowed by sublimation in vdifc). c We might not want albedodat to be modified because it is used to write c restart files. Instead, albedo is directly modified when needed (i.e. c if we have watercaptag and no co2 ice), below and in albedocaps.F90 c "Dryness coefficient" controlling the evaporation and c sublimation from the ground water ice (close to 1) c HERE, the goal is to correct for the fact c that the simulated permanent water ice polar caps c is larger than the actual cap and the atmospheric c opacity not always realistic. alternate = 0 if (ngridmx .ne. 1) then watercaptag(:) = .false. longwatercaptag(:) = .false. endif write(*,*) "Ice dryness ?" icedryness=1. ! default value call getin("icedryness",icedryness) write(*,*) " icedryness = ",icedryness #ifdef MESOSCALE do ig=1,ngridmx !write(*,*) "all qsurf to zero. dirty." do iq=1,nqmx qsurf(ig,iq)=0. !! on jette les inputs GCM !! on regle juste watercaptag !! il faudrait garder les inputs GCM !! si elles sont consequentes enddo if ( ( lati(ig)*180./pi .gt. 70. ) .and. . ( albedodat(ig) .ge. 0.26 ) ) then write(*,*)"outlier ",ig watercaptag(ig) = .true. dryness(ig) = 1. albedodat(ig) = albedo_h2o_ice !! pour output else watercaptag(ig) = .false. dryness(ig) = 1. endif enddo #else if (improvedicelocation) then if (ngridmx .eq. 738) then ! hopefully 32x24 print*,'water ice caps distribution for 32x24 resolution:' longwatercaptag(1:9) = .true. ! central cap - core longwatercaptag(26:33) = .true. ! central cap else if (ngridmx .eq. 3010) then ! hopefully 64x48 ! For latitudes: ! [ 67.5 71.25 75. 78.75 82.5 86.25] ! The water ice caps should be (according to TES temperatures): ! [ 8.63e-03 1.02e+00 5.99e+00 2.66e+01 5.65e+01] print*,'water ice caps distribution for 64x48 resolution:' longwatercaptag(1:65) = .true. ! central cap - core longwatercaptag(75:85) = .true. ! central cap longwatercaptag(93:114) = .true. ! central cap !--------------------- OUTLIERS ---------------------------- longwatercaptag(136) = .true. ! outlier, lat = 78.75 longwatercaptag(138) = .true. ! outlier, lat = 78.75 longwatercaptag(140) = .true. ! outlier, lat = 78.75 longwatercaptag(142) = .true. ! outlier, lat = 78.75 longwatercaptag(161) = .true. ! outlier, lat = 78.75 longwatercaptag(163) = .true. ! outlier, lat = 78.75 longwatercaptag(165) = .true. ! outlier, lat = 78.75 longwatercaptag(183) = .true. ! outlier, lat = 78.75 longwatercaptag(185) = .true. ! outlier, lat = 78.75 longwatercaptag(187) = .true. ! outlier, lat = 78.75 longwatercaptag(189) = .true. ! outlier, lat = 78.75 longwatercaptag(191) = .true. ! outlier, lat = 78.75 longwatercaptag(193) = .true. ! outlier, lat = 78.75 longwatercaptag(194) = .true. ! outlier, lat = 75 longwatercaptag(203) = .true. ! outlier, lat = 75 longwatercaptag(207) = .true. ! outlier, lat = 75 longwatercaptag(242) = .true. ! outlier, lat = 75 longwatercaptag(244) = .true. ! outlier, lat = 75 longwatercaptag(246) = .true. ! outlier, lat = 75 longwatercaptag(248) = .true. ! outlier, lat = 75 longwatercaptag(250) = .true. ! outlier, lat = 75 longwatercaptag(252) = .true. ! outlier, lat = 75 longwatercaptag(254) = .true. ! outlier, lat = 75 longwatercaptag(256:257)= .true. ! outlier, lat = 75 !!------------------------- OLD ---------------------------- !! longwatercaptag(83:85) = .true. !! longwatercaptag(135:142) = .true. !! longwatercaptag(181:193) = .true. !! longwatercaptag(242:257) = .true. else if (ngridmx .ne. 1) then print*,'No improved ice location for this resolution!' print*,'Please set improvedicelocation to false in surfini.' print*,'And check lat and lon conditions for ice caps in code.' call exit(1) endif ! print caps locations print*,'latitude | longitude | ig' do ig=1,ngridmx dryness (ig) = icedryness if (longwatercaptag(ig)) then watercaptag(ig) = .True. print*,'ice water cap', lati(ig)*180./pi, . long(ig)*180./pi, ig endif enddo else do ig=1,ngridmx dryness (ig) = icedryness !!c Towards olympia planitia water caps ('relatively' low latitude ones) !!c---------------- proposition par AS -------------------- !!c-------------------------------------------------------- !!c if ( ( lati(ig)*180./pi .ge. 75 ) .and. !!c . ( lati(ig)*180./pi .le. 77 ) .and. !!c . ( ( ( long(ig)*180./pi .ge. 000. ) .and. !!c . ( long(ig)*180./pi .le. 120. ) ) !!c . .or. !!c . ( ( long(ig)*180./pi .ge. -130. ) .and. !!c . ( long(ig)*180./pi .le. -115. ) ) ) ) then !!c---------------- proposition par TN -------------------- !!c---------------- HIGHLY EXPERIMENTAL ------------------- !!c-------------------------------------------------------- !! if ( ( ( lati(ig)*180./pi .ge. 73. ) .and. ! cap #1 !! . ( lati(ig)*180./pi .le. 75.1 ) .and. !! . ( long(ig)*180./pi .ge. 95. ) .and. !! . ( long(ig)*180./pi .le. 110. ) ) !! . .or. !! . ( ( lati(ig)*180./pi .ge. 77. ) .and. ! cap #2 !! . ( lati(ig)*180./pi .le. 80. ) .and. !! . ( long(ig)*180./pi .ge. 110. ) .and. !! . ( long(ig)*180./pi .le. 140. ) ) !! . .or. !! . ( ( lati(ig)*180./pi .ge. 74.9 ) .and. ! cap #3 !! . ( lati(ig)*180./pi .le. 78. ) .and. !! . ( long(ig)*180./pi .ge. 155. ) .and. !! . ( long(ig)*180./pi .le. 180. ) ) !! . .or. !! . ( ( lati(ig)*180./pi .ge. 71. ) .and. ! cap #4 (Korolev crater) !! . ( lati(ig)*180./pi .le. 72. ) .and. !! . ( long(ig)*180./pi .ge. 163. ) .and. !! . ( long(ig)*180./pi .le. 168. ) ) !! . .or. !! . ( ( lati(ig)*180./pi .ge. 74.9 ) .and. ! cap #5 !! . ( lati(ig)*180./pi .le. 78. ) .and. !! . ( long(ig)*180./pi .ge. -160.) .and. !! . ( long(ig)*180./pi .le. -120.) ) ) !! . then if ( ( ( lati(ig)*180./pi .ge. 77. ) .and. ! cap #2 . ( lati(ig)*180./pi .le. 80. ) .and. . ( long(ig)*180./pi .ge. 110. ) .and. . ( long(ig)*180./pi .le. 181. ) ) . .or. . ( ( lati(ig)*180./pi .ge. 75. ) .and. ! cap #4 (Korolev crater) . ( lati(ig)*180./pi .le. 76. ) .and. . ( long(ig)*180./pi .ge. 150. ) .and. . ( long(ig)*180./pi .le. 168. ) ) . .or. . ( ( lati(ig)*180./pi .ge. 77 ) .and. ! cap #5 . ( lati(ig)*180./pi .le. 80. ) .and. . ( long(ig)*180./pi .ge. -150.) .and. . ( long(ig)*180./pi .le. -110.) ) ) . then if ((alternate .eq. 0)) then ! 1/2 en 64x48 sinon trop large en lat if (ngridmx.ne.1) watercaptag(ig)=.true. print*,'ice water cap', lati(ig)*180./pi, . long(ig)*180./pi, ig dryness(ig) = 1. alternate = 1 else alternate = 0 endif !end if alternate = 0 endif c Opposite olympia planitia water cap c---------------- proposition par AS -------------------- c-------------------------------------------------------- c if ( ( lati(ig)*180./pi .ge. 82 ) .and. c . ( lati(ig)*180./pi .le. 84 ) .and. c . ( ( long(ig)*180./pi .gt. -030. ) .and. c . ( long(ig)*180./pi .lt. 090. ) ) ) then c---------------- proposition par TN -------------------- c-------------------------------------------------------- if ( ( ( lati(ig)*180./pi .ge. 80 ) .and. . ( lati(ig)*180./pi .le. 84 ) ) . .and. . ( ( long(ig)*180./pi .lt. -95. ) .or. !!! 32x24 . ( long(ig)*180./pi .gt. 85. ) ) ) then !!! 32x24 ! . ( ( ( long(ig)*180./pi .ge. -29. ) .and. !!! 64x48 ! . ( long(ig)*180./pi .le. 90. ) ) .or. !!! 64x48 ! . ( ( long(ig)*180./pi .ge. -77. ) .and. !!! 64x48 ! . ( long(ig)*180./pi .le. -70. ) ) ) ) then !!! 64x48 if (ngridmx.ne.1) then watercaptag(ig)=.true. print*,'ice water cap', lati(ig)*180./pi, . long(ig)*180./pi, ig endif dryness(ig) = 1. endif c Central cap c---------------- anciens reglages FF ------------------- c-------------------------------------------------------- if (abs(lati(ig)*180./pi).gt.80) then print*,'ice water cap', lati(ig)*180./pi, . long(ig)*180./pi, ig if (ngridmx.ne.1) watercaptag(ig)=.true. !dryness(ig) = 1. c Use the following cap definition for high spatial resolution (latitudinal bin <= 5 deg) c if (lati(ig)*180./pi.lt.85.and.long(ig).ge.0) then c if (ngridmx.ne.1) watercaptag(ig)=.true. c dryness(ig) = 1. c endif c if (lati(ig)*180./pi.ge.85) then c if (ngridmx.ne.1) watercaptag(ig)=.true. c dryness(ig) = 1. c endif endif ! (lati>80 deg) end do ! (ngridmx) endif ! of if (improvedicelocation) #endif ENDIF ! (caps & water) c =============================================================== c INITIAL ALBEDO c =============================================================== write(*,*)"surfini: water frost thickness", s frost_albedo_threshold write(*,*)"surfini: water ice albedo:", albedo_h2o_ice write(*,*)"surfini: water ice TI:", inert_h2o_ice c To start with : Initial albedo = observed dataset c ------------------------------------------------- DO ig=1,ngrid psolaralb(ig,1)=albedodat(ig) psolaralb(ig,2)=albedodat(ig) END DO PRINT*,'minimum albedo sans water caps', s albedodat(ISMIN(ngrid,albedodat,1)) PRINT*,'maximum albedo sans water caps', s albedodat(ISMAX(ngrid,albedodat,1)) c initial albedo if permanent H2O ice is present c ------------------------------------------------ IF ((water) .and. (caps)) THEN DO ig=1,ngrid IF (watercaptag(ig)) THEN psolaralb(ig,1) = albedo_h2o_ice psolaralb(ig,2) = albedo_h2o_ice ENDIF END DO PRINT*,'minimum albedo avec water caps', s psolaralb(ISMIN(ngrid,psolaralb,1),1) PRINT*,'maximum albedo avec water caps', s psolaralb(ISMAX(ngrid,psolaralb,1),1) ENDIF c changing initial albedo if CO2 ice is present c ------------------------------------------- DO ig=1,ngrid IF (piceco2(ig) .GT. 0.) THEN IF(ig.GT.ngrid/2+1) THEN icap=2 ELSE icap=1 ENDIF psolaralb(ig,1) = albedice(icap) psolaralb(ig,2) = albedice(icap) END IF END DO c changing initial albedo if water ice frost is present c ------------------------------------------- IF (water) THEN do iq=1,nqmx c if there is frost and surface albedo is set to albedo_h2o_ice if(noms(iq).eq."h2o_ice") then do ig=1,ngrid if ((piceco2(ig) .eq. 0.).and. & (qsurf(ig,iq).gt.frost_albedo_threshold)) then psolaralb(ig,1) = albedo_h2o_ice psolaralb(ig,2) = albedo_h2o_ice c PRINT*,'surfini.F frost', c & lati(ig)*180./pi, long(ig)*180./pi endif enddo endif end do PRINT*,'minimum albedo avec givre et co2', s psolaralb(ISMIN(ngrid,psolaralb,1),1) PRINT*,'maximum albedo avec givre et co2', s psolaralb(ISMAX(ngrid,psolaralb,1),1) END IF RETURN END