Ignore:
Timestamp:
May 24, 2012, 5:02:41 PM (13 years ago)
Author:
tnavarro
Message:

Ice effective radius for output is weighted by ice total surface instead of ice total mass in physiq.F

  • New distribution for permanent ice reservoirs in surfini.F

icelocationmode = 1 allows for realistic ice distribution, no matter what resolution.
icelocationmode = 2 is predefined 32x24 or 64x48 and should be USED BY DEFAULT. It currently overestimates ice.
icelocationmode = 3 is good for quick logical definitions (paleoclimates,stability studies,etc...)

  • Possibility to change perihelion date in solar longitude as well as in sol in tabfi.F for newstarts
  • Improved latent heat release when ground ice sublimates: now a implicit scheme for ground temperature
File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/LMDZ.MARS/libf/phymars/surfini.F

    r643 r669  
    1818#include "comgeomfi.h"
    1919#include "comcstfi.h"
    20 c
     20
     21#include "datafile.h"
     22#include "netcdf.inc"
     23
    2124      INTEGER ngrid,ig,icap,iq,alternate
    2225      REAL  piceco2(ngrid),psolaralb(ngrid,2)
    … …  
    3033      INTEGER ISMIN,ISMAX
    3134     
    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.
     35! There are 3 different modes for ice distribution:
     36! icelocationmode = 1 ---> based on data from surface.nc
     37! icelocationmode = 2 ---> directly predefined for GCM resolutions 32x24 or 64x48
     38! icelocationmode = 3 ---> based on logical relations for latitude and longitude
     39! For visualisation : > /u/tnalmd/bin/watercaps gcm_txt_output_file
     40      INTEGER,SAVE :: icelocationmode = 2
     41       
     42       
     43      !in case icelocationmode == 1
     44      INTEGER i,j
     45      INTEGER     imd,jmd
     46      PARAMETER   (imd=360,jmd=180)
     47      REAL        zdata(imd*jmd)
     48      REAL        zelat,zelon
     49! on a data(i,j) = zdata((j-1)*jmd + i) avec i longitude de 1 a imd (E vers O) et j lat de 1 a jmd (N vers S)
     50
     51      INTEGER nb_ice(ngrid,2)              ! number of counts | detected ice for GCM grid
     52      INTEGER latice(jjm,2),lonice (iim,2) ! number of counts | detected ice along lat & lon axis
     53
     54      REAL step,count,ratiolat
     55
     56      INTEGER   ierr,nid,nvarid
     57     
     58      REAL,SAVE :: min_icevalue = 500.
     59      PARAMETER string = 'thermal'
     60     
     61      character (len=100) :: zedatafile
    3762c
    3863c=======================================================================
    … …  
    6893         call getin("icedryness",icedryness)
    6994         write(*,*) " icedryness = ",icedryness
     95         dryness (:) = icedryness
    7096         
    7197       
    … …  
    93119         
    94120      enddo
    95 
    96121#else
    97122
    98       if (improvedicelocation) then
    99      
    100         if (ngridmx .eq. 738) then ! hopefully 32x24
     123
     124
     125      IF (ngridmx .eq. 1) THEN ! special case for 1d --> do nothing
     126     
     127         print*, 'ngridmx = 1, do no put ice caps in surfini.F'
     128
     129      ELSE IF (icelocationmode .eq. 1) THEN
     130     
     131         print*,'Surfini: ice caps defined from surface.nc'
     132           
     133! This method detects ice as gridded value above min_icevalue in the field "string" from surface.nc
     134! Typically, it is for thermal inertia above 500 tiu.
     135! Two conditions are verified:
     136! 1. GCM ice caps are defined such as area is conserved for a given latitude
     137! (the approximation is that all points within the GCM latitude resolution have the same area).
     138! 2. caps are placed to fill the GCM points with the most detected ice first.
     139     
     140
    101141           
    102           print*,'water ice caps distribution for 32x24 resolution:'
     142        zedatafile = trim(datafile)
     143 
     144       
     145        ierr =nf_open (trim(zedatafile)//'/surface.nc',
     146     &  NF_NOWRITE,nid)
     147     
     148      IF (ierr.NE.nf_noerr) THEN
     149       write(*,*)'Error : cannot open file surface.nc '
     150       write(*,*)'(in phymars/surfini.F)'
     151       write(*,*)'It should be in :',trim(zedatafile),'/'
     152       write(*,*)'1) You can set this path in the callphys.def file:'
     153       write(*,*)'   datadir=/path/to/the/datafiles'
     154       write(*,*)'2) If necessary, surface.nc (and other datafiles)'
     155       write(*,*)'   can be obtained online on:'
     156       write(*,*)' http://www.lmd.jussieu.fr/~forget/datagcm/datafile'
     157       CALL ABORT
     158      ENDIF
     159     
     160     
     161      ierr = NF_INQ_VARID (nid, string, nvarid)
     162      if (ierr.ne.nf_noerr) then
     163        write(*,*) 'datareadnc error, cannot find ',trim(string)
     164        write(*,*) ' in file ',trim(zedatafile),'/surface.nc'
     165        stop
     166      endif
     167#ifdef NC_DOUBLE
     168      ierr = NF_GET_VAR_DOUBLE(nid, nvarid, zdata)
     169#else
     170      ierr = NF_GET_VAR_REAL(nid, nvarid, zdata)
     171#endif
     172      if (ierr.ne.nf_noerr) then
     173        write(*,*) 'datareadnc: error failed loading ',trim(string)
     174        stop
     175      endif
     176 
     177                     
     178      ierr=nf_close(nid)
     179 
     180
     181      nb_ice(:,1) = 1 ! default: there is no ice
     182      latice(:,1) = 1
     183      lonice(:,1) = 1
     184      nb_ice(:,2) = 0
     185      latice(:,2) = 0
     186      lonice(:,2) = 0
     187      !print*,'jjm,iim',jjm,iim ! jjm =  nb lati , iim = nb longi
     188
     189      ! loop over the GCM grid - except for poles (ig=1 and ngridmx)
     190      do ig=2,ngridmx-1
     191     
     192        ! loop over the surface file grid
     193        do i=1,imd*jmd
     194         zelat = 90. - 90./real(jmd) - (i-1)/imd      ! surface file latitude
     195         zelon = mod(i-1,imd) - (180.-180./real(imd)) ! surface file longitude
     196                 
     197          if ((abs(lati(ig)*180./pi-zelat) .le. 90./real(jjm)) .and.
     198     &        (abs(long(ig)*180./pi-zelon) .le. 180./real(iim))) then
     199              ! count all points in that GCM grid point
     200              nb_ice(ig,1) = nb_ice(ig,1) + 1
     201              if (zdata(i) > min_icevalue)
     202                 ! count all detected points in that GCM grid point
     203     &           nb_ice(ig,2) = nb_ice(ig,2) + 1
     204          endif
     205         
     206        enddo ! of do i=1,imd*jmd
     207       
     208        ! projection of nb_ice on GCM lat and lon axes
     209        latice(1+(ig-2)/iim,:) =
     210     &     latice(1+(ig-2)/iim,:) + nb_ice(ig,:)
     211        lonice(1+mod(ig-2,iim),:) =
     212     &     lonice(1+mod(ig-2,iim),:) + nb_ice(ig,:) ! lonice is USELESS ...
     213
     214      enddo ! of do ig=2,ngridmx-1
     215     
     216
     217     
     218      ! special case for poles
     219      nb_ice(1,2)   = 1  ! ice prescribed on north pole
     220      latice(1,:)   = nb_ice(1,:)
     221      lonice(1,:)   = nb_ice(1,:)
     222      latice(jjm,:) = nb_ice(ngridmx,:)
     223      lonice(iim,:) = nb_ice(ngridmx,:)
     224     
     225     
     226  !    print*, 'latice TOT', latice(:,1)
     227  !    print*, 'latice FOUND', latice(:,2)
     228  !    print*, 'lonice TOT', lonice(:,1)
     229  !    print*, 'lonice FOUND', lonice(:,2)
     230     
     231  !    print*, 'lat ratio', int(real(latice(:,2))/real(latice(:,1))*iim)
     232  !    print*, 'lon ratio', int(real(lonice(:,2))/real(lonice(:,1))*jjm)
     233     
     234  !    print*,''
     235  !    print*,'sum lat', sum(latice(:,1)), sum(lonice(:,1))
     236  !    print*,'sum lon', sum(latice(:,2)), sum(lonice(:,2))
     237     
     238   
     239      ! loop over GCM latitudes. CONSIDER ONLY NORTHERN HEMISPHERE
     240      do i=1,jjm/2
     241      step  = 1. ! threshold to add ice cap
     242      count = 0. ! number of ice GCM caps at this latitude
     243      ! ratiolat is the ratio of area covered by ice within this GCM latitude range
     244      ratiolat  = real(latice(i,2))/real(latice(i,1))
     245      !print*,'i',i,(i-1)*iim+2,i*iim+1
     246     
     247        ! put ice caps while there is not enough ice,
     248        ! as long as the threshold is above 20%
     249        do while ( (count .le. ratiolat*iim ) .and. (step .ge. 0.2))
     250          count = 0.
     251          ! loop over GCM longitudes
     252          do j=1,iim
     253            ! if the detected ice ratio in the GCM grid point
     254            ! is more than 'step', then add ice
     255            if (real(nb_ice((i-1)*iim+1+j,2))
     256     &        / real(nb_ice((i-1)*iim+1+j,1)) .ge. step) then
     257                  watercaptag((i-1)*iim+1+j) = .true.
     258                  count = count + 1
     259            endif
     260          enddo ! of do j=1,iim
     261          !print*, 'step',step,count,ratiolat*iim
     262          step = step - 0.01
     263        enddo ! of do while
     264      !print*, 'step',step,count,ratiolat*iim
     265
     266      enddo ! of do i=1,jjm/2
     267           
     268
     269      ELSE IF (icelocationmode .eq. 2) THEN
     270     
     271        print*,'Surfini: predefined ice caps'
     272     
     273        if ((iim .eq. 32) .and. (jjm .eq. 24)) then ! 32x24
     274           
     275          print*,'water ice caps distribution for 32x24 resolution'
    103276          longwatercaptag(1:9)    = .true. ! central cap - core
    104277          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:'
     278          longwatercaptag(1:33)  = .true. ! central cap
     279          longwatercaptag(56)  = .true. ! central cap
     280          longwatercaptag(58)  = .true. ! central cap
     281          longwatercaptag(60)  = .true. ! central cap
     282          longwatercaptag(62)  = .true. ! central cap
     283          longwatercaptag(64)  = .true. ! central cap
     284!---------------------   OUTLIERS  ----------------------------
     285
     286        else if ((iim .eq. 64) .and. (jjm .eq. 48)) then ! 64x48
     287
     288          print*,'water ice caps distribution for 64x48 resolution'
    114289          longwatercaptag(1:65)   = .true. ! central cap - core
    115290          longwatercaptag(75:85)  = .true. ! central cap
    116291          longwatercaptag(93:114) = .true. ! central cap
    117 !---------------------   OUTLIERS  ----------------------------         
     292!---------------------   OUTLIERS  ----------------------------
     293          if (.true.) then
    118294          longwatercaptag(136)    = .true. ! outlier, lat = 78.75
    119295          longwatercaptag(138)    = .true. ! outlier, lat = 78.75
    … …  
    140316          longwatercaptag(254)    = .true. ! outlier, lat = 75
    141317          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.
     318          endif
     319!--------------------------------------------------------------       
     320
    147321
    148322           
    149323        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)
     324       
     325      print*,'No predefined ice location for this resolution :',iim,jjm
     326      print*,'Please change icelocationmode in surfini.F'
     327      print*,'Or add some new definitions ...'
     328      call abort
    154329         
    155330        endif
    156        
    157         ! print caps locations
    158         print*,'latitude | longitude | ig'
     331
    159332        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
     333          if (longwatercaptag(ig)) watercaptag(ig) = .true.
    167334        enddo
    168335
    169336
    170       else
     337      ELSE IF (icelocationmode .eq. 3) THEN
     338     
     339        print*,'Surfini: ice caps defined by lat and lon values'
    171340
    172341         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
     342         
     343c-------- Towards olympia planitia water caps -----------
     344c--------------------------------------------------------
     345
    214346       if ( ( ( lati(ig)*180./pi .ge. 77.  ) .and. ! cap #2
    215347     .           ( lati(ig)*180./pi .le. 80.  ) .and.
    … …  
    228360     .           ( long(ig)*180./pi .le. -110.) ) )
    229361     .         then
    230      
    231362             
    232363               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,
     364              !    watercaptag(ig)=.true.
     365                  alternate = 1
     366               else
     367                  alternate = 0
     368               endif !end if alternate = 0
     369               
     370       endif
     371
     372c----------- Opposite olympia planitia water cap --------
     373c--------------------------------------------------------
     374
     375        if ( ( ( lati(ig)*180./pi     .ge.  80 ) .and.
     376     .         ( lati(ig)*180./pi     .le.  84 ) )
     377     .         .and.
     378     .       ( ( long(ig)*180./pi .lt. -95. ) .or.       !!! 32x24
     379     .         ( long(ig)*180./pi .gt.  85. ) ) ) then   !!! 32x24
     380!     .     ( ( ( long(ig)*180./pi .ge. -29. ) .and.       !!! 64x48
     381!     .         ( long(ig)*180./pi .le.  90. ) ) .or.      !!! 64x48
     382!     .       ( ( long(ig)*180./pi .ge. -77. ) .and.       !!! 64x48
     383!     .         ( long(ig)*180./pi .le. -70. ) ) ) ) then  !!! 64x48
     384        !   watercaptag(ig)=.true.
     385        endif
     386
     387
     388c -------------------- Central cap ----------------------
     389c--------------------------------------------------------
     390
     391      if (abs(lati(ig)*180./pi).gt.80)
     392     .       watercaptag(ig)=.true.
     393           
     394c--------------------------------------------------------
     395c--------------------------------------------------------
     396      end do ! of (ngridmx)
     397
     398
     399       ELSE
     400     
     401         print*, 'In surfini.F, icelocationmode is ', icelocationmode
     402         print*, 'It should be 1, 2 or 3.'
     403         call abort
     404
     405       ENDIF ! of if (icelocation)
     406       
     407       
     408        ! print caps locations - useful for plots too
     409        print*,'latitude | longitude | ig'
     410        do ig=1,ngridmx
     411          dryness (ig) = icedryness
     412
     413          if (watercaptag(ig)) then
     414             print*,'ice water cap', lati(ig)*180./pi,
    235415     .              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)
     416          endif
     417        enddo
     418       
    295419#endif     
    296420
    297421       ENDIF ! (caps & water)
     422       
    298423
    299424c ===============================================================
    … …  
    371496       END IF
    372497         
    373 
    374498      RETURN
    375499      END
Note: See TracChangeset for help on using the changeset viewer.