Changeset 283 for trunk/LMDZ.MARS/libf


Ignore:
Timestamp:
Sep 7, 2011, 12:57:29 PM (13 years ago)
Author:
aslmd
Message:

LMDZ.MARS:

AS: J'ai fait ce commit et fait a la main un merge avec les modifications entre temps

24/08/11 == TN

Attempts to tune the water cycle by adding outliers

+ A few structural changes !!

  • watercap.h is now obsolete and removed -- all is in surfdat.h
  • watercaptag initialized in surfini.F (up to 7 areas defined) instead of initracer.F
    • settings proposed by AS commented
    • experiments by TN decommented. use with caution.
  • water ice albedo and thermal inertia in callphys.def and inifis.F
  • water ice albedo in surfini.F
  • water ice albedo computation in albedocaps.F90
  • alb_surfice is now obsolete in physiq.F, albedo_h2o_ice is used instead
  • frost_albedo_threshold defined in surfdat.h
  • water ice thermal inertia in soil.F

TODO: * calibrate thermal inertia and ice albedo

  • have a look at subgrid-scale ice with dryness ?
Location:
trunk/LMDZ.MARS/libf/phymars
Files:
1 deleted
12 edited

Legend:

Unmodified
Added
Removed

  • trunk/LMDZ.MARS/libf/phymars/albedocaps.F90

    r38 r283  
    1212#include"dimphys.h"
    1313#include"surfdat.h"
     14#include"callkeys.h"
    1415
    1516! arguments:
    … …  
    7071      psolaralb(ig,2)=albedice(icap)
    7172    endif
     73  else if (watercaptag(ig) .and. water) then
     74    ! there is a water ice cap: set the surface albedo to the water ice one
     75    ! to do : emissivity
     76    !write(*,*) "watercaptag in albedocaps:",ig
     77    emisref(ig) = 1
     78    psolaralb(ig,1)=albedo_h2o_ice
     79    psolaralb(ig,2)=albedo_h2o_ice
    7280  else
    7381    ! set emissivity of surface to be bare ground emissivity
    … …  
    7886  endif ! of if (piceco2(ig).gt.0)
    7987enddo ! of ig=1,ngrid
    80 
    8188end subroutine albedocaps
    8289

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

    r234 r283  
    412412         call getin("caps",caps)
    413413         write(*,*) " caps = ",caps
     414
     415
     416         write(*,*) "water ice albedo ?"
     417         albedo_h2o_ice=0.45
     418         call getin("albedo_h2o_ice",albedo_h2o_ice)
     419         write(*,*) " albedo_h2o_ice = ",albedo_h2o_ice
     420         
     421         write(*,*) "water ice thermal inertia ?"
     422         inert_h2o_ice=2400 ! (J.m^-2.K^-1.s^-1/2)
     423         call getin("inert_h2o_ice",inert_h2o_ice)
     424         write(*,*) " inert_h2o_ice = ",inert_h2o_ice
     425         
     426         write(*,*) "frost thickness threshold for albedo ?"
     427         frost_albedo_threshold=0.005 ! 5.4 µm (i.e 0.005 kg/m²)
     428         call getin("frost_albedo_threshold",
     429     &    frost_albedo_threshold)
     430         write(*,*) " frost_albedo_threshold = ",
     431     &            frost_albedo_threshold
     432
     433!!!!!!!!!!!!!!!! TEMP !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     434!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!     
     435         write(*,*) "temp tag for water caps ?"
     436         temptag=.false.
     437         call getin("temptag",temptag)
     438         write(*,*) " temptag = ",temptag
     439!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     440!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!             
    414441
    415442         write(*,*) "photochemistry: include chemical species"

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

    r171 r283  
    4343#include "advtrac.h"
    4444#include "comgeomfi.h"
    45 #include "watercap.h"
    4645#include "chimiedata.h"
    4746
    … …  
    5049      real qsurf(ngridmx,nqmx)       ! tracer on surface (e.g.  kg.m-2)
    5150      real co2ice(ngridmx)           ! co2 ice mass on surface (e.g.  kg.m-2)
    52       integer iq,ig,count, yeyey
     51      integer iq,ig,count
    5352      real r0_lift , reff_lift, nueff_lift
    5453c     Ratio of small over large dust particles (used when both
    … …  
    436435       
    437436
    438 
    439437      enddo ! of do iq=1,nqmx
    440438!      count=count+nbqchem
    … …  
    587585         alpha_lift(igcm_h2o_vap) =0.
    588586         alpha_devil(igcm_h2o_vap)=0.
    589 
    590 c       "Dryness coefficient" controlling the evaporation and
    591 c        sublimation from the ground water ice (close to 1)
    592 c        HERE, the goal is to correct for the fact
    593 c        that the simulated permanent water ice polar caps
    594 c        is larger than the actual cap and the atmospheric
    595 c        opacity not always realistic.
    596 
    597          do ig=1,ngridmx
    598            if (ngridmx.ne.1) watercaptag(ig)=.false.
    599            dryness(ig) = 1.
    600          enddo
    601 
    602          IF (caps) THEN
    603 c Perennial H20 north cap defined by watercaptag=true (allows surface to be
    604 c hollowed by sublimation in vdifc).
    605          yeyey = 0
    606          do ig=1,ngridmx
    607 !          !!! TESTS TESTS outliers
    608 !          !!! TESTS TESTS outliers
    609 !          if ( ( lati(ig)*180./pi      .ge.  75 ) .and.
    610 !     .         ( lati(ig)*180./pi      .le.  77 ) .and.
    611 !     .         ( ( ( long(ig)*180./pi .ge. 000. ) .and.
    612 !     .              ( long(ig)*180./pi .le. 120. ) )
    613 !     .             .or.
    614 !     .             ( ( long(ig)*180./pi .ge. -130. ) .and.
    615 !     .             ( long(ig)*180./pi .le. -115. ) ) ) ) then
    616 !             if (yeyey .eq. 0) then  !!! 1/2 en 64x48 sinon trop large en lat
    617 !              write(*,*) "outliers ", lati(ig)*180./pi, long(ig)*180./pi
    618 !              if (ngridmx.ne.1) watercaptag(ig)=.true.
    619 !              dryness(ig) = 1.
    620 !              albedodat(ig) = 0.45 !! comme alb_surfice
    621 !              yeyey = 1
    622 !             else
    623 !              yeyey = 0
    624 !             endif
    625 !          endif
    626 !          !!! TESTS TESTS outliers
    627 !          !!! TESTS TESTS outliers
    628 !
    629 !          !!! TESTS TESTS addcap
    630 !          !!! TESTS TESTS addcap     
    631 !          if ( ( lati(ig)*180./pi      .ge.  82 ) .and.
    632 !     .         ( lati(ig)*180./pi      .le.  84 ) .and.
    633 !     .         ( ( long(ig)*180./pi .gt. -030. ) .and.
    634 !     .              ( long(ig)*180./pi .lt. 090. ) ) ) then
    635 !              write(*,*) "capadd ", lati(ig)*180./pi, long(ig)*180./pi
    636 !              if (ngridmx.ne.1) watercaptag(ig)=.true.
    637 !              albedodat(ig) = 0.45 !! comme alb_surfice
    638 !              dryness(ig) = 1.
    639 !          endif
    640 !          !!! TESTS TESTS addcap
    641 !          !!! TESTS TESTS addcap
    642    
    643            if (lati(ig)*180./pi.gt.84) then
    644              if (ngridmx.ne.1) watercaptag(ig)=.true.
    645              dryness(ig) = 1.
    646 c Use the following cap definition for high spatial resolution (latitudinal bin <= 5 deg)
    647 c             if (lati(ig)*180./pi.lt.85.and.long(ig).ge.0) then
    648 c               if (ngridmx.ne.1) watercaptag(ig)=.true.
    649 c               dryness(ig) = 1.
    650 c             endif
    651 c             if (lati(ig)*180./pi.ge.85) then
    652 c               if (ngridmx.ne.1) watercaptag(ig)=.true.
    653 c               dryness(ig) = 1.
    654 c             endif
    655            endif  ! (lati>80 deg)
    656          end do ! (ngridmx)
    657         ENDIF ! (caps)
    658 
    659587         if(water.and.(nqmx.ge.2)) then
    660588           radius(igcm_h2o_ice)=3.e-6

  • trunk/LMDZ.MARS/libf/phymars/meso_inc/meso_inc_caps.F

    r226 r283  
     1
     2      !!!! see meso_inc_ini
     3      ! alb_lim = 0.26
     4      ! lat_lim = 70.
     5      ! inertie_lim = 800.
     6
     7      !!! This has to go after initracer which change dryness and watercaptag
     8
    19      !!!!!! MARS MESOSCALE MODELING
    210      !!!!!! TEST TEST TEST TEST  AS+JBM 28/02/11
    … …  
    1422      !!!!           du GCM lorsqu'ils sont consequents
    1523      !!!!
    16       IF ( caps .and. (igcm_h2o_ice .ne. 0) ) THEN
     24      !IF ( caps .and. (igcm_h2o_ice .ne. 0) ) THEN
     25      IF ( caps .and. water ) THEN
    1726          PRINT *, 'OVERWRITING watercaptag DEFINITION in INITRACER'
    18           PRINT *, 'lat>70 et alb>0.26 => watercaptag=T'
     27          PRINT *, 'lat > lat_lim et alb > alb_lim => watercaptag=T'
     28          PRINT *, 'ind for water ice: ', igcm_h2o_ice
    1929          !! Perennial H20 north cap defined by watercaptag=true (allows surface to be
    2030          !! hollowed by sublimation in vdifc).
    2131          do ig=1,ngridmx
    2232            qsurf(ig,igcm_h2o_ice)=0.  !! on jette les inputs GCM
    23             if ( (lati(ig)*180./pi.gt.70.) .and.
    24      .           (albedodat(ig).ge.0.26) )  then
    25                     watercaptag(ig)=.true.
    26                     dryness(ig) = 1.
     33            if ( ( lati(ig)*180./pi .gt. lat_lim ) .and.
     34     .           ( albedodat(ig) .ge. alb_lim    ) )  then
     35                    watercaptag(ig)  = .true.
     36                    dryness(ig)      = 1.
    2737            else
    28                     watercaptag(ig)=.false.
    29                     dryness(ig) = 1.
     38                    watercaptag(ig)  = .false.
     39                    dryness(ig)      = 1.
    3040            endif  ! (lati, albedodat)
    3141          end do ! (ngridmx)

  • trunk/LMDZ.MARS/libf/phymars/meso_inc/meso_inc_ini.F

    r226 r283  
    6868c
    6969ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
     70
     71      !!! see meso_inc_caps
     72      !!! this is a test to change outliers' albedo and thermal inertia
     73      alb_lim = 0.26
     74      lat_lim = 70.
     75      inertie_lim = 800.
     76      PRINT *, 'lat_lim ',lat_lim
     77      PRINT *, 'alb_lim ',alb_lim
     78      PRINT *, 'inertie_lim ',inertie_lim
     79      !!!
     80      !!!
     81      IF ( caps .and. water ) THEN
     82          do ig=1,ngridmx
     83            if ( lati(ig)*180./pi .gt. lat_lim ) then
     84               if ( albedodat(ig) .ge. alb_lim ) then
     85                    albedodat(ig) = alb_surfice
     86                    inertiedat(ig,1) = inertie_lim 
     87               endif
     88               if (inertiedat(ig,1) .ge. inertie_lim ) then
     89                    inertiedat(ig,1) = inertie_lim
     90               endif
     91            endif  ! (lati, albedodat)
     92          end do ! (ngridmx)
     93      ENDIF ! (caps)

  • trunk/LMDZ.MARS/libf/phymars/meso_inc/meso_inc_var.F

    r234 r283  
    1818      INTEGER tracerset    !!! this corresponds to config%mars
    1919      CHARACTER (len=20) :: wtnom(nqmx) ! tracer name
     20
     21      REAL alb_lim
     22      REAL lat_lim
     23      REAL inertie_lim

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

    r277 r283  
    128128
    129129#include "chimiedata.h"
    130 #include "watercap.h"
    131130#include "param.h"
    132131#include "param_v3.h"
    … …  
    191190      REAL fluxgrd(ngridmx)          ! surface conduction flux (W.m-2)
    192191      REAL qsurf(ngridmx,nqmx)       ! tracer on surface (e.g. kg.m-2)
    193       REAL q2(ngridmx,nlayermx+1)    ! Turbulent Kinetic Energy
     192      REAL q2(ngridmx,nlayermx+1)    ! Turbulent Kinetic Energy
     193     
     194      REAL watercapflag(ngridmx)     ! water cap flag
    194195
    195196c     Variables used by the water ice microphysical scheme:
    … …  
    197198      REAL nuice(ngridmx,nlayermx)   ! Estimated effective variance
    198199                                     !   of the size distribution
    199 c     Albedo of deposited surface ice
    200       !!REAL, PARAMETER :: alb_surfice = 0.4 ! 0.45
    201       REAL, PARAMETER :: alb_surfice = 0.45 !!TESTS_JB
    202200
    203201c     Variables used by the slope model
    … …  
    422420
    423421        IF (tracer.AND.water.AND.(ngridmx.NE.1)) THEN
    424           write(*,*)"physiq: water_param Surface ice alb:",alb_surfice
     422          write(*,*)"physiq: water_param Surface water ice albedo:",
     423     .                  albedo_h2o_ice
    425424        ENDIF
    426425                   
    … …  
    10951094#endif
    10961095c       -------------------------------------------------------------
    1097 c       Change of surface albedo (set to 0.4) in case of ground frost
     1096c       Change of surface albedo in case of ground frost
    10981097c       everywhere except on the north permanent cap and in regions
    10991098c       covered by dry ice.
    … …  
    11021101         do ig=1,ngrid
    11031102           if ((co2ice(ig).eq.0).and.
    1104      &        (qsurf(ig,igcm_h2o_ice).gt.0.005)) then
    1105               albedo(ig,1) = alb_surfice
    1106               albedo(ig,2) = alb_surfice
     1103     &        (qsurf(ig,igcm_h2o_ice).gt.frost_albedo_threshold)) then
     1104              albedo(ig,1) = albedo_h2o_ice
     1105              albedo(ig,2) = albedo_h2o_ice
     1106c              write(*,*) "frost thickness", qsurf(ig,igcm_h2o_ice)
     1107c              write(*,*) "physiq.F frost :"
     1108c     &        ,lati(ig)*180./pi, long(ig)*180./pi
    11071109           endif
    11081110         enddo  ! of do ig=1,ngrid
    … …  
    15121514     &                       'surface h2o_ice',
    15131515     &                       'kg.m-2',2,qsurf(1,igcm_h2o_ice))
     1516
     1517            if (caps) then
     1518             do ig=1,ngridmx
     1519                if (watercaptag(ig)) watercapflag(ig) = 1
     1520             enddo
     1521             CALL WRITEDIAGFI(ngridmx,'watercaptag',
     1522     &                         'Ice water caps',
     1523     &                         '',2,watercapflag)
     1524            endif
     1525            CALL WRITEDIAGFI(ngridmx,'albedo',
     1526     &                         'albedo',
     1527     &                         '',2,albedo(1:ngridmx,1))
    15141528           endif !(water)
    15151529

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

    r38 r283  
    1818
    1919#include"comsoil.h"
     20
     21#include"surfdat.h"
     22#include"callkeys.h"
    2023
    2124c-----------------------------------------------------------------------
    … …  
    5356! 0.1 Build mthermdiff(:), the mid-layer thermal diffusivities
    5457      do ig=1,ngrid
    55         do ik=0,nsoil-1
    56           mthermdiff(ig,ik)=therm_i(ig,ik+1)*therm_i(ig,ik+1)/volcapa
    57 !         write(*,*),'soil: ik: ',ik,' mthermdiff:',mthermdiff(ig,ik)
    58         enddo
     58        if (watercaptag(ig)) then
     59          do ik=0,nsoil-1
     60! If we have permanent ice, we use the water ice thermal inertia from ground to last layer.
     61              mthermdiff(ig,ik)=inert_h2o_ice*inert_h2o_ice/volcapa
     62          enddo
     63        else
     64          do ik=0,nsoil-1
     65           mthermdiff(ig,ik)=therm_i(ig,ik+1)*therm_i(ig,ik+1)/volcapa
     66          enddo
     67        endif
    5968      enddo
    6069

  • trunk/LMDZ.MARS/libf/phymars/surfdat.h

    r226 r283  
    55     &   zmea(ngridmx),zstd(ngridmx),                                   &
    66     &   zsig(ngridmx),zgam(ngridmx),zthe(ngridmx),                     &
    7      &   z0(ngridmx),z0_default
     7     &   z0(ngridmx),z0_default, albedo_h2o_ice, inert_h2o_ice,         &
     8     &   frost_albedo_threshold
    89
    9       COMMON/surfdatl/TESicealbedo
     10      COMMON/surfdatl/TESicealbedo,watercaptag,temptag
     11
    1012
    1113      real albedodat ! albedo of bare ground
    … …  
    1517      real emissiv ! emissivity of bare ground
    1618      logical TESicealbedo ! use TES ice cap albedoes (if set to .true.)
     19      logical watercaptag(ngridmx) ! flag for water ice surface
     20     
     21      logical temptag !temp tag for water caps
     22     
     23      real albedo_h2o_ice ! water ice albedo
     24      real inert_h2o_ice ! water ice thermal inertia
     25      real frost_albedo_threshold ! water frost thickness on the ground (kg.m^-2, ie mm)
    1726      real TESice_Ncoef ! coefficient for TES ice albedo in Northern hemisphere
    1827      real TESice_Scoef ! coefficient for TES ice albedo in Southern hemisphere

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

    r38 r283  
    1414#include "callkeys.h"
    1515#include "tracer.h"
    16 c
    17       INTEGER ngrid,ig,icap
     16#include "comgeomfi.h"
     17#include "comconst.h"
     18
     19c
     20      INTEGER ngrid,ig,icap,iq,alternate
    1821      REAL  piceco2(ngrid),psolaralb(ngrid,2)
    1922      REAL qsurf(ngrid,nqmx) !tracer on surface (kg/m2)
    … …  
    2427c=======================================================================
    2528
    26 c
    27 c     calcul de piceco2 (kg/m2) a l'etat initial
     29
     30c     water ice outliers
    2831c     ------------------------------------------
    2932
    30       DO 100 ig=1,ngrid
    31          psolaralb(ig,1)=albedodat(ig)
    32          psolaralb(ig,2)=albedodat(ig)
    33 100   CONTINUE
    34 
    35       PRINT*,'minimum des donnees albedo',
     33      IF ((water) .and. (caps)) THEN
     34     
     35c Perennial H20 north cap defined by watercaptag=true (allows surface to be
     36c hollowed by sublimation in vdifc).
     37
     38c We might not want albedodat to be modified because it is used to write
     39c restart files. Instead, albedo is directly modified when needed (i.e.
     40c if we have watercaptag and no co2 ice), below and in albedocaps.F90
     41
     42c       "Dryness coefficient" controlling the evaporation and
     43c        sublimation from the ground water ice (close to 1)
     44c        HERE, the goal is to correct for the fact
     45c        that the simulated permanent water ice polar caps
     46c        is larger than the actual cap and the atmospheric
     47c        opacity not always realistic.
     48
     49         alternate = 0
     50         
     51         do ig=1,ngridmx
     52         
     53         dryness (ig) = 1
     54         
     55
     56c Towards olympia planitia water caps ('relatively' low latitude ones)
     57c---------------- proposition par AS --------------------
     58c--------------------------------------------------------
     59c          if ( ( lati(ig)*180./pi      .ge.  75 ) .and.
     60c     .         ( lati(ig)*180./pi      .le.  77 ) .and.
     61c     .         ( ( ( long(ig)*180./pi .ge. 000. ) .and.
     62c     .             ( long(ig)*180./pi .le. 120. ) )
     63c     .           .or.
     64c     .           ( ( long(ig)*180./pi .ge. -130. ) .and.
     65c     .             ( long(ig)*180./pi .le. -115. ) ) ) ) then
     66c---------------- proposition par TN --------------------
     67c---------------- HIGHLY EXPERIMENTAL -------------------
     68c--------------------------------------------------------     
     69          if ( ( ( lati(ig)*180./pi .ge. 73.  ) .and. ! cap #1
     70     .           ( lati(ig)*180./pi .le. 75.1 ) .and.
     71     .           ( long(ig)*180./pi .ge. 95.  ) .and.
     72     .           ( long(ig)*180./pi .le. 110. ) )
     73     .          .or.
     74     .         ( ( lati(ig)*180./pi .ge. 77.  ) .and. ! cap #2
     75     .           ( lati(ig)*180./pi .le. 80.  ) .and.
     76     .           ( long(ig)*180./pi .ge. 110. ) .and.
     77     .           ( long(ig)*180./pi .le. 140. ) )
     78     .          .or.
     79     .         ( ( lati(ig)*180./pi .ge. 74.9 ) .and. ! cap #3
     80     .           ( lati(ig)*180./pi .le. 78.  ) .and.
     81     .           ( long(ig)*180./pi .ge. 155. ) .and.
     82     .           ( long(ig)*180./pi .le. 180. ) )
     83c     .         .or.
     84c     .         ( ( lati(ig)*180./pi .ge. 71.  ) .and. ! cap #4 (Korolev crater)
     85c     .           ( lati(ig)*180./pi .le. 72.  ) .and.
     86c     .           ( long(ig)*180./pi .ge. 163. ) .and.
     87c     .           ( long(ig)*180./pi .le. 168. ) )
     88     .          .or.
     89     .         ( ( lati(ig)*180./pi .ge. 74.9 ) .and. ! cap #5
     90     .           ( lati(ig)*180./pi .le. 78.  ) .and.
     91     .           ( long(ig)*180./pi .ge. -160.) .and.
     92     .           ( long(ig)*180./pi .le. -120.) ) )
     93     .          then
     94     
     95             if (temptag) then
     96             
     97               if ((alternate .eq. 0)) then  !!! 1/2 en 64x48 sinon trop large en lat
     98                 if (ngridmx.ne.1) watercaptag(ig)=.true.
     99                  write(*,*) "outliers ", lati(ig)*180./pi,
     100     .              long(ig)*180./pi
     101                !dryness(ig) = 1.
     102                 alternate = 1
     103               else
     104                 alternate = 0
     105               endif !end if alternate = 0
     106               
     107             else
     108             
     109               if (ngridmx.ne.1) watercaptag(ig)=.true.
     110                  write(*,*) "outliers ", lati(ig)*180./pi,
     111     .              long(ig)*180./pi
     112     
     113             endif ! end if temptag
     114             
     115           endif
     116
     117
     118c Opposite olympia planitia water cap
     119c---------------- proposition par AS --------------------
     120c--------------------------------------------------------
     121c           if ( ( lati(ig)*180./pi      .ge.  82 ) .and.
     122c     .          ( lati(ig)*180./pi      .le.  84 ) .and.
     123c     .          ( ( long(ig)*180./pi .gt. -030. ) .and.
     124c     .          ( long(ig)*180./pi .lt.  090. ) ) ) then
     125c---------------- proposition par TN --------------------
     126c--------------------------------------------------------
     127           if ( ( lati(ig)*180./pi     .ge.  80 ) .and.
     128     .          ( lati(ig)*180./pi     .le.  84 ) .and.
     129     .          ( ( long(ig)*180./pi .gt. -030. ) .and.
     130     .            ( long(ig)*180./pi .lt.  090. ) ) ) then
     131              if (ngridmx.ne.1) then
     132                watercaptag(ig)=.true.
     133                write(*,*) "central cap add ", lati(ig)*180./pi,
     134     .            long(ig)*180./pi
     135              endif
     136              !dryness(ig) = 1.
     137           endif
     138
     139c Central cap
     140c---------------- proposition par TN --------------------
     141c--------------------------------------------------------
     142           if (lati(ig)*180./pi.gt.84) then
     143             PRINT*,'central cap', lati(ig)*180./pi,
     144     .         long(ig)*180./pi
     145             if (ngridmx.ne.1) watercaptag(ig)=.true.
     146             !dryness(ig) = 1.
     147c Use the following cap definition for high spatial resolution (latitudinal bin <= 5 deg)
     148c             if (lati(ig)*180./pi.lt.85.and.long(ig).ge.0) then
     149c               if (ngridmx.ne.1) watercaptag(ig)=.true.
     150c               dryness(ig) = 1.
     151c             endif
     152c             if (lati(ig)*180./pi.ge.85) then
     153c               if (ngridmx.ne.1) watercaptag(ig)=.true.
     154c               dryness(ig) = 1.
     155c             endif
     156           endif  ! (lati>84 deg)
     157           
     158         end do ! (ngridmx)
     159        ENDIF ! (caps & water)
     160
     161c ===============================================================
     162c      INITIAL ALBEDO
     163c ===============================================================
     164
     165         write(*,*)"surfini: water frost thickness",
     166     s     frost_albedo_threshold
     167         write(*,*)"surfini: water ice albedo:", albedo_h2o_ice
     168         write(*,*)"surfini: water ice TI:", inert_h2o_ice
     169
     170c        To start with : Initial albedo = observed dataset
     171c        -------------------------------------------------
     172         DO ig=1,ngrid
     173              psolaralb(ig,1)=albedodat(ig)
     174              psolaralb(ig,2)=albedodat(ig)
     175         END DO
     176         PRINT*,'minimum albedo sans water caps',
    36177     s     albedodat(ISMIN(ngrid,albedodat,1))
    37       PRINT*,'maximum des donnees albedo',
     178         PRINT*,'maximum albedo sans water caps',
    38179     s     albedodat(ISMAX(ngrid,albedodat,1))
    39 c      calcul de psolaralb
    40 c      -------------------
    41       DO 115 ig=1,ngrid
    42 
    43 c        IF (water) THEN
    44 c          if (qsurf(ig,nqmx).gt.0.005) then
    45 c             psolaralb(ig,1) = 0.4
    46 c             psolaralb(ig,2) = 0.4
    47 c           endif
    48 c         ENDIF
    49 c IF there is more than 5 pr. um of h2o ice but no C02 ice, surface albedo is set to 0.4.
     180
     181c        initial albedo if permanent H2O ice is present
     182c        ------------------------------------------------
     183         IF ((water) .and. (caps)) THEN
     184           DO ig=1,ngrid
     185            IF (watercaptag(ig)) THEN
     186              psolaralb(ig,1) = albedo_h2o_ice
     187              psolaralb(ig,2) = albedo_h2o_ice
     188            ENDIF
     189           END DO
     190           PRINT*,'minimum albedo avec water caps',
     191     s     psolaralb(ISMIN(ngrid,psolaralb,1),1)
     192           PRINT*,'maximum albedo avec water caps',
     193     s     psolaralb(ISMAX(ngrid,psolaralb,1),1)
     194         ENDIF
     195
     196c      changing initial albedo if CO2 ice is present
     197c      -------------------------------------------
     198
     199       DO ig=1,ngrid
    50200         IF (piceco2(ig) .GT. 0.) THEN
    51201             IF(ig.GT.ngrid/2+1) THEN
    … …  
    55205             ENDIF
    56206             psolaralb(ig,1) = albedice(icap)
    57              psolaralb(ig,2) =  albedice(icap)
     207             psolaralb(ig,2) = albedice(icap)   
    58208         END IF
    59 115   CONTINUE     
    60 
    61       PRINT*,'minimum des donnees albedo',
     209       END DO
     210
     211c      changing initial albedo if water ice frost is present
     212c      -------------------------------------------
     213       IF (water) THEN
     214          do iq=1,nqmx
     215c          if there is frost and surface albedo is set to albedo_h2o_ice
     216           if(noms(iq).eq."h2o_ice") then
     217             do ig=1,ngrid
     218              if ((piceco2(ig) .eq. 0.).and.
     219     &          (qsurf(ig,iq).gt.frost_albedo_threshold)) then
     220                     psolaralb(ig,1) = albedo_h2o_ice
     221                     psolaralb(ig,2) = albedo_h2o_ice
     222c                     PRINT*,'surfini.F frost',
     223c     &                  lati(ig)*180./pi, long(ig)*180./pi
     224               endif
     225              enddo
     226           endif
     227          end do
     228          PRINT*,'minimum albedo avec givre et co2',
    62229     s     psolaralb(ISMIN(ngrid,psolaralb,1),1)
    63       PRINT*,'maximum des donnees albedo',
     230          PRINT*,'maximum albedo avec givre et co2',
    64231     s     psolaralb(ISMAX(ngrid,psolaralb,1),1)
     232       END IF
     233         
    65234
    66235      RETURN

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

    r246 r283  
    4040#include "netcdf.inc"
    4141#include "comg1d.h"
    42 #include "watercap.h"
    4342#include "logic.h"
    4443#include "advtrac.h"
    … …  
    572571
    573572      call physdem1("startfi.nc",long,lati,nsoilmx,nqmx,
    574      .              dtphys,float(day0),
    575      .              time,tsurf,tsoil,co2ice,emis,q2,qsurf,
    576      .              area,albedodat,inertiedat,zmea,zstd,zsig,zgam,zthe)
     573     .              dtphys,float(day0),time,tsurf,
     574     .              tsoil,co2ice,emis,q2,qsurf,area,albedodat,
     575     .              inertiedat,zmea,zstd,zsig,zgam,zthe)
    577576c=======================================================================
    578577c  BOUCLE TEMPORELLE DU MODELE 1D

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

    r268 r283  
    3636#include "tracer.h"
    3737
    38 #include "watercap.h"
    3938c
    4039c   arguments:
Note: See TracChangeset for help on using the changeset viewer.