Changeset 3329 for trunk/LMDZ.PLUTO


Ignore:
Timestamp:
May 15, 2024, 8:29:23 PM (6 months ago)
Author:
afalco
Message:

Pluto PCM:
Include cooling, hazes in radiative module
AF

Location:
trunk/LMDZ.PLUTO
Files:
6 added
11 edited

Legend:

Unmodified
Added
Removed

  • trunk/LMDZ.PLUTO/deftank/callphys.def

    r3197 r3329  
    8888
    8989naerkind = 1
    90 hazeprop = optprop_tholins_fractal100nm.dat
     90hazeprop_file = optprop_tholins_fractal100nm.dat
    9191
    9292## Other physics options

  • trunk/LMDZ.PLUTO/deftank/z2sig.def

    r3193 r3329  
    11  18.00000
    2   0.010   
    3   0.015   
     2  0.010
     3  0.015
    44  0.025
    55  0.040

  • trunk/LMDZ.PLUTO/libf/phypluto/aerosol_mod.F90

    r3195 r3329  
    7474      IF (firstcall) then
    7575        firstcall=.false.
    76         file_path=trim(datadir)//'/haze_prop/hazemmr.txt'
     76        file_path=trim(datadir)//'/haze_prop/'//hazemmr_file
    7777        open(224,file=file_path,form='formatted')
    7878        do ifine=1,Nfine

  • trunk/LMDZ.PLUTO/libf/phypluto/callcorrk.F90

    r3275 r3329  
    77      subroutine callcorrk(ngrid,nlayer,pq,nq,qsurf,           &
    88          albedo,albedo_equivalent,emis,mu0,pplev,pplay,pt,    &
    9           zzlay,tsurf,fract,dist_star,aerosol,muvar,                 &
     9          zzlay,tsurf,fract,dist_star,aerosol,muvar,           &
    1010          dtlw,dtsw,fluxsurf_lw,                               &
    1111          fluxsurf_sw,fluxsurfabs_sw,fluxtop_lw,               &
    … …  
    2626      use ioipsl_getin_p_mod, only: getin_p
    2727      use gases_h, only: ngasmx
    28       use radii_mod, only : su_aer_radii
     28      use radii_mod, only : su_aer_radii, haze_reffrad_fix
    2929      use aerosol_mod, only : iaero_haze
    3030      use aeropacity_mod, only: aeropacity
    … …  
    225225      parameter(Nfine=701)
    226226      real,save :: levdat(Nfine),vmrdat(Nfine)
     227      REAL dtlw_hcn_c2h2(ngrid, nlayer) ! cooling rate (K/s) due to C2H2/HCN (diagnostic)
    227228      real :: vmrch4(ngrid,nlayer)              ! vmr ch4 from vmrch4_proffix
    228229
    … …  
    506507
    507508      if (aerohaze) then
    508       !    if (haze_radproffix) then
    509       !       call haze_reffrad_fix(ngrid,nlayer,zzlay,&
    510       !           reffrad,nueffrad)
    511       !     endif
     509         if (haze_radproffix) then
     510            call haze_reffrad_fix(ngrid,nlayer,zzlay,&
     511                reffrad,nueffrad)
     512          endif
    512513
    513514         ! Get 3D aerosol optical properties.
    … …  
    554555    !   dtlw_hcn_c2h2=0.
    555556      if (cooling) then
    556         ! CALL cooling_hcn_c2h2(ngrid,nlayer,pplay,&
    557                                 !   pt,dtlw_hcn_c2h2)
     557        CALL cooling_hcn_c2h2(ngrid,nlayer,pplay,&
     558                                  pt,dtlw_hcn_c2h2)
    558559      endif
    559560
    … …  
    575576
    576577
    577          if (aerohaze) then
    578             do iaer=1,naerkind
    579                ! Shortwave.
    580                do nw=1,L_NSPECTV
    581 
    582                   do l=1,nlayer
    583 
    584                      temp1=QVISsQREF3d(ig,nlayer+1-l,nw,iaer)         &
    585                          *QREFvis3d(ig,nlayer+1-l,iaer)
    586 
    587                      temp2=QVISsQREF3d(ig,max(nlayer-l,1),nw,iaer)    &
    588                          *QREFvis3d(ig,max(nlayer-l,1),iaer)
    589 
    590                      qxvaer(2*l,nw,iaer)  = temp1
    591                      qxvaer(2*l+1,nw,iaer)=(temp1+temp2)/2
    592 
    593                      temp1=temp1*omegavis3d(ig,nlayer+1-l,nw,iaer)
    594                      temp2=temp2*omegavis3d(ig,max(nlayer-l,1),nw,iaer)
    595 
    596                      qsvaer(2*l,nw,iaer)  = temp1
    597                      qsvaer(2*l+1,nw,iaer)=(temp1+temp2)/2
    598 
    599                      temp1=gvis3d(ig,nlayer+1-l,nw,iaer)
    600                      temp2=gvis3d(ig,max(nlayer-l,1),nw,iaer)
    601 
    602                      gvaer(2*l,nw,iaer)  = temp1
    603                      gvaer(2*l+1,nw,iaer)=(temp1+temp2)/2
    604 
    605                   end do ! nlayer
    606 
    607                   qxvaer(1,nw,iaer)=qxvaer(2,nw,iaer)
    608                   qxvaer(2*nlayer+1,nw,iaer)=0.
    609 
    610                   qsvaer(1,nw,iaer)=qsvaer(2,nw,iaer)
    611                   qsvaer(2*nlayer+1,nw,iaer)=0.
    612 
    613                   gvaer(1,nw,iaer)=gvaer(2,nw,iaer)
    614                   gvaer(2*nlayer+1,nw,iaer)=0.
    615 
    616                end do ! L_NSPECTV
    617 
    618                do nw=1,L_NSPECTI
    619                   ! Longwave
    620                   do l=1,nlayer
    621 
    622                      temp1=QIRsQREF3d(ig,nlayer+1-l,nw,iaer)         &
    623                           *QREFir3d(ig,nlayer+1-l,iaer)
    624 
    625                      temp2=QIRsQREF3d(ig,max(nlayer-l,1),nw,iaer)    &
    626                           *QREFir3d(ig,max(nlayer-l,1),iaer)
    627 
    628                      qxiaer(2*l,nw,iaer)  = temp1
    629                      qxiaer(2*l+1,nw,iaer)=(temp1+temp2)/2
    630 
    631                      temp1=temp1*omegair3d(ig,nlayer+1-l,nw,iaer)
    632                      temp2=temp2*omegair3d(ig,max(nlayer-l,1),nw,iaer)
    633 
    634                      qsiaer(2*l,nw,iaer)  = temp1
    635                      qsiaer(2*l+1,nw,iaer)=(temp1+temp2)/2
    636 
    637                      temp1=gir3d(ig,nlayer+1-l,nw,iaer)
    638                      temp2=gir3d(ig,max(nlayer-l,1),nw,iaer)
    639 
    640                      giaer(2*l,nw,iaer)  = temp1
    641                      giaer(2*l+1,nw,iaer)=(temp1+temp2)/2
    642 
    643                   end do ! nlayer
    644 
    645                   qxiaer(1,nw,iaer)=qxiaer(2,nw,iaer)
    646                   qxiaer(2*nlayer+1,nw,iaer)=0.
    647 
    648                   qsiaer(1,nw,iaer)=qsiaer(2,nw,iaer)
    649                   qsiaer(2*nlayer+1,nw,iaer)=0.
    650 
    651                   giaer(1,nw,iaer)=giaer(2,nw,iaer)
    652                   giaer(2*nlayer+1,nw,iaer)=0.
    653 
    654                end do ! L_NSPECTI
    655 
    656             end do ! naerkind
    657 
    658             ! Test / Correct for freaky s. s. albedo values.
    659             do iaer=1,naerkind
    660                do k=1,L_LEVELS
    661 
    662                   do nw=1,L_NSPECTV
    663                      if(qsvaer(k,nw,iaer).gt.1.05*qxvaer(k,nw,iaer))then
    664                         message='Serious problems with qsvaer values'
    665                         call abort_physic(subname,message,1)
    666                      endif
    667                      if(qsvaer(k,nw,iaer).gt.qxvaer(k,nw,iaer))then
    668                         qsvaer(k,nw,iaer)=qxvaer(k,nw,iaer)
    669                      endif
    670                   end do
    671 
    672                   do nw=1,L_NSPECTI
    673                      if(qsiaer(k,nw,iaer).gt.1.05*qxiaer(k,nw,iaer))then
    674                         message='Serious problems with qsvaer values'
    675                         call abort_physic(subname,message,1)
    676                      endif
    677                      if(qsiaer(k,nw,iaer).gt.qxiaer(k,nw,iaer))then
    678                         qsiaer(k,nw,iaer)=qxiaer(k,nw,iaer)
    679                      endif
    680                   end do
    681 
    682                end do ! L_LEVELS
    683             end do ! naerkind
    684 
    685          endif ! aerohaze
     578        do iaer=1,naerkind
     579            ! Shortwave.
     580            do nw=1,L_NSPECTV
     581
     582                do l=1,nlayer
     583
     584                    temp1=QVISsQREF3d(ig,nlayer+1-l,nw,iaer)         &
     585                        *QREFvis3d(ig,nlayer+1-l,iaer)
     586
     587                    temp2=QVISsQREF3d(ig,max(nlayer-l,1),nw,iaer)    &
     588                        *QREFvis3d(ig,max(nlayer-l,1),iaer)
     589
     590                    qxvaer(2*l,nw,iaer)  = temp1
     591                    qxvaer(2*l+1,nw,iaer)=(temp1+temp2)/2
     592
     593                    temp1=temp1*omegavis3d(ig,nlayer+1-l,nw,iaer)
     594                    temp2=temp2*omegavis3d(ig,max(nlayer-l,1),nw,iaer)
     595
     596                    qsvaer(2*l,nw,iaer)  = temp1
     597                    qsvaer(2*l+1,nw,iaer)=(temp1+temp2)/2
     598
     599                    temp1=gvis3d(ig,nlayer+1-l,nw,iaer)
     600                    temp2=gvis3d(ig,max(nlayer-l,1),nw,iaer)
     601
     602                    gvaer(2*l,nw,iaer)  = temp1
     603                    gvaer(2*l+1,nw,iaer)=(temp1+temp2)/2
     604
     605                end do ! nlayer
     606
     607                qxvaer(1,nw,iaer)=qxvaer(2,nw,iaer)
     608                qxvaer(2*nlayer+1,nw,iaer)=0.
     609
     610                qsvaer(1,nw,iaer)=qsvaer(2,nw,iaer)
     611                qsvaer(2*nlayer+1,nw,iaer)=0.
     612
     613                gvaer(1,nw,iaer)=gvaer(2,nw,iaer)
     614                gvaer(2*nlayer+1,nw,iaer)=0.
     615
     616            end do ! L_NSPECTV
     617
     618            do nw=1,L_NSPECTI
     619                ! Longwave
     620                do l=1,nlayer
     621
     622                    temp1=QIRsQREF3d(ig,nlayer+1-l,nw,iaer)         &
     623                        *QREFir3d(ig,nlayer+1-l,iaer)
     624
     625                    temp2=QIRsQREF3d(ig,max(nlayer-l,1),nw,iaer)    &
     626                        *QREFir3d(ig,max(nlayer-l,1),iaer)
     627
     628                    qxiaer(2*l,nw,iaer)  = temp1
     629                    qxiaer(2*l+1,nw,iaer)=(temp1+temp2)/2
     630
     631                    temp1=temp1*omegair3d(ig,nlayer+1-l,nw,iaer)
     632                    temp2=temp2*omegair3d(ig,max(nlayer-l,1),nw,iaer)
     633
     634                    qsiaer(2*l,nw,iaer)  = temp1
     635                    qsiaer(2*l+1,nw,iaer)=(temp1+temp2)/2
     636
     637                    temp1=gir3d(ig,nlayer+1-l,nw,iaer)
     638                    temp2=gir3d(ig,max(nlayer-l,1),nw,iaer)
     639
     640                    giaer(2*l,nw,iaer)  = temp1
     641                    giaer(2*l+1,nw,iaer)=(temp1+temp2)/2
     642
     643                end do ! nlayer
     644
     645                qxiaer(1,nw,iaer)=qxiaer(2,nw,iaer)
     646                qxiaer(2*nlayer+1,nw,iaer)=0.
     647
     648                qsiaer(1,nw,iaer)=qsiaer(2,nw,iaer)
     649                qsiaer(2*nlayer+1,nw,iaer)=0.
     650
     651                giaer(1,nw,iaer)=giaer(2,nw,iaer)
     652                giaer(2*nlayer+1,nw,iaer)=0.
     653
     654            end do ! L_NSPECTI
     655
     656        end do ! naerkind
     657
     658        ! Test / Correct for freaky s. s. albedo values.
     659        do iaer=1,naerkind
     660            do k=1,L_LEVELS
     661
     662                do nw=1,L_NSPECTV
     663                    if(qsvaer(k,nw,iaer).gt.1.05*qxvaer(k,nw,iaer))then
     664                    message='Serious problems with qsvaer values'
     665                    call abort_physic(subname,message,1)
     666                    endif
     667                    if(qsvaer(k,nw,iaer).gt.qxvaer(k,nw,iaer))then
     668                    qsvaer(k,nw,iaer)=qxvaer(k,nw,iaer)
     669                    endif
     670                end do
     671
     672                do nw=1,L_NSPECTI
     673                    if(qsiaer(k,nw,iaer).gt.1.05*qxiaer(k,nw,iaer))then
     674                    message='Serious problems with qsvaer values'
     675                    call abort_physic(subname,message,1)
     676                    endif
     677                    if(qsiaer(k,nw,iaer).gt.qxiaer(k,nw,iaer))then
     678                    qsiaer(k,nw,iaer)=qxiaer(k,nw,iaer)
     679                    endif
     680                end do
     681
     682            end do ! L_LEVELS
     683        end do ! naerkind
     684
    686685!-----------------------------------------------------------------------
    687686!     Aerosol optical depths
    688687!-----------------------------------------------------------------------
    689         IF (aerohaze) THEN
    690          do iaer=1,naerkind     ! a bug was here
     688        do iaer=1,naerkind     ! a bug was here
    691689            do k=0,nlayer-1
    692690
    … …  
    706704
    707705         end do ! naerkind
    708         ELSE
    709          tauaero(:,:)=0
    710         ENDIF ! aerohaze
    711706
    712707         ! Albedo and Emissivity.
    … …  
    727722
    728723      if (generic_condensation .or. methane .or. carbox) then
     724
     725        ! IF (methane) then
     726
     727        !    do l=1,nlayer
     728        !        qvar(2*l)   = vmrch4(ig,nlayer+1-l)/100.*    &
     729        !                             mmol(igcm_ch4_gas)/mmol(igcm_n2)
     730        !        qvar(2*l+1) = ((vmrch4(ig,nlayer+1-l)+vmrch4(ig, &
     731        !                     max(nlayer-l,1)))/2.)/100.* &
     732        !                     mmol(igcm_ch4_gas)/mmol(igcm_n2)
     733        !    end do
     734        !    qvar(1)=qvar(2)
     735
     736        ! ELSE
    729737
    730738         ! For now, only one GCS tracer can be both variable and radiatively active

  • trunk/LMDZ.PLUTO/libf/phypluto/datafile_mod.F90

    r3193 r3329  
    1212      character(len=300),save :: datadir='/u/lmdz/WWW/planets/LMDZ.GENERIC/datagcm'
    1313!$OMP THREADPRIVATE(datadir)
    14      
     14
    1515      ! Subdirectories of 'datadir':
    16      
     16
    1717      ! surfdir stores planetary topography, albedo, etc. (surface.nc files)
    1818      character(len=12),parameter :: surfdir="surface_data"
    19      
     19
    2020      ! aerdir stores aerosol properties files (optprop_*dat files)
    21       character(LEN=18),parameter :: aerdir="aerosol_properties" 
     21      character(LEN=18),parameter :: aerdir="aerosol_properties"
    2222
    2323      ! Data haze properties
    24       character(len=300),save :: hazeprop
    25      
     24      character(len=300),save :: hazeprop_file
     25      character(len=300),save :: hazerad_file
     26      character(len=300),save :: hazemmr_file
     27
    2628      end module datafile_mod
    2729!-----------------------------------------------------------------------

  • trunk/LMDZ.PLUTO/libf/phypluto/dyn1d/rcm1d.F

    r3258 r3329  
    2222      use time_phylmdz_mod, only: daysec, dtphys, day_step, ecritphy,
    2323     &                            nday, iphysiq
    24       use callkeys_mod, only: tracer, specOLR,pceil
     24      use callkeys_mod, only: tracer, specOLR,pceil,haze
    2525      USE comvert_mod, ONLY: ap,bp,aps,bps,pa,preff, sig,
    2626     &                       presnivs,pseudoalt,scaleheight
    … …  
    7878c
    7979      INTEGER day0          ! date initial (sol ; =0 a Ls=0)
     80      INTEGER lecttsoil     ! lecture of tsoil from proftsoil
     81      INTEGER lecthaze      ! lecture of haze from profhaze
    8082      REAL day              ! date durant le run
    8183      REAL time             ! time (0<time<1 ; time=0.5 a midi)
    … …  
    329331            endif
    330332          enddo !of do iq=1,nq
    331 ! check for n2 / h2o_ice tracers:
     333! check for n2 / ch4_ice tracers:
    332334         i_n2=0
    333335         do iq=1,nq
    … …  
    342344           elseif (tname(iq)=="co_gas") then
    343345             i_co_gas=iq
     346           elseif (tname(iq)=="haze") then
     347             i_haze=iq
    344348           endif
    345349         enddo
    … …  
    657661                q(ilayer,iq) = 0.01*coref*28./28.
    658662                ENDDO
    659             ! else if ((iq.eq.i_haze10).or.(iq.eq.i_haze30).or.(iq.eq.i_haze).or.(iq.eq.i_haze50).or.(iq.eq.i_haze100)) then
    660             !     hazeref=0. ! default value for haze mmr
    661             !     PRINT *,'qhaze (kg/kg) ?'
    662             !     call getin("hazeref",hazeref)
    663             !     write(*,*) " haze ref (kg/kg) = ",hazeref
    664             !     DO ilayer=1,nlayer
    665             !         q(ilayer,iq) = hazeref
    666             !     ENDDO
    667             ! else if (iq.eq.i_prec_haze) then
    668             !     prechazeref=0. ! default value for vmr ch4
    669             !     PRINT *,'qprechaze (kg/kg) ?'
    670             !     call getin("prechazeref",prechazeref)
    671             !     write(*,*) " prec haze ref (kg/kg) = ",prechazeref
    672             !     DO ilayer=1,nlayer
    673             !         q(ilayer,iq) = prechazeref
    674             !     ENDDO
    675 
     663            else if ((iq.eq.i_haze10).or.(iq.eq.i_haze30)
     664     &               .or.(iq.eq.i_haze).or.(iq.eq.i_haze50)
     665     &               .or.(iq.eq.i_haze100)) then
     666                hazeref=0. ! default value for haze mmr
     667                PRINT *,'qhaze (kg/kg) ?'
     668                call getin("hazeref",hazeref)
     669                write(*,*) " haze ref (kg/kg) = ",hazeref
     670                DO ilayer=1,nlayer
     671                    q(ilayer,iq) = hazeref
     672                ENDDO
     673            else if (iq.eq.i_prec_haze) then
     674                prechazeref=0. ! default value for vmr ch4
     675                PRINT *,'qprechaze (kg/kg) ?'
     676                call getin("prechazeref",prechazeref)
     677                write(*,*) " prec haze ref (kg/kg) = ",prechazeref
     678                DO ilayer=1,nlayer
     679                    q(ilayer,iq) = prechazeref
     680                ENDDO
    676681            else
    677682                DO ilayer=1,nlayer
    … …  
    856861! ------------------------------------------
    857862      volcapa=1.e6 ! volumetric heat capacity
    858       DO isoil=1,nsoil
     863      lecttsoil=0 ! default value for lecttsoil
     864      call getin("lecttsoil",lecttsoil)
     865      if (lecttsoil==1) then
     866         OPEN(14,file='proftsoil',status='old',form='formatted',err=101)
     867         DO isoil=1,nsoil
     868            READ (14,*) tsoil(isoil)
     869            inertiedat(1,isoil)=inertiedat(1,1) ! soil thermal inertia
     870         ENDDO
     871         GOTO 201
     872101      STOP'fichier proftsoil inexistant'
     873201      CONTINUE
     874         CLOSE(14)
     875
     876      else
     877        DO isoil=1,nsoil
    859878         inertiedat(1,isoil)=inertiedat(1,1) ! soil thermal inertia
    860879         tsoil(isoil)=tsurf(1)  ! soil temperature
    861       ENDDO
     880        ENDDO
     881      endif
     882
    862883
    863884! Initialize depths
    … …  
    870891      enddo
    871892
     893!     Initialize haze profile
     894!     ------------------------------------------
     895      if (haze) then
     896
     897      lecthaze=0 ! default value for lecthaze
     898      call getin("lecthaze",lecthaze)
     899      if (lecthaze==1) then
     900         OPEN(15,file='profhaze',status='old',form='formatted',err=301)
     901         DO iq=1,nq
     902            if (iq.eq.i_haze) then
     903            DO ilayer=1,nlayer
     904               READ (15,*) q(ilayer,iq)
     905            ENDDO
     906            endif
     907         ENDDO
     908         GOTO 401
     909301      STOP'fichier profhaze inexistant'
     910401      CONTINUE
     911         CLOSE(15)
     912      endif
     913      endif
    872914! Initialize cloud fraction and oceanic ice !AF24: removed
    873915! Initialize slab ocean !AF24: removed
    … …  
    10221064         CLOSE(12)
    10231065      endif
    1024 
     1066! save haze profile
     1067      if (haze.and.lecthaze.eq.1) then
     1068            OPEN(16,file='profhaze.out',form='formatted')
     1069            DO iq=1,nq
     1070             if (iq.eq.i_haze) then
     1071               DO ilayer=1,nlayer
     1072                 write(16,*) q(ilayer,iq)
     1073               ENDDO
     1074             endif
     1075            ENDDO
     1076            CLOSE(16)
     1077         endif
    10251078
    10261079      ENDDO   ! fin de la boucle temporelle

  • trunk/LMDZ.PLUTO/libf/phypluto/inifis_mod.F90

    r3275 r3329  
    1313  use radcommon_h, only: ini_radcommon_h
    1414  use radii_mod, only: radfixed, Nmix_n2
    15   use datafile_mod, only: datadir, hazeprop
     15  use datafile_mod, only: datadir, hazeprop_file, hazerad_file, hazemmr_file
    1616  use comdiurn_h, only: sinlat, coslat, sinlon, coslon
    1717  use comgeomfi_h, only: totarea, totarea_planet
    … …  
    138138
    139139     if (is_master) write(*,*) "Haze optical properties datafile"
    140      hazeprop="../../../datadir"  ! default path
    141      call getin_p("hazeprop",hazeprop)
    142      if (is_master) write(*,*) trim(rname)//" hazeprop = ",trim(hazeprop)
     140     hazeprop_file="optprop_tholins_fractal050nm.dat"  ! default file
     141     call getin_p("hazeprop_file",hazeprop_file)
     142     if (is_master) write(*,*) trim(rname)//" hazeprop_file = ",trim(hazeprop_file)
     143
     144     if (is_master) write(*,*) "Haze radii datafile"
     145     hazerad_file="hazerad.txt"  ! default file
     146     call getin_p("hazerad_file",hazerad_file)
     147     if (is_master) write(*,*) trim(rname)//" hazerad_file = ",trim(hazerad_file)
     148
     149     if (is_master) write(*,*) "Haze mmr datafile"
     150     hazemmr_file="hazemmr.txt"  ! default file
     151     call getin_p("hazemmr_file",hazemmr_file)
     152     if (is_master) write(*,*) trim(rname)//" hazemmr_file = ",trim(hazemmr_file)
    143153
    144154

  • trunk/LMDZ.PLUTO/libf/phypluto/physiq_mod.F90

    r3275 r3329  
    5454                              fast,fasthaze,haze,metcloud,monoxcloud,&
    5555                              n2cond,nearn2cond,noseason_day,conservn2, &
    56                               kbo,triton,paleo,paleoyears, &
     56                              convergeps,kbo,triton,paleo,paleoyears, &
    5757                              carbox, methane, oldplutovdifc, oldplutocorrk, &
    5858                              aerohaze,haze_proffix,source_haze, &
    … …  
    6767      use phys_state_var_mod
    6868      use callcorrk_mod, only: callcorrk
    69     !   use callcorrk_pluto_mod, only: callcorrk_pluto
     69      use callcorrk_pluto_mod, only: callcorrk_pluto
    7070      use vdifc_mod, only: vdifc
    7171      use vdifc_pluto_mod, only: vdifc_pluto
    … …  
    260260      REAL,save  :: tcond1p4Pa
    261261      DATA tcond1p4Pa/38/
     262      real :: tidat(ngrid,nsoilmx)    ! thermal inertia soil
    262263
    263264
    … …  
    306307      real zlss                      ! Sub solar point longitude (radians).
    307308      real zday                      ! Date (time since Ls=0, calculated in sols).
     309      REAL,save :: saveday           ! saved date
     310      REAL,save :: savedeclin        ! saved declin
    308311      real zzlay(ngrid,nlayer)       ! Altitude at the middle of the atmospheric layers.
    309312      real zzlev(ngrid,nlayer+1)     ! Altitude at the atmospheric layer boundaries.
    … …  
    619622         call iniorbit(apoastr,periastr,year_day,peri_day,obliquit)
    620623
     624         savedeclin=0.
     625         saveday=pday
     626         !adjust=0. ! albedo adjustment for convergeps
    621627
    622628!        Initialize soil.
    … …  
    871877         if ((ngrid.eq.1).and.(global1d)) then ! Fixed zenith angle 'szangle' in 1D simulations w/ globally-averaged sunlight.
    872878            mu0 = cos(pi*szangle/180.0)
     879            fract= 1/(4*mu0) ! AF24: from pluto.old
    873880         endif
    874881
    875882      endif
    876883
    877       ! AF24: TODO insert surfprop for pluto & triton around here
     884
     885!     Pluto albedo /IT changes depending on surface ices (only in 3D)
     886!     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     887      if (ngrid.ne.1) then
     888
     889         !! Specific to change albedo of N2 so that Psurf
     890         !! converges toward 1.4 Pa in "1989" seasons for Triton
     891         !! converges toward 1.1 Pa in "2015" seasons for Pluto
     892         if (convergeps) then
     893            if (triton) then
     894               ! 1989 declination
     895               if (declin*180./pi.gt.-46..and.declin*180./pi.lt.-45.   &
     896            .and.zday.gt.saveday+1000.   &
     897            .and.declin.lt.savedeclin) then
     898               call globalaverage2d(ngrid,pplev(:,1),globave)
     899               if (globave.gt.1.5) then
     900                     adjust=adjust+0.005
     901               else if (globave.lt.1.3) then
     902                     adjust=adjust-0.005
     903               endif
     904               saveday=zday
     905               endif
     906            else
     907               ! Pluto : 2015 declination current epoch
     908               if (declin*180./pi.gt.50.and.declin*180./pi.lt.51.  &
     909            .and.zday.gt.saveday+10000.  &
     910            .and.declin.gt.savedeclin) then
     911               call globalaverage2d(ngrid,pplev(:,1),globave)
     912               if (globave.gt.1.2) then
     913                     adjust=adjust+0.005
     914               else if (globave.lt.1.) then
     915                     adjust=adjust-0.005
     916               endif
     917               saveday=zday
     918               endif
     919            endif
     920         end if
     921      end if ! if ngrid ne 1
     922
     923      call surfprop(ngrid,nq,fract,qsurf,tsurf,tidat,   &
     924      capcal,adjust,dist_star,albedo,emis,flusurfold,ptimestep,   &
     925      zls)
     926      ! do i=2,ngrid
     927      !    albedo(i,:) = albedo(1,:)
     928      ! enddo
    878929
    879930      if (callrad) then
    … …  
    920971
    921972               ! standard callcorrk
    922                ! clearsky=.false.
    923             !    if (oldplutocorrk) then
    924             !       call callcorrk_pluto(icount,ngrid,nlayer,pq,nq,qsurf,          &
    925             !                    albedo,emis,mu0,pplev,pplay,pt,                   &
    926             !                    tsurf,fract,dist_star,aerosol, &
    927             !                    zdtlw,zdtsw,fluxsurf_lw,fluxsurf_sw,fluxtop_lw,  &
    928             !                    fluxabs_sw,fluxtop_dn,reffrad,tau_col,ptime,pday, &
    929             !                    firstcall,lastcall,zzlay)
    930             !    else
    931                 call callcorrk(ngrid,nlayer,pq,nq,qsurf,                          &
     973               if (oldplutocorrk) then
     974                  call callcorrk_pluto(icount,ngrid,nlayer,pq,nq,qsurf,          &
     975                               albedo(:,1),emis,mu0,pplev,pplay,pt,                   &
     976                               zzlay,tsurf,fract,dist_star,aerosol,              &
     977                               zdtlw,zdtsw,fluxsurf_lw,fluxsurf_sw,fluxtop_lw,  &
     978                               fluxabs_sw,fluxtop_dn,reffrad,tau_col,ptime,pday, &
     979                               cloudfrac,totcloudfrac,.false.,                  &
     980                               firstcall,lastcall)
     981               else
     982                call callcorrk(ngrid,nlayer,pq,nq,qsurf,  &
    932983                              albedo,albedo_equivalent,emis,mu0,pplev,pplay,pt,   &
    933984                              zzlay,tsurf,fract,dist_star,aerosol,muvar,                &
    … …  
    938989                              tau_col,cloudfrac,totcloudfrac,                     &
    939990                              .false.,firstcall,lastcall)
    940             !    endif ! oldplutocorrk
     991               endif ! oldplutocorrk
    941992                !GG (feb2021): Option to "artificially" decrease the raditive time scale in
    942993                !the deep atmosphere  press > 0.1 bar. Suggested by MT.
    … …  
    9561007
    9571008                ! AF24: removed CLFvarying Option
    958 
    959                ! if(ok_slab_ocean) then
    960                !    tsurf(:)=tsurf2(:)
    961                ! endif
    962 
    9631009
    9641010               ! Radiative flux from the sky absorbed by the surface (W.m-2).
    … …  
    10521098
    10531099         if (oldplutovdifc) then
    1054             zflubid(1:ngrid)=fluxrad(1:ngrid)+fluxgrd(1:ngrid)
    10551100            zdum1(:,:) = 0
    10561101            zdum2(:,:) = 0
    … …  
    10661111                    zdqdif,zdqsdif,qsat_ch4,qsat_ch4_l1) !,zq1temp_ch4,qsat_ch4)
    10671112
    1068             pdv(1:ngrid,1:nlayer)=pdv(1:ngrid,1:nlayer)+zdvdif(1:ngrid,1:nlayer)
    1069             pdu(1:ngrid,1:nlayer)=pdu(1:ngrid,1:nlayer)+zdudif(1:ngrid,1:nlayer)
    1070             pdt(1:ngrid,1:nlayer)=pdt(1:ngrid,1:nlayer)+zdhdif(1:ngrid,1:nlayer)*zpopsk(1:ngrid,1:nlayer)
    1071 
    10721113            zdtdif(1:ngrid,1:nlayer)=zdhdif(1:ngrid,1:nlayer)*zpopsk(1:ngrid,1:nlayer) ! for diagnostic only
    1073             zdtsurf(1:ngrid)=zdtsurf(1:ngrid)+zdtsdif(1:ngrid)
    10741114
    10751115            bcond=1./tcond1p4Pa
    10761116            acond=r/lw_n2
    1077 
    1078             if (tracer) then
    1079                pdq(1:ngrid,1:nlayer,1:nq)=pdq(1:ngrid,1:nlayer,1:nq)+ zdqdif(1:ngrid,1:nlayer,1:nq)
    1080                dqsurf(1:ngrid,1:nq)=dqsurf(1:ngrid,1:nq) + zdqsdif(1:ngrid,1:nq)
    1081             end if ! of if (tracer)
    10821117
    10831118            !------------------------------------------------------------------
    … …  
    13791414         zdqphot_prec(:,:)=0.
    13801415         zdqphot_ch4(:,:)=0.
     1416         zdqhaze(:,:,:)=0
    13811417         ! Forcing to a fixed haze profile if haze_proffix
    13821418         if (haze_proffix.and.i_haze.gt.0.) then
    1383          call haze_prof(ngrid,nlayer,zzlay,pplay,pt,  &
    1384                                   reffrad,profmmr)
    1385          zdqhaze(:,:,i_haze)=(profmmr(:,:)-pq(:,:,igcm_haze))  &
    1386                                                   /ptimestep
     1419            call haze_prof(ngrid,nlayer,zzlay,pplay,pt,  &
     1420                           reffrad,profmmr)
     1421            zdqhaze(:,:,i_haze)=(profmmr(:,:)-pq(:,:,igcm_haze))  &
     1422                                 /ptimestep
    13871423         else
    1388          call hazecloud(ngrid,nlayer,nq,ptimestep, &
    1389             pplay,pplev,pq,pdq,dist_star,mu0,zfluxuv,zdqhaze,   &
    1390             zdqphot_prec,zdqphot_ch4,zdqconv_prec,declin)
     1424            call hazecloud(ngrid,nlayer,nq,ptimestep, &
     1425               pplay,pplev,pq,pdq,dist_star,mu0,zfluxuv,zdqhaze,   &
     1426               zdqphot_prec,zdqphot_ch4,zdqconv_prec,declin)
    13911427         endif
    13921428
    … …  
    15891625      ! todo: should be placed in haze and use tendency of n2 instead of flusurf
    15901626      IF (source_haze) THEN
     1627         write(*,*) "Source haze not supported yet."
     1628         stop
    15911629            !  call hazesource(ngrid,nlayer,nq,ptimestep,  &
    15921630            !                 pplev,flusurf,mu0,zdq_source)

  • trunk/LMDZ.PLUTO/libf/phypluto/radii_mod.F90

    r3195 r3329  
    2828      use ioipsl_getin_p_mod, only: getin_p
    2929      use radinc_h, only: naerkind
     30      use datafile_mod, only: hazerad_file
    3031      use aerosol_mod, only: iaero_haze, i_haze, &
    3132                              iaero_generic, i_rgcs_ice
    … …  
    5051           !reffrad(1:ngrid,1:nlayer,iaer) = 2.e-6 ! haze
    5152           nueffrad(1:ngrid,1:nlayer,iaer) = 0.02 ! haze
    52            print*, 'TB22 Hello2',radius(i_haze)*nmono**(1./3.)
    5353         endif
    5454      enddo
    … …  
    100100      IF (firstcall) then
    101101         firstcall=.false.
    102          file_path=trim(datadir)//'/haze_prop/hazerad.txt'
     102         file_path=trim(datadir)//'/haze_prop/'//hazerad_file
    103103         open(223,file=file_path,form='formatted')
    104104         do ifine=1,Nfine

  • trunk/LMDZ.PLUTO/libf/phypluto/suaer_corrk.F90

    r3195 r3329  
    1010      use radinc_h,    only: L_NSPECTI,L_NSPECTV,nsizemax,naerkind
    1111      use radcommon_h, only: blamv,blami,lamrefir,lamrefvis
    12       use datafile_mod, only: datadir, aerdir, hazeprop
     12      use datafile_mod, only: datadir, aerdir, hazeprop_file
    1313
    1414      ! outputs
    … …  
    142142         ! Only one table of optical properties :
    143143         write(*,*)'Suaer haze optical properties, using: ', &
    144                            TRIM(hazeprop)
     144                           TRIM(hazeprop_file)
    145145
    146146         ! visible
    147          file_id(iaer,1)=TRIM(hazeprop)
     147         file_id(iaer,1)=TRIM(hazeprop_file)
    148148         ! infrared
    149149         file_id(iaer,2)=file_id(iaer,1)
    … …  
    155155      enddo
    156156
    157       IF (naerkind .gt. 1) THEN
    158          print*,'naerkind = ',naerkind
    159          print*,'but we only have data for 1 type, exiting.'
    160          call abort
    161       ENDIF
     157      ! IF (naerkind .gt. 1) THEN ! AF24: ?
     158      !    print*,'naerkind = ',naerkind
     159      !    print*,'but we only have data for 1 type, exiting.'
     160      !    call abort
     161      ! ENDIF
    162162
    163163!------------------------------------------------------------------
    … …  
    314314
    315315      jfile = jfile+1
    316       IF ((idomain.NE.iaero_haze).OR.(iaer.NE.iaero_haze)) THEN
     316      print *, idomain, " ", iaer, " ", iaero_haze
     317      IF ((idomain.NE.1).OR.(iaer.NE.1).OR.(jfile.EQ.6)) THEN
     318         ! AF24: what does it mean? :-O
    317319         endwhile = .true.
    318320      ENDIF

  • trunk/LMDZ.PLUTO/libf/phypluto/surfprop.F90

    r3275 r3329  
    44
    55    !   use comgeomfi_h, only:
     6      use radinc_h, only: L_NSPECTV, L_NSPECTI
    67      use comcstfi_mod, only: pi
    78      use comsoil_h, only: nsoilmx,mlayer,volcapa,inertiedat
    … …  
    6465      REAL,INTENT(IN) :: adjust
    6566      REAL,INTENT(IN) :: dist
    66       REAL,INTENT(OUT) :: albedo(ngrid)
     67      REAL,INTENT(OUT) :: albedo(ngrid,L_NSPECTV)
    6768      REAL,INTENT(OUT) :: emis(ngrid)
    6869      REAL,INTENT(OUT) :: tidat(ngrid,nsoilmx)
    … …  
    169170
    170171                 CASE (0) ! fixed albedo
    171                      albedo(ig)=min(max(alb_n2b+adjust,0.),1.)
     172                     albedo(ig,:)=min(max(alb_n2b+adjust,0.),1.)
    172173
    173174                 CASE (1) ! Albedo decreases with thickness
    174                      albedo(ig)=(alb_n2b-deltab)/(1.-10000.)*qsurf(ig,igcm_n2)+alb_n2b
    175                      albedo(ig)=min(max(alb_n2b-deltab,albedo(ig)),alb_n2b)   ! should not be higher than alb_n2b and not lower than alb_n2b-deltab
     175                     albedo(ig,:)=(alb_n2b-deltab)/(1.-10000.)*qsurf(ig,igcm_n2)+alb_n2b
     176                     albedo(ig,:)=min(max(alb_n2b-deltab,albedo(ig,:)),alb_n2b)   ! should not be higher than alb_n2b and not lower than alb_n2b-deltab
    176177                 CASE (2) ! Special Sputnik differences of albedo
    177178                     if (qsurf(ig,igcm_n2).ge.1.e4.and.(longitude(ig)*180./pi.le.-161.or.longitude(ig)*180./pi.ge.146)) then
    … …  
    179180                                ((longitude(ig)*180./pi.gt.140.).and. &
    180181                                 (longitude(ig)*180./pi.lt.165.)) ) then
    181                                          albedo(ig)=alb_n2b-deltab
     182                                         albedo(ig,:)=alb_n2b-deltab
    182183                         else
    183                                          albedo(ig)=alb_n2b+deltab
     184                                         albedo(ig,:)=alb_n2b+deltab
    184185                         endif
    185186                     else
    186                              albedo(ig)=alb_n2b
     187                             albedo(ig,:)=alb_n2b
    187188                     endif
    188189
    189190                 CASE (3) ! Albedo increases (delta neg) or decreases (delta pos) with sublimation rates
    190                      albedo(ig)=deltab/1.e-4*fluold(ig,igcm_n2)+albedo(ig)
    191                      albedo(ig)=min(max(alb_n2b-deltab,albedo(ig)),alb_n2b+deltab)
     191                     albedo(ig,:)=deltab/1.e-4*fluold(ig,igcm_n2)+albedo(ig,:)
     192                     albedo(ig,:)=min(max(alb_n2b-deltab,albedo(ig,:)),alb_n2b+deltab)
    192193                     !! Triton:
    193                      !albedo(ig)=deltab/1.e-4*fluold(ig,igcm_n2)+albedo(ig)
    194                      !albedo(ig)=min(max(alb_n2b-deltab,albedo(ig)),alb_n2b+deltab)
     194                     !albedo(ig,:)=deltab/1.e-4*fluold(ig,igcm_n2)+albedo(ig,:)
     195                     !albedo(ig,:)=min(max(alb_n2b-deltab,albedo(ig,:)),alb_n2b+deltab)
    195196
    196197                 CASE (4) ! Albedo Difference in N/S
    197198                     if (latitude(ig)*180./pi.ge.0.) then
    198                         albedo(ig)=min(max(alb_n2b-deltab+adjust,0.),1.)
     199                        albedo(ig,:)=min(max(alb_n2b-deltab+adjust,0.),1.)
    199200                     else
    200                         albedo(ig)=min(max(alb_n2b+deltab+adjust,0.),1.)
     201                        albedo(ig,:)=min(max(alb_n2b+deltab+adjust,0.),1.)
    201202                     endif
    202203
    … …  
    210211                     ! .or.  ((latitude(ig)*180./pi.le.6.).and.(latitude(ig)*180./pi.ge.2.))   &
    211212                     !       ) then
    212                      !                    albedo(ig)=alb_n2b-deltab
     213                     !                    albedo(ig,:)=alb_n2b-deltab
    213214
    214215                     !! Patches GCM
    … …  
    218219                         .or. ((latitude(ig)*180./pi.lt.30.5).and.(latitude(ig)*180./pi.gt.29.5).and. &
    219220                              (longitude(ig)*180./pi.lt.169.).and.(longitude(ig)*180./pi.gt.168.))) then
    220                                          albedo(ig)=alb_n2b-deltab
     221                                         albedo(ig,:)=alb_n2b-deltab
    221222                         else if (((latitude(ig)*180./pi.lt.30.5).and.(latitude(ig)*180./pi.gt.29.).and. &
    222223                                   (longitude(ig)*180./pi.lt.165.5).and.(longitude(ig)*180./pi.gt.164.5)) &
    223224                             .or. ((latitude(ig)*180./pi.lt.33.).and.(latitude(ig)*180./pi.gt.32.).and.  &
    224225                                   (longitude(ig)*180./pi.lt.169.).and.(longitude(ig)*180./pi.gt.168.))) then
    225                                          albedo(ig)=alb_n2b+deltab
     226                                         albedo(ig,:)=alb_n2b+deltab
    226227                         else
    227                              albedo(ig)=alb_n2b
     228                             albedo(ig,:)=alb_n2b
    228229                         endif
    229230                     else
    230                              albedo(ig)=alb_n2b
     231                             albedo(ig,:)=alb_n2b
    231232                     endif
    232233
    233234                 CASE (7) ! Albedo from albedodat and adjusted emissivity
    234                    albedo(ig)=albedodat(ig)
     235                   albedo(ig,:)=albedodat(ig)
    235236                   ! adjust emissivity if convergeps = true
    236237                   emis(ig)=min(max(emis(ig)+adjust,0.),1.)
    … …  
    247248
    248249        else if (ice_case.eq.2) then
    249                albedo(ig)=alb_co
     250               albedo(ig,:)=alb_co
    250251               emis(ig)=emis_co
    251252
    … …  
    261262                   emis(ig)=emis_ch4
    262263                   if (latitude(ig)*180./pi.le.metlateq.and.latitude(ig)*180./pi.gt.-metlateq) then
    263                       albedo(ig)=alb_ch4_eq
     264                      albedo(ig,:)=alb_ch4_eq
    264265                   else
    265                       albedo(ig)=alb_ch4
     266                      albedo(ig,:)=alb_ch4
    266267                   endif
    267268
    … …  
    269270                   emis(ig)=emis_ch4
    270271                   if (latitude(ig)*180./pi.le.metlateq.and.latitude(ig)*180./pi.gt.-metlateq) then
    271                       albedo(ig)=alb_ch4_eq
     272                      albedo(ig,:)=alb_ch4_eq
    272273                   else if (latitude(ig)*180./pi.le.-metlateq) then
    273                       albedo(ig)=alb_ch4_s
     274                      albedo(ig,:)=alb_ch4_s
    274275                   else
    275                       albedo(ig)=alb_ch4
     276                      albedo(ig,:)=alb_ch4
    276277                   endif
    277278
    … …  
    279280
    280281                   emis(ig)=emis_ch4
    281                    albedo(ig)=alb_ch4
     282                   albedo(ig,:)=alb_ch4
    282283
    283284                   if (latitude(ig)*180./pi.le.metlateq.and.latitude(ig)*180./pi.gt.-metlateq) then
    284                       albedo(ig)=alb_ch4_eq
     285                      albedo(ig,:)=alb_ch4_eq
    285286                   endif
    286287
    … …  
    291292                          SELECT CASE (feedback_met)
    292293                            CASE (0) ! Default (linear from alb_ch4_eq)
    293                               albedo(ig)=(1.-alb_ch4_eq)/0.002*qsurf(ig,igcm_ch4_ice)+alb_ch4_eq
     294                              albedo(ig,:)=(1.-alb_ch4_eq)/0.002*qsurf(ig,igcm_ch4_ice)+alb_ch4_eq
    294295                              !emis(ig)=(1.-emis_ch4)/0.002*qsurf(ig,igcm_ch4_ice)+emis_ch4
    295                               if (albedo(ig).gt.fdch4_maxalb) albedo(ig)=fdch4_maxalb
     296                              if (albedo(ig,1).gt.fdch4_maxalb) albedo(ig,:)=fdch4_maxalb
    296297                              !if (emis(ig).gt.1.) emis(ig)=1.
    297298                            CASE (1) ! Hyperbolic tangent old
    298                               albedo(ig)=0.6*0.5*(1.+tanh(6.*(0.001+qsurf(ig,igcm_ch4_ice))/0.5))+0.3 !
     299                              albedo(ig,:)=0.6*0.5*(1.+tanh(6.*(0.001+qsurf(ig,igcm_ch4_ice))/0.5))+0.3 !
    299300                              !emis(ig)=0.2*0.5*(1.+tanh(6.*(0.001+qsurf(ig,igcm_ch4_ice))/0.5))+0.8 !
    300                               if (albedo(ig).gt.fdch4_maxalb) albedo(ig)=fdch4_maxalb
     301                              if (albedo(ig,1).gt.fdch4_maxalb) albedo(ig,:)=fdch4_maxalb
    301302                              !if (emis(ig).gt.1.) emis(ig)=1.
    302303                            CASE (2) ! hyperbolic tangent old
    303                               albedo(ig)=0.5*0.6*(1.+tanh(1000.*(qsurf(ig,igcm_ch4_ice)-0.002)))+0.3 !
     304                              albedo(ig,:)=0.5*0.6*(1.+tanh(1000.*(qsurf(ig,igcm_ch4_ice)-0.002)))+0.3 !
    304305                              !emis(ig)=0.2*0.5*(1.+tanh(1000.*(qsurf(ig,igcm_ch4_ice)-0.002)))+0.8 !
    305                               if (albedo(ig).gt.fdch4_maxalb) albedo(ig)=fdch4_maxalb
     306                              if (albedo(ig,1).gt.fdch4_maxalb) albedo(ig,:)=fdch4_maxalb
    306307                              !if (emis(ig).gt.1.) emis(ig)=1.
    307308                            CASE (3) ! hyperbolic tangent equation with parameters
    308                               albedo(ig)=aa*(-1+tanh(fdch4_ampl*qsurf(ig,igcm_ch4_ice)))+bb
    309                               if (albedo(ig).gt.fdch4_maxalb) albedo(ig)=fdch4_maxalb
     309                              albedo(ig,:)=aa*(-1+tanh(fdch4_ampl*qsurf(ig,igcm_ch4_ice)))+bb
     310                              if (albedo(ig,1).gt.fdch4_maxalb) albedo(ig,:)=fdch4_maxalb
    310311                            CASE DEFAULT
    311312                              write(*,*) 'STOP surfprop.F90: fd wrong'
    … …  
    319320                   emis(ig)=emis_ch4
    320321                   if (latitude(ig)*180./pi.le.metlateq.and.latitude(ig)*180./pi.gt.-metlateq) then
    321                       albedo(ig)=alb_ch4_eq     ! Pure methane ice
     322                      albedo(ig,:)=alb_ch4_eq     ! Pure methane ice
    322323                   else if (latitude(ig)*180./pi.le.-metlateq) then
    323                       albedo(ig)=alb_ch4_s     ! Pure methane ice
     324                      albedo(ig,:)=alb_ch4_s     ! Pure methane ice
    324325                      if (pls.le.metls2.and.pls.gt.metls1) then
    325                         albedo(ig)=alb_ch4_s+deltab  ! Also used for N2, careful
     326                        albedo(ig,:)=alb_ch4_s+deltab  ! Also used for N2, careful
    326327                      endif
    327328                   else
    328                       albedo(ig)=alb_ch4     ! Pure methane ice
     329                      albedo(ig,:)=alb_ch4     ! Pure methane ice
    329330                   endif
    330331
    331332                 CASE (4) ! Albedo from albedodat
    332333                   emis(ig)=emis_ch4
    333                    albedo(ig)=albedodat(ig)
     334                   albedo(ig,:)=albedodat(ig)
    334335
    335336                 CASE DEFAULT
    … …  
    349350
    350351                   if (latitude(ig)*180./pi.le.tholateq.and.latitude(ig)*180./pi.gt.-tholateq) then
    351                       albedo(ig)=alb_tho_eq
     352                      albedo(ig,:)=alb_tho_eq
    352353                      emis(ig)=emis_tho_eq
    353354                   else
    354                       albedo(ig)=alb_tho      ! default = 0.1
     355                      albedo(ig,:)=alb_tho      ! default = 0.1
    355356                      emis(ig)=emis_tho
    356357                   endif
    … …  
    359360
    360361                   if (latitude(ig)*180./pi.le.tholateq.and.latitude(ig)*180./pi.gt.-tholateq) then
    361                       albedo(ig)=alb_tho_eq
     362                      albedo(ig,:)=alb_tho_eq
    362363                      emis(ig)=emis_tho_eq
    363364                   else
    364                       albedo(ig)=alb_tho      ! default = 0.1
     365                      albedo(ig,:)=alb_tho      ! default = 0.1
    365366                      emis(ig)=emis_tho
    366367                   endif
    … …  
    368369                   if (latitude(ig)*180./pi.le.tholatn.and.latitude(ig)*180./pi.ge.tholats &
    369370                  .and.longitude(ig)*180./pi.ge.tholone.and.longitude(ig)*180./pi.le.tholonw) then
    370                       albedo(ig)=alb_tho_spe
     371                      albedo(ig,:)=alb_tho_spe
    371372                      emis(ig)=emis_tho_spe
    372373                   endif
    … …  
    374375                 CASE (2) ! Depends on the albedo map
    375376
    376                    albedo(ig)=albedodat(ig)
    377                    if (albedo(ig).gt.alb_ch4) then
     377                   albedo(ig,:)=albedodat(ig)
     378                   if (albedo(ig,1).gt.alb_ch4) then
    378379                      emis(ig)=emis_ch4
    379380                   else
    … …  
    391392
    392393                if (specalb) then
    393                   albedo(ig)=albedodat(ig)        ! Specific ground properties
     394                  albedo(ig,:)=albedodat(ig)        ! Specific ground properties
    394395                  !if (albedodat(ig).lt.0.25) then
    395                   !    albedo(ig)=alb_tho_eq
     396                  !    albedo(ig,:)=alb_tho_eq
    396397                  !else if (albedodat(ig).lt.0.40) then
    397                   !    albedo(ig)=0.35
     398                  !    albedo(ig,:)=0.35
    398399                  !else if (albedodat(ig).lt.0.70) then
    399                   !    albedo(ig)=0.51
     400                  !    albedo(ig,:)=0.51
    400401                  !endif
    401402                endif
Note: See TracChangeset for help on using the changeset viewer.