Changeset 2618 for LMDZ5/trunk


Ignore:
Timestamp:
Sep 1, 2016, 10:47:39 AM (8 years ago)
Author:
oboucher
Message:

I have moved the computation of cloud optical properties
into the radiation if .. .endif sequence, after the aerosol
optical properties. This follows up on rev 2614 and should
fix the 1+1=2 problem that was introduced back then.

File:
1 edited

Legend:

Unmodified
Added
Removed

  • LMDZ5/trunk/libf/phylmd/physiq_mod.F90

    r2614 r2618  
    33823382#endif
    33833383    END IF !type_trac = inca
    3384     !     
    3385     ! Calculer les parametres optiques des nuages et quelques
    3386     ! parametres pour diagnostiques:
    3387     !
    3388 
    3389     IF (aerosol_couple.AND.config_inca=='aero') THEN
    3390        mass_solu_aero(:,:)    = ccm(:,:,1)
    3391        mass_solu_aero_pi(:,:) = ccm(:,:,2)
    3392     END IF
    3393 
    3394     if (ok_newmicro) then
    3395        IF (iflag_rrtm.NE.0) THEN
     3384
     3385
     3386    !
     3387    ! Appeler le rayonnement mais calculer tout d'abord l'albedo du sol.
     3388    !
     3389    IF (MOD(itaprad,radpas).EQ.0) THEN
     3390
     3391       !
     3392       !jq - introduce the aerosol direct and first indirect radiative forcings
     3393       !jq - Johannes Quaas, 27/11/2003 (quaas@lmd.jussieu.fr)
     3394       IF (flag_aerosol .gt. 0) THEN
     3395          IF (iflag_rrtm .EQ. 0) THEN !--old radiation
     3396             IF (.NOT. aerosol_couple) THEN
     3397                !
     3398                CALL readaerosol_optic( &
     3399                     debut, new_aod, flag_aerosol, itap, jD_cur-jD_ref, &
     3400                     pdtphys, pplay, paprs, t_seri, rhcl, presnivs,  &
     3401                     mass_solu_aero, mass_solu_aero_pi,  &
     3402                     tau_aero, piz_aero, cg_aero,  &
     3403                     tausum_aero, tau3d_aero)
     3404             ENDIF
     3405          ELSE                       ! RRTM radiation
     3406             IF (aerosol_couple .AND. config_inca == 'aero' ) THEN
     3407                abort_message='config_inca=aero et rrtm=1 impossible'
     3408                call abort_physic(modname,abort_message,1)
     3409             ELSE
     3410                !
    33963411#ifdef CPP_RRTM
    3397           IF (ok_cdnc.AND.NRADLP.NE.3) THEN
     3412                IF (NSW.EQ.6) THEN
     3413                   !--new aerosol properties
     3414                   !
     3415                   CALL readaerosol_optic_rrtm( debut, aerosol_couple, &
     3416                        new_aod, flag_aerosol, itap, jD_cur-jD_ref, &
     3417                        pdtphys, pplay, paprs, t_seri, rhcl, presnivs,  &
     3418                        tr_seri, mass_solu_aero, mass_solu_aero_pi,  &
     3419                        tau_aero_sw_rrtm, piz_aero_sw_rrtm, cg_aero_sw_rrtm,  &
     3420                        tausum_aero, tau3d_aero)
     3421
     3422                ELSE IF (NSW.EQ.2) THEN
     3423                   !--for now we use the old aerosol properties
     3424                   !
     3425                   CALL readaerosol_optic( &
     3426                        debut, new_aod, flag_aerosol, itap, jD_cur-jD_ref, &
     3427                        pdtphys, pplay, paprs, t_seri, rhcl, presnivs,  &
     3428                        mass_solu_aero, mass_solu_aero_pi,  &
     3429                        tau_aero, piz_aero, cg_aero,  &
     3430                        tausum_aero, tau3d_aero)
     3431                   !
     3432                   !--natural aerosols
     3433                   tau_aero_sw_rrtm(:,:,1,:)=tau_aero(:,:,3,:)
     3434                   piz_aero_sw_rrtm(:,:,1,:)=piz_aero(:,:,3,:)
     3435                   cg_aero_sw_rrtm (:,:,1,:)=cg_aero (:,:,3,:)
     3436                   !--all aerosols
     3437                   tau_aero_sw_rrtm(:,:,2,:)=tau_aero(:,:,2,:)
     3438                   piz_aero_sw_rrtm(:,:,2,:)=piz_aero(:,:,2,:)
     3439                   cg_aero_sw_rrtm (:,:,2,:)=cg_aero (:,:,2,:)
     3440                ELSE
     3441                   abort_message='Only NSW=2 or 6 are possible with ' &
     3442                        // 'aerosols and iflag_rrtm=1'
     3443                   call abort_physic(modname,abort_message,1)
     3444                ENDIF
     3445
     3446                !--call LW optical properties for tropospheric aerosols
     3447                !--only works for INCA aerosol (aerosol_couple = TRUE)
     3448                CALL aeropt_lw_rrtm(aerosol_couple,paprs,tr_seri)
     3449                !
     3450#else
     3451                abort_message='You should compile with -rrtm if running ' &
     3452                     // 'with iflag_rrtm=1'
     3453                call abort_physic(modname,abort_message,1)
     3454#endif
     3455                !
     3456             ENDIF
     3457          ENDIF
     3458       ELSE
     3459          tausum_aero(:,:,:) = 0.
     3460          IF (iflag_rrtm .EQ. 0) THEN !--old radiation
     3461             tau_aero(:,:,:,:) = 1.e-15
     3462             piz_aero(:,:,:,:) = 1.
     3463             cg_aero(:,:,:,:)  = 0.
     3464          ELSE
     3465             tau_aero_sw_rrtm(:,:,:,:) = 1.e-15
     3466             tau_aero_lw_rrtm(:,:,:,:) = 1.e-15
     3467             piz_aero_sw_rrtm(:,:,:,:) = 1.0
     3468             cg_aero_sw_rrtm(:,:,:,:)  = 0.0
     3469          ENDIF
     3470       ENDIF
     3471       !
     3472       !--STRAT AEROSOL
     3473       !--updates tausum_aero,tau_aero,piz_aero,cg_aero
     3474       IF (flag_aerosol_strat.GT.0) THEN
     3475          IF (prt_level .GE.10) THEN
     3476             PRINT *,'appel a readaerosolstrat', mth_cur
     3477          ENDIF
     3478          IF (iflag_rrtm.EQ.0) THEN
     3479           IF (flag_aerosol_strat.EQ.1) THEN
     3480             CALL readaerosolstrato(debut)
     3481           ELSE
     3482             abort_message='flag_aerosol_strat must equal 1 for rrtm=0'
     3483             CALL abort_physic(modname,abort_message,1)
     3484           ENDIF
     3485          ELSE
     3486#ifdef CPP_RRTM
     3487            IF (flag_aerosol_strat.EQ.1) THEN
     3488             CALL readaerosolstrato1_rrtm(debut)
     3489            ELSEIF (flag_aerosol_strat.EQ.2) THEN
     3490             CALL stratosphere_mask(t_seri, pplay, latitude_deg)
     3491             CALL readaerosolstrato2_rrtm(debut)
     3492            ELSE
     3493             abort_message='flag_aerosol_strat must equal 1 or 2 for rrtm=1'
     3494             CALL abort_physic(modname,abort_message,1)
     3495            ENDIF
     3496#else
     3497             abort_message='You should compile with -rrtm if running ' &
     3498                  // 'with iflag_rrtm=1'
     3499             CALL abort_physic(modname,abort_message,1)
     3500#endif
     3501          ENDIF
     3502       ENDIF
     3503       !--fin STRAT AEROSOL
     3504       !     
     3505
     3506       ! Calculer les parametres optiques des nuages et quelques
     3507       ! parametres pour diagnostiques:
     3508       !
     3509       IF (aerosol_couple.AND.config_inca=='aero') THEN
     3510          mass_solu_aero(:,:)    = ccm(:,:,1)
     3511          mass_solu_aero_pi(:,:) = ccm(:,:,2)
     3512       END IF
     3513
     3514       IF (ok_newmicro) then
     3515          IF (iflag_rrtm.NE.0) THEN
     3516#ifdef CPP_RRTM
     3517             IF (ok_cdnc.AND.NRADLP.NE.3) THEN
    33983518             abort_message='RRTM choix incoherent NRADLP doit etre egal a 3 ' &
    33993519                  // 'pour ok_cdnc'
    3400              call abort_physic(modname,abort_message,1)
    3401           endif
     3520             CALL abort_physic(modname,abort_message,1)
     3521             ENDIF
    34023522#else
    34033523
    3404           abort_message='You should compile with -rrtm if running with ' &
    3405                // 'iflag_rrtm=1'
    3406           call abort_physic(modname,abort_message,1)
     3524             abort_message='You should compile with -rrtm if running with '//'iflag_rrtm=1'
     3525             CALL abort_physic(modname,abort_message,1)
    34073526#endif
     3527          ENDIF
     3528          CALL newmicro (ok_cdnc, bl95_b0, bl95_b1, &
     3529               paprs, pplay, t_seri, cldliq, cldfra, &
     3530               cldtau, cldemi, cldh, cldl, cldm, cldt, cldq, &
     3531               flwp, fiwp, flwc, fiwc, &
     3532               mass_solu_aero, mass_solu_aero_pi, &
     3533               cldtaupi, re, fl, ref_liq, ref_ice, &
     3534               ref_liq_pi, ref_ice_pi)
     3535       ELSE
     3536          CALL nuage (paprs, pplay, &
     3537               t_seri, cldliq, cldfra, cldtau, cldemi, &
     3538               cldh, cldl, cldm, cldt, cldq, &
     3539               ok_aie, &
     3540               mass_solu_aero, mass_solu_aero_pi, &
     3541               bl95_b0, bl95_b1, &
     3542               cldtaupi, re, fl)
    34083543       ENDIF
    3409        CALL newmicro (ok_cdnc, bl95_b0, bl95_b1, &
    3410             paprs, pplay, t_seri, cldliq, cldfra, &
    3411             cldtau, cldemi, cldh, cldl, cldm, cldt, cldq, &
    3412             flwp, fiwp, flwc, fiwc, &
    3413             mass_solu_aero, mass_solu_aero_pi, &
    3414             cldtaupi, re, fl, ref_liq, ref_ice, &
    3415             ref_liq_pi, ref_ice_pi)
    3416     else
    3417        CALL nuage (paprs, pplay, &
    3418             t_seri, cldliq, cldfra, cldtau, cldemi, &
    3419             cldh, cldl, cldm, cldt, cldq, &
    3420             ok_aie, &
    3421             mass_solu_aero, mass_solu_aero_pi, &
    3422             bl95_b0, bl95_b1, &
    3423             cldtaupi, re, fl)
    3424     endif
    3425     !
    3426     !IM betaCRF
    3427     !
    3428     cldtaurad   = cldtau
    3429     cldtaupirad = cldtaupi
    3430     cldemirad   = cldemi
    3431     cldfrarad   = cldfra
    3432 
    3433     !
    3434     if(lon1_beta.EQ.-180..AND.lon2_beta.EQ.180..AND. &
    3435          lat1_beta.EQ.90..AND.lat2_beta.EQ.-90.) THEN
    3436        !
    3437        ! global
    3438        !
    3439        DO k=1, klev
    3440           DO i=1, klon
    3441              if (pplay(i,k).GE.pfree) THEN
    3442                 beta(i,k) = beta_pbl
    3443              else
    3444                 beta(i,k) = beta_free
    3445              endif
    3446              if (mskocean_beta) THEN
    3447                 beta(i,k) = beta(i,k) * pctsrf(i,is_oce)
    3448              endif
    3449              cldtaurad(i,k)   = cldtau(i,k) * beta(i,k)
    3450              cldtaupirad(i,k) = cldtaupi(i,k) * beta(i,k)
    3451              cldemirad(i,k)   = cldemi(i,k) * beta(i,k)
    3452              cldfrarad(i,k)   = cldfra(i,k) * beta(i,k)
    3453           ENDDO
    3454        ENDDO
    3455        !
    3456     else
    3457        !
    3458        ! regional
    3459        !
    3460        DO k=1, klev
    3461           DO i=1,klon
    3462              !
    3463              if (longitude_deg(i).ge.lon1_beta.AND. &
    3464                   longitude_deg(i).le.lon2_beta.AND. &
    3465                   latitude_deg(i).le.lat1_beta.AND. &
    3466                   latitude_deg(i).ge.lat2_beta) THEN
    3467                 if (pplay(i,k).GE.pfree) THEN
     3544       !
     3545       !IM betaCRF
     3546       !
     3547       cldtaurad   = cldtau
     3548       cldtaupirad = cldtaupi
     3549       cldemirad   = cldemi
     3550       cldfrarad   = cldfra
     3551
     3552       !
     3553       IF (lon1_beta.EQ.-180..AND.lon2_beta.EQ.180..AND. &
     3554           lat1_beta.EQ.90..AND.lat2_beta.EQ.-90.) THEN
     3555          !
     3556          ! global
     3557          !
     3558          DO k=1, klev
     3559             DO i=1, klon
     3560                IF (pplay(i,k).GE.pfree) THEN
    34683561                   beta(i,k) = beta_pbl
    3469                 else
     3562                ELSE
    34703563                   beta(i,k) = beta_free
    3471                 endif
    3472                 if (mskocean_beta) THEN
     3564                ENDIF
     3565                IF (mskocean_beta) THEN
    34733566                   beta(i,k) = beta(i,k) * pctsrf(i,is_oce)
    3474                 endif
     3567                ENDIF
    34753568                cldtaurad(i,k)   = cldtau(i,k) * beta(i,k)
    34763569                cldtaupirad(i,k) = cldtaupi(i,k) * beta(i,k)
    34773570                cldemirad(i,k)   = cldemi(i,k) * beta(i,k)
    34783571                cldfrarad(i,k)   = cldfra(i,k) * beta(i,k)
    3479              endif
     3572             ENDDO
     3573          ENDDO
     3574          !
     3575       ELSE
     3576          !
     3577          ! regional
     3578          !
     3579          DO k=1, klev
     3580             DO i=1,klon
     3581                !
     3582                IF (longitude_deg(i).ge.lon1_beta.AND. &
     3583                    longitude_deg(i).le.lon2_beta.AND. &
     3584                    latitude_deg(i).le.lat1_beta.AND.  &
     3585                    latitude_deg(i).ge.lat2_beta) THEN
     3586                   IF (pplay(i,k).GE.pfree) THEN
     3587                      beta(i,k) = beta_pbl
     3588                   ELSE
     3589                      beta(i,k) = beta_free
     3590                   ENDIF
     3591                   IF (mskocean_beta) THEN
     3592                      beta(i,k) = beta(i,k) * pctsrf(i,is_oce)
     3593                   ENDIF
     3594                   cldtaurad(i,k)   = cldtau(i,k) * beta(i,k)
     3595                   cldtaupirad(i,k) = cldtaupi(i,k) * beta(i,k)
     3596                   cldemirad(i,k)   = cldemi(i,k) * beta(i,k)
     3597                   cldfrarad(i,k)   = cldfra(i,k) * beta(i,k)
     3598                ENDIF
    34803599             !
     3600             ENDDO
    34813601          ENDDO
    3482        ENDDO
    3483        !
    3484     endif
    3485     !
    3486     ! Appeler le rayonnement mais calculer tout d'abord l'albedo du sol.
    3487     !
    3488     IF (MOD(itaprad,radpas).EQ.0) THEN
    3489 
    3490        !albedo SB >>> 
    3491        if(ok_chlorophyll)then
     3602       !
     3603       ENDIF
     3604
     3605       !lecture de la chlorophylle pour le nouvel albedo de Sunghye Baek
     3606       IF (ok_chlorophyll) THEN
    34923607          print*,"-- reading chlorophyll"
    3493           call readchlorophyll(debut)
    3494        endif
    3495        !do i=1,klon
    3496        !if(chl_con(i)>1.) print*,i,chl_con(i),pctsrf(i,is_ter)
    3497        !enddo
    3498        !albedo SB <<<
     3608          CALL readchlorophyll(debut)
     3609       ENDIF
    34993610
    35003611       !
    … …  
    36593770               cldtaupirad, &
    36603771               topswai_aero, solswai_aero)
    3661 
    36623772#endif
    36633773       ELSE
    … …  
    36743784             print *,' ->radlwsw, number 1 '
    36753785          ENDIF
    3676 
    36773786          !
    36783787          CALL radlwsw &
    … …  
    38143923    !
    38153924    radsol=solsw*swradcorr+sollw
     3925
    38163926    if (ok_4xCO2atm) then
    38173927       radsolp=solswp*swradcorr+sollwp
Note: See TracChangeset for help on using the changeset viewer.