Ignore:
Timestamp:
Jan 11, 2017, 3:33:51 PM (8 years ago)
Author:
jvatant
Message:

+ Major clean of the new LMDZ.TITAN from too-generic options and routines (water, co2, ocean, surface type ...)
+ From this revision LMDZ.TITAN begins to be really separated from LMDZ.GENERIC
+ Partial desactivation of aerosols, only the dummy case is still enabled to keep the code running ( new aerosol routines to come in followings commits )

JVO

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/LMDZ.TITAN/libf/phytitan/aeropacity.F90

    r1542 r1647  
    11      Subroutine aeropacity(ngrid,nlayer,nq,pplay,pplev,pq, &
    2          aerosol,reffrad,QREFvis3d,QREFir3d,tau_col,cloudfrac,totcloudfrac,clearsky)
     2         aerosol,reffrad,QREFvis3d,QREFir3d,tau_col)
    33
    44       use radinc_h, only : L_TAUMAX,naerkind
    55       use aerosol_mod
    6        USE tracer_h, only: noms,rho_co2,rho_ice
     6       USE tracer_h, only: noms
    77       use comcstfi_mod, only: g
    8        use callkeys_mod, only: aerofixco2,aerofixh2o,kastprof,cloudlvl, &
    9                 CLFvarying,CLFfixval,dusttau,                           &
    10                 pres_bottom_tropo,pres_top_tropo,obs_tau_col_tropo,     &
     8       use callkeys_mod, only: pres_bottom_tropo,pres_top_tropo,obs_tau_col_tropo,  &
    119                pres_bottom_strato,pres_top_strato,obs_tau_col_strato
    1210                 
    … …  
    5553      REAL,INTENT(IN) :: QREFir3d(ngrid,nlayer,naerkind)
    5654      REAL,INTENT(OUT):: tau_col(ngrid) !column integrated visible optical depth
    57       ! BENJAMIN MODIFS
    58       real,intent(in) :: cloudfrac(ngrid,nlayer) ! cloud fraction
    59       real,intent(out) :: totcloudfrac(ngrid) ! total cloud fraction
    60       logical,intent(in) :: clearsky
    6155
    6256      real aerosol0
    … …  
    6963      EXTERNAL CBRT
    7064
    71       INTEGER,SAVE :: i_co2ice=0      ! co2 ice
    72       INTEGER,SAVE :: i_h2oice=0      ! water ice
    73 !$OMP THREADPRIVATE(i_co2ice,i_h2oice)
    74       CHARACTER(LEN=20) :: tracername ! to temporarily store text
    75 
    7665      ! for fixed dust profiles
    77       real topdust, expfactor, zp
    78       REAL taudusttmp(ngrid) ! Temporary dust opacity used before scaling
    79       REAL tauh2so4tmp(ngrid) ! Temporary h2so4 opacity used before scaling
    80 
    81       real CLFtot
     66      real expfactor, zp
    8267
    8368      ! identify tracers
    … …  
    8570
    8671        write(*,*) "Tracers found in aeropacity:"
    87         do iq=1,nq
    88           tracername=noms(iq)
    89           if (tracername.eq."co2_ice") then
    90             i_co2ice=iq
    91           write(*,*) "i_co2ice=",i_co2ice
    92 
    93           endif
    94           if (tracername.eq."h2o_ice") then
    95             i_h2oice=iq
    96             write(*,*) "i_h2oice=",i_h2oice
    97           endif
    98         enddo
    9972
    10073        if (noaero) then
    … …  
    10780        endif
    10881
    109         if ((iaero_co2.ne.0).and.(.not.noaero)) then
    110           print*, 'iaero_co2=  ',iaero_co2
    111         endif
    112         if (iaero_h2o.ne.0) then
    113           print*,'iaero_h2o=  ',iaero_h2o   
    114         endif
    115         if (iaero_dust.ne.0) then
    116           print*,'iaero_dust= ',iaero_dust
    117         endif
    118         if (iaero_h2so4.ne.0) then
    119           print*,'iaero_h2so4= ',iaero_h2so4
    120         endif
    12182        if (iaero_back2lay.ne.0) then
    12283          print*,'iaero_back2lay= ',iaero_back2lay
    … …  
    12788
    12889
    129 !     ---------------------------------------------------------
    130 !==================================================================
    131 !    CO2 ice aerosols
    132 !==================================================================
     90!     ---------------------------------------------------------   
    13391
    134       if (iaero_co2.ne.0) then
    135            iaer=iaero_co2
    136 !       1. Initialization
    137             aerosol(1:ngrid,1:nlayer,iaer)=0.0
    138 !       2. Opacity calculation
    139             if (noaero) then ! aerosol set to zero
    140              aerosol(1:ngrid,1:nlayer,iaer)=0.0
    141             elseif (aerofixco2.or.(i_co2ice.eq.0)) then !  CO2 ice cloud prescribed
    142                aerosol(1:ngrid,1:nlayer,iaer)=1.e-9
    143                !aerosol(1:ngrid,12,iaer)=4.0 ! single cloud layer option
    144             else
    145                DO ig=1, ngrid
    146                   DO l=1,nlayer-1 ! to stop the rad tran bug
    147 
    148                      aerosol0 =                         &
    149                           (  0.75 * QREFvis3d(ig,l,iaer) /        &
    150                           ( rho_co2 * reffrad(ig,l,iaer) )  ) *   &
    151                           ( pq(ig,l,i_co2ice) + 1.E-9 ) *         &
    152                           ( pplev(ig,l) - pplev(ig,l+1) ) / g
    153                      aerosol0           = max(aerosol0,1.e-9)
    154                      aerosol0           = min(aerosol0,L_TAUMAX)
    155                      aerosol(ig,l,iaer) = aerosol0
    156 !                     aerosol(ig,l,iaer) = 0.0
    157 !                     print*, aerosol(ig,l,iaer)
    158 !        using cloud fraction
    159 !                     aerosol(ig,l,iaer) = -log(1 - CLF + CLF*exp(-aerosol0/CLF))
    160 !                     aerosol(ig,l,iaer) = min(aerosol(ig,l,iaer),L_TAUMAX)
    161 
    162 
    163                   ENDDO
    164                ENDDO
    165             end if ! if fixed or varying
    166       end if ! if CO2 aerosols   
    167 !==================================================================
    168 !     Water ice / liquid
    169 !==================================================================
    170 
    171       if (iaero_h2o.ne.0) then
    172            iaer=iaero_h2o
    173 !       1. Initialization
    174             aerosol(1:ngrid,1:nlayer,iaer)=0.0
    175 !       2. Opacity calculation
    176             if (aerofixh2o.or.(i_h2oice.eq.0).or.clearsky) then
    177                aerosol(1:ngrid,1:nlayer,iaer) =1.e-9
    178 
    179                ! put cloud at cloudlvl
    180                if(kastprof.and.(cloudlvl.ne.0.0))then
    181                   ig=1
    182                   do l=1,nlayer
    183                      if(int(cloudlvl).eq.l)then
    184                      !if(cloudlvl.gt.(pplay(ig,l)/pplev(ig,1)))then
    185                         print*,'Inserting cloud at level ',l
    186                         !aerosol(ig,l,iaer)=10.0
    187 
    188                         rho_ice=920.0
    189 
    190                         ! the Kasting approximation
    191                         aerosol(ig,l,iaer) =                      &
    192                           (  0.75 * QREFvis3d(ig,l,iaer) /        &
    193                           ( rho_ice * reffrad(ig,l,iaer) )  ) *   &
    194                           !( pq(ig,l,i_h2oice) + 1.E-9 ) *         &
    195                           ( 4.0e-4 + 1.E-9 ) *         &
    196                           ( pplev(ig,l) - pplev(ig,l+1) ) / g
    197 
    198 
    199                         open(115,file='clouds.out',form='formatted')
    200                         write(115,*) l,aerosol(ig,l,iaer)
    201                         close(115)
    202 
    203                         return
    204                      endif
    205                   end do
    206 
    207                   call abort
    208                endif
    209 
    210             else
    211 
    212                do ig=1, ngrid
    213                   do l=1,nlayer-1 ! to stop the rad tran bug
    214 
    215                      aerosol(ig,l,iaer) =                                    & !modification by BC
    216                           (  0.75 * QREFvis3d(ig,l,iaer) /        &
    217                           ( rho_ice * reffrad(ig,l,iaer) )  ) *   &
    218                           !  pq(ig,l,i_h2oice) *                   & !JL I dropped the +1e-9 here to have the same
    219                           !( pplev(ig,l) - pplev(ig,l+1) ) / g       !   opacity in the clearsky=true and the
    220                                                                      !   clear=false/pq=0 case
    221                           ( pq(ig,l,i_h2oice) + 1.E-9 ) *         & ! Doing this makes the code unstable, so I have restored it (RW)
    222                           ( pplev(ig,l) - pplev(ig,l+1) ) / g
    223 
    224                   enddo
    225                enddo
    226 
    227                if(CLFvarying)then
    228                   call totalcloudfrac(ngrid,nlayer,nq,cloudfrac,totcloudfrac,pplev,pq,aerosol(1,1,iaer))
    229                   do ig=1, ngrid
    230                      do l=1,nlayer-1 ! to stop the rad tran bug
    231                         CLFtot  = max(totcloudfrac(ig),0.01)
    232                         aerosol(ig,l,iaer)=aerosol(ig,l,iaer)/CLFtot
    233                         aerosol(ig,l,iaer) = max(aerosol(ig,l,iaer),1.e-9)
    234                      enddo
    235                   enddo
    236                else
    237                   do ig=1, ngrid
    238                      do l=1,nlayer-1 ! to stop the rad tran bug
    239                         CLFtot  = CLFfixval
    240                         aerosol(ig,l,iaer)=aerosol(ig,l,iaer)/CLFtot
    241                         aerosol(ig,l,iaer) = max(aerosol(ig,l,iaer),1.e-9)
    242                      enddo
    243                   enddo
    244               end if!(CLFvarying)
    245             endif !(aerofixed.or.(i_h2oice.eq.0).or.clearsky)
    246              
    247       end if ! End if h2o aerosol
    248 
    249 !==================================================================
    250 !             Dust
    251 !==================================================================
    252       if (iaero_dust.ne.0) then
    253           iaer=iaero_dust
    254 !         1. Initialization
    255           aerosol(1:ngrid,1:nlayer,iaer)=0.0
    256          
    257           topdust=30.0 ! km  (used to be 10.0 km) LK
    258 
    259 !       2. Opacity calculation
    260 
    261 !           expfactor=0.
    262            DO l=1,nlayer-1
    263              DO ig=1,ngrid
    264 !             Typical mixing ratio profile
    265 
    266                  zp=(pplev(ig,1)/pplay(ig,l))**(70./topdust)
    267                  expfactor=max(exp(0.007*(1.-max(zp,1.))),1.e-3)
    268 
    269 !             Vertical scaling function
    270               aerosol(ig,l,iaer)= (pplev(ig,l)-pplev(ig,l+1)) &
    271                *expfactor
    272 
    273 
    274              ENDDO
    275            ENDDO
    276 
    277 !          Rescaling each layer to reproduce the choosen (or assimilated)
    278 !          dust extinction opacity at visible reference wavelength, which
    279 !          is scaled to the surface pressure pplev(ig,1)
    280 
    281             taudusttmp(1:ngrid)=0.
    282               DO l=1,nlayer
    283                 DO ig=1,ngrid
    284                    taudusttmp(ig) = taudusttmp(ig) &
    285                           +  aerosol(ig,l,iaer)
    286                 ENDDO
    287               ENDDO
    288             DO l=1,nlayer-1
    289                DO ig=1,ngrid
    290                   aerosol(ig,l,iaer) = max(1E-20, &
    291                           dusttau &
    292                        *  pplev(ig,1) / pplev(ig,1) &
    293                        *  aerosol(ig,l,iaer) &
    294                        /  taudusttmp(ig))
    295 
    296               ENDDO
    297             ENDDO
    298       end if ! If dust aerosol   
    299 
    300 !==================================================================
    301 !           H2SO4
    302 !==================================================================
    303 ! added by LK
    304       if (iaero_h2so4.ne.0) then
    305          iaer=iaero_h2so4
    306 
    307 !       1. Initialization
    308          aerosol(1:ngrid,1:nlayer,iaer)=0.0
    309 
    310 
    311 !       2. Opacity calculation
    312 
    313 !           expfactor=0.
    314          DO l=1,nlayer-1
    315             DO ig=1,ngrid
    316 !              Typical mixing ratio profile
    317 
    318                zp=(pplev(ig,1)/pplay(ig,l))**(70./30) !emulating topdust
    319                expfactor=max(exp(0.007*(1.-max(zp,1.))),1.e-3)
    320 
    321 !             Vertical scaling function
    322                aerosol(ig,l,iaer)= (pplev(ig,l)-pplev(ig,l+1))*expfactor
    323 
    324             ENDDO
    325          ENDDO
    326          tauh2so4tmp(1:ngrid)=0.
    327          DO l=1,nlayer
    328             DO ig=1,ngrid
    329                tauh2so4tmp(ig) = tauh2so4tmp(ig) + aerosol(ig,l,iaer)
    330             ENDDO
    331          ENDDO
    332          DO l=1,nlayer-1
    333             DO ig=1,ngrid
    334                aerosol(ig,l,iaer) = max(1E-20, &
    335                           1 &
    336                        *  pplev(ig,1) / pplev(ig,1) &
    337                        *  aerosol(ig,l,iaer) &
    338                        /  tauh2so4tmp(ig))
    339 
    340             ENDDO
    341          ENDDO
    342 
    343 ! 1/700. is assuming a "sulfurtau" of 1
    344 ! Sulfur aerosol routine to be improved.
    345 !                     aerosol0 =                         &
    346 !                          (  0.75 * QREFvis3d(ig,l,iaer) /        &
    347 !                          ( rho_h2so4 * reffrad(ig,l,iaer) )  ) *   &
    348 !                          ( pq(ig,l,i_h2so4) + 1.E-9 ) *         &
    349 !                          ( pplev(ig,l) - pplev(ig,l+1) ) / g
    350 !                     aerosol0           = max(aerosol0,1.e-9)
    351 !                     aerosol0           = min(aerosol0,L_TAUMAX)
    352 !                     aerosol(ig,l,iaer) = aerosol0
    353 
    354 !                  ENDDO
    355 !               ENDDO
    356       end if
    357  
    358            
    359 !     ---------------------------------------------------------
    36092!==================================================================
    36193!    Two-layer aerosols (unknown composition)
Note: See TracChangeset for help on using the changeset viewer.