Changeset 669 for trunk/LMDZ.MARS/libf


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
Location:
trunk/LMDZ.MARS/libf/phymars
Files:
4 edited

Legend:

Unmodified
Added
Removed

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

    r668 r669  
    311311      REAL mtot(ngridmx)          ! Total mass of water vapor (kg/m2)
    312312      REAL icetot(ngridmx)        ! Total mass of water ice (kg/m2)
    313       REAL ccntot(ngridmx)        ! Total number of ccn (nbr/m2)
     313      REAL Nccntot(ngridmx)       ! Total number of ccn (nbr/m2)
     314      REAL Mccntot(ngridmx)       ! Total mass of ccn (kg/m2)
    314315      REAL rave(ngridmx)          ! Mean water ice effective radius (m)
    315316      REAL opTES(ngridmx,nlayermx)! abs optical depth at 825 cm-1
    … …  
    14451446         if (tracer) then
    14461447           if (water) then
    1447            
    1448              if (scavenging) then
    1449                ccntot(:)= 0
    1450                do ig=1,ngrid
    1451                  do l=1,nlayermx
    1452                    ccntot(ig) = ccntot(ig) +
    1453      &              zq(ig,l,igcm_ccn_number)*tauscaling(ig)
    1454      &              *(pplev(ig,l) - pplev(ig,l+1)) / g
    1455                  enddo
    1456                enddo
    1457              endif
    14581448
    14591449             mtot(:)=0
    … …  
    14691459     &                        zq(ig,l,igcm_h2o_ice) *
    14701460     &                        (pplev(ig,l) - pplev(ig,l+1)) / g
    1471                  rave(ig) = rave(ig) +
    1472      &                      zq(ig,l,igcm_h2o_ice) *
    1473      &                      (pplev(ig,l) - pplev(ig,l+1)) / g *
    1474      &                      rice(ig,l) * (1.+nuice_ref)
     1461cccc Column integrated effective ice radius
     1462cccc is weighted by total ice mass (LESS GOOD than total ice surface area)
     1463c                 rave(ig) = rave(ig) +
     1464c     &                      zq(ig,l,igcm_h2o_ice) *
     1465c     &                      (pplev(ig,l) - pplev(ig,l+1)) / g *
     1466c     &                      rice(ig,l) * (1.+nuice_ref)
    14751467c                Computing abs optical depth at 825 cm-1 in each
    14761468c                  layer to simulate NEW TES retrieval
    … …  
    14841476                 tauTES(ig)=tauTES(ig)+ opTES(ig,l)
    14851477               enddo
    1486                rave(ig)=rave(ig)/max(icetot(ig),1.e-30)
     1478c              rave(ig)=rave(ig)/max(icetot(ig),1.e-30)       ! mass weight
     1479c               if (icetot(ig)*1e3.lt.0.01) rave(ig)=0.
     1480             enddo
     1481             call watersat(ngridmx*nlayermx,zt,pplay,zqsat)
     1482             satu(:,:) = zq(:,:,igcm_h2o_vap)/zqsat(:,:)
     1483
     1484             if (scavenging) then
     1485               Nccntot(:)= 0
     1486               Mccntot(:)= 0
     1487               rave(:)=0
     1488               do ig=1,ngrid
     1489                 do l=1,nlayermx
     1490                    Nccntot(ig) = Nccntot(ig) +
     1491     &              zq(ig,l,igcm_ccn_number)*tauscaling(ig)
     1492     &              *(pplev(ig,l) - pplev(ig,l+1)) / g
     1493                    Mccntot(ig) = Mccntot(ig) +
     1494     &              zq(ig,l,igcm_ccn_mass)*tauscaling(ig)
     1495     &              *(pplev(ig,l) - pplev(ig,l+1)) / g
     1496cccc Column integrated effective ice radius
     1497cccc is weighted by total ice surface area (BETTER than total ice mass)
     1498                    rave(ig) = rave(ig) +
     1499     &                      tauscaling(ig) *
     1500     &                      zq(ig,l,igcm_ccn_number) *
     1501     &                      (pplev(ig,l) - pplev(ig,l+1)) / g *
     1502     &                      rice(ig,l) * rice(ig,l)*  (1.+nuice_ref)
     1503                 enddo
     1504               rave(ig)=(icetot(ig)/rho_ice+Mccntot(ig)/rho_dust)*0.75
     1505     &                  /max(pi*rave(ig),1.e-30) ! surface weight
    14871506               if (icetot(ig)*1e3.lt.0.01) rave(ig)=0.
    1488              enddo
     1507               enddo
     1508             endif ! of if (scavenging)
    14891509
    14901510           endif ! of if (water)
    … …  
    15481568     &                    "H2O ice volume mixing ratio","mol/mol",
    15491569     &                    3,vmr)
     1570               vmr=zqsat(1:ngridmx,1:nlayermx)
     1571     &                  *mugaz/mmol(igcm_h2o_vap)
     1572               call wstats(ngrid,"vmr_h2osat",
     1573     &                    "saturation volume mixing ratio","mol/mol",
     1574     &                    3,vmr)
    15501575               call wstats(ngrid,"h2o_ice_s",
    15511576     &                    "surface h2o_ice","kg/m2",
    … …  
    15631588     &                    "Mean reff","m",
    15641589     &                    2,rave)
    1565               call wstats(ngrid,"ccntot",
     1590              call wstats(ngrid,"Nccntot",
    15661591     &                    "condensation nuclei","Nbr/m2",
    1567      &                    2,ccntot)
     1592     &                    2,Nccntot)
     1593              call wstats(ngrid,"Mccntot",
     1594     &                    "condensation nuclei mass","kg/m2",
     1595     &                    2,Mccntot)
    15681596              call wstats(ngrid,"rice",
    15691597     &                    "Ice particle size","m",
    … …  
    17951823c            call WRITEDIAGFI(ngridmx,'vmr_h2ovap','h2o vap vmr',
    17961824c     &                       'mol/mol',3,vmr)
    1797 c             CALL WRITEDIAGFI(ngridmx,'reffice',
    1798 c     &                       'Mean reff',
    1799 c     &                       'm',2,rave)
    1800 c             CALL WRITEDIAGFI(ngrid,"ccntot",
    1801 c     &                    "condensation nuclei","Nbr/m2",
    1802 c     &                    2,ccntot)
     1825             CALL WRITEDIAGFI(ngridmx,'reffice',
     1826     &                       'Mean reff',
     1827     &                       'm',2,rave)
     1828             CALL WRITEDIAGFI(ngrid,"Nccntot",
     1829     &                    "condensation nuclei","Nbr/m2",
     1830     &                    2,Nccntot)
     1831c             CALL WRITEDIAGFI(ngrid,"Mccntot",
     1832c     &                    "mass condensation nuclei","kg/m2",
     1833c     &                    2,Mccntot)
    18031834c            call WRITEDIAGFI(ngridmx,'rice','Ice particle size',
    18041835c     &                       'm',3,rice)
    … …  
    20542085             icetot = icetot +  zq(1,l,igcm_h2o_ice)
    20552086     &                 * (pplev(1,l) - pplev(1,l+1)) / g
    2056              rave = rave + zq(1,l,igcm_h2o_ice)
    2057      &                 * (pplev(1,l) - pplev(1,l+1)) / g
    2058      &                 *  rice(1,l) * (1.+nuice_ref)
     2087cccc Column integrated effective ice radius
     2088cccc is weighted by total ice surface area (BETTER)
     2089              rave = rave + tauscaling(ig) *
     2090     &               zq(1,l,igcm_ccn_number) *
     2091     &               (pplev(1,l) - pplev(1,l+1)) / g *
     2092     &               rice(1,l) * rice(1,l)*  (1.+nuice_ref)
     2093cccc Column integrated effective ice radius
     2094cccc is weighted by total ice mass         (LESS GOOD)
     2095c             rave = rave + zq(1,l,igcm_h2o_ice)
     2096c     &                 * (pplev(1,l) - pplev(1,l+1)) / g
     2097c     &                 *  rice(1,l) * (1.+nuice_ref)
    20592098           end do
    2060              rave=rave/max(icetot,1.e-30)
     2099           rave=icetot*0.75/max(rave*pi*rho_ice,1.e-30) ! surface weight
     2100c           rave=rave/max(icetot,1.e-30)    ! mass weight
    20612101            h2otot = h2otot+mtot+icetot
    20622102           
    20632103           
    20642104             if (scavenging) then
    2065                ccntot= 0
     2105               Nccntot= 0
    20662106              call watersat(ngridmx*nlayermx,zt,pplay,zqsat)
    20672107               do l=1,nlayermx
    2068                    ccntot = ccntot +
     2108                   Nccntot = Nccntot +
    20692109     &              zq(1,l,igcm_ccn_number)*tauscaling(1)
    20702110     &              *(pplev(1,l) - pplev(1,l+1)) / g
    … …  
    20752115               enddo
    20762116
    2077                CALL WRITEDIAGFI(ngridmx,'ccntot',
    2078      &                         'ccntot',
    2079      &                         'nbr/m2',0,ccntot)
     2117               CALL WRITEDIAGFI(ngridmx,'Nccntot',
     2118     &                         'Nccntot',
     2119     &                         'nbr/m2',0,Nccntot)
    20802120             endif
    20812121             

  • 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

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

    r552 r669  
    6767      INTEGER Lmodif
    6868      INTEGER lmax
    69       REAL p_rad,p_omeg,p_g,p_mugaz,p_daysec,time
     69      REAL p_rad,p_omeg,p_g,p_mugaz,p_daysec,time,peri_ls
    7070
    7171c Variables
    … …  
    289289      write(*,*) '(18)        obliquit : planet obliquity (deg)'
    290290      write(*,*) '(17)      peri_day : perihelion date (sol since Ls=0)'
     291      write(*,*) '(  )      peri_ls : perihelion date (Ls since Ls=0)'
    291292      write(*,*) '(15)      periheli : min. sun-mars dist (Mkm)'
    292293      write(*,*) '(16)      aphelie  : max. sun-mars dist (Mkm)'
    … …  
    429430        else if (trim(modif) .eq. 'peri_day') then
    430431          write(*,*) 'current value:',peri_day
    431           write(*,*) 'peri_day should be 485 on current Mars'
     432          write(*,*) 'peri_day should be 485 sols on current Mars'
    432433          write(*,*) 'enter new value:'
    433434 116      read(*,*,iostat=ierr) peri_day
    … …  
    435436          write(*,*)
    436437          write(*,*) ' peri_day (new value):',peri_day
     438         
     439        else if (trim(modif) .eq. 'peri_ls') then
     440          write(*,*) 'peri_ls value is not stored in start files,'
     441          write(*,*) 'but it should be 251 degrees on current Mars'
     442          write(*,*) '(peri_day should be 485 sols on current Mars)'
     443          write(*,*) 'enter new value:'
     444 1160     read(*,*,iostat=ierr) peri_ls
     445          if(ierr.ne.0) goto 1160
     446          write(*,*)
     447          write(*,*) ' peri_ls asked :',peri_ls
     448          call lsp2solp(peri_ls,peri_day)
     449          write(*,*) ' peri_day (new value):',peri_day
     450
    437451
    438452        else if (trim(modif) .eq. 'periheli') then
    … …  
    529543c-----------------------------------------------------------------------
    530544      end
     545
     546
     547
     548     
     549!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     550! gives sol at perihelion for ls at perihelion (for precession cycles)
     551      subroutine lsp2solp(lsp,solp)
     552
     553      implicit none
     554!  Arguments:
     555      real lsp     ! Input: ls at perihelion
     556      real solp    ! Output: sol at perihelion
     557     
     558!  Local:
     559      double precision zx0 ! eccentric anomaly at Ls=0
     560      double precision year_day
     561      double precision e_elips
     562      double precision pi,degrad
     563     
     564      parameter (year_day=668.6d0) ! number of sols in a martian year
     565      parameter (e_elips=0.0934d0)  ! eccentricity of orbit
     566      parameter (pi=3.14159265358979d0)
     567      parameter (degrad=57.2957795130823d0)
     568     
     569      zx0 = -2.0*datan(dtan(0.5*lsp/degrad)
     570     .          *dsqrt((1.-e_elips)/(1.+e_elips)))
     571      if (zx0 .le. 0.) zx0 = zx0 + 2.*pi
     572     
     573      solp  = year_day*(1.-(zx0-e_elips*dsin(zx0))/(2.*pi))
     574
     575
     576      end subroutine lsp2solp
     577
     578
     579

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

    r660 r669  
    102102      SAVE firstcall
    103103
    104       REAL lw
    105104
    106105c     variable added for CO2 condensation:
    … …  
    111110      SAVE acond,bcond
    112111      DATA latcond,tcond1mb/5.9e5,136.27/
     112     
     113c     For latent heat release from ground ice sublimation   
     114      REAL tsrf_lw(ngridmx)
     115      REAL T1,T2
     116      SAVE T1,T2
     117      DATA T1,T2/-877.5,807.5/ ! zeros of latent heat equation for ice
    113118
    114119c    Tracers :
    … …  
    651656
    652657      pdtsrf(:)=(ztsrf2(:)-ptsrf(:))/ptimestep
    653 
     658     
    654659      DO ig=1,ngrid  ! computing sensible heat flux (atm => surface)
    655660         sensibFlux(ig)=cpp*zb(ig,1)/ptimestep*(zhs(ig,1)-ztsrf2(ig))
    … …  
    804809              endif ! if (.not.watercaptag(ig))
    805810c             Starting upward calculations for water :
    806                zq(ig,1,igcm_h2o_vap)=zq1temp(ig)
    807 c             Take into account H2o latent heat in surface energy budget         
    808                lw = (2834.3-0.28*(ptsrf(ig)-To)
    809      &              -0.004*(ptsrf(ig)-To)**2)*1.e+3
    810                pdtsrf(ig) = pdtsrf(ig)
    811      &                    + pdqsdif(ig,igcm_h2o_ice)*lw/pcapcal(ig)
     811              zq(ig,1,igcm_h2o_vap)=zq1temp(ig)
     812             
     813c             Take into account H2O latent heat in surface energy budget
     814c             We solve dT/dt = (2834.3-0.28*(T-To)-0.004*(T-To)^2)*1e3*iceflux/cpp
     815              tsrf_lw(ig) = ptsrf(ig) + pdtsrf(ig) *ptimestep
     816             
     817              tsrf_lw(ig) = (T1+T2)*(T1+T2)
     818     &       -  4*(T2*T1 - (tsrf_lw(ig)-T1)*(tsrf_lw(ig)-T2)*
     819     &            exp( -0.25*abs(T1-T2)*pdqsdif(ig,igcm_h2o_ice)
     820     &            *ptimestep/pcapcal(ig)) )       
     821              tsrf_lw(ig) = (T1+T2)/2 +  sqrt(tsrf_lw(ig))/2 ! surface temperature at t+1
     822
     823              pdtsrf(ig) = (tsrf_lw(ig)-ptsrf(ig))/ptimestep
     824
     825               if(pqsurf(ig,igcm_h2o_ice)
     826     &           +pdqsdif(ig,igcm_h2o_ice)*ptimestep
     827     &           .gt.frost_albedo_threshold) ! if there is still ice, T cannot exceed To
     828     &           pdtsrf(ig) = min(pdtsrf(ig),(To-ptsrf(ig))/ptimestep) ! ice melt case
     829     
    812830            ENDDO ! of DO ig=1,ngrid
    813831          ELSE
    … …  
    825843        enddo ! of do iq=1,nq
    826844      end if ! of if(tracer)
     845     
    827846
    828847c-----------------------------------------------------------------------
Note: See TracChangeset for help on using the changeset viewer.