Changeset 5888 for LMDZ6/branches


Ignore:
Timestamp:
Nov 25, 2025, 9:59:23 AM (7 weeks ago)
Author:
evignon
Message:

continue the work on water isotopes implementation in blowing snow.

Location:
LMDZ6/branches/lmdz-snow/libf
Files:
2 added
8 edited

Legend:

Unmodified
Added
Removed

  • LMDZ6/branches/lmdz-snow/libf/phylmd/lmdz_blowing_snow_sublim_sedim.f90

    r5872 r5888  
    66!==============================================================================
    77! Routine that calculates the evaporation and sedimentation of blowing snow
    8 ! inspired by what is done in lscp_mod
     8! inspired by what is done in the cloud routines lscp
     9! Details are provided in Vignon et al. 2025, GMD
     10!
    911! Etienne Vignon, October 2022
    1012!==============================================================================

  • LMDZ6/branches/lmdz-snow/libf/phylmd/phys_state_var_mod.F90

    r5776 r5888  
    111111      REAL,ALLOCATABLE,SAVE :: xtsol(:,:),fxtevap(:,:,:)
    112112!$OMP THREADPRIVATE(xtsol,fxtevap)
    113       REAL, ALLOCATABLE, SAVE :: xtrain_fall(:,:), xtsnow_fall(:,:)
    114 !$OMP THREADPRIVATE(xtrain_fall,xtsnow_fall)
     113      REAL, ALLOCATABLE, SAVE :: xtrain_fall(:,:), xtsnow_fall(:,:),xtbs_fall(:,:)
     114!$OMP THREADPRIVATE(xtrain_fall,xtsnow_fall,xtbs_fall)
    115115#endif
    116116
    … …  
    872872      ALLOCATE(xtrain_fall(ntraciso,klon))
    873873      ALLOCATE(xtsnow_fall(ntraciso,klon))
     874      ALLOCATE(xtbs_fall(ntraciso,klon))
    874875      ALLOCATE(xtrain_con(ntraciso,klon))
    875876      ALLOCATE(xtsnow_con(ntraciso,klon))
    … …  
    10921093      DEALLOCATE(xtsol,fxtevap) 
    10931094      DEALLOCATE(xt_ancien,xtl_ancien,xts_ancien,xtbs_ancien, fxtd, wake_deltaxt)
    1094       DEALLOCATE(xtrain_fall, xtsnow_fall, xtrain_con, xtsnow_con)
     1095      DEALLOCATE(xtrain_fall, xtsnow_fall, xtbs_fall, xtrain_con, xtsnow_con)
    10951096#ifdef ISOTRAC
    10961097      DEALLOCATE(bassin_map,boite_map)

  • LMDZ6/branches/lmdz-snow/libf/phylmdiso/add_phys_tend_mod.F90

    r5285 r5888  
    1818SUBROUTINE add_pbl_tend(zdu, zdv, zdt, zdq, zdql, zdqi, zdqbs, paprs, text,abortphy,flag_inhib_tend, itap &
    1919#ifdef ISO
    20         ,zdxt,zdxtl,zdxti &
     20        ,zdxt,zdxtl,zdxti,zdxtbs &
    2121#endif
    2222        )
    … …  
    6868  INTEGER itap ! time step number
    6969#ifdef ISO
    70   REAL zdxt(ntraciso,klon, klev), zdxtl(ntraciso,klon, klev), zdxti(ntraciso,klon, klev)
     70  REAL zdxt(ntraciso,klon, klev), zdxtl(ntraciso,klon, klev), zdxti(ntraciso,klon, klev), zdxtbs(ntraciso,klon, klev)
    7171#endif
    7272
    … …  
    106106    CALL add_phys_tend(zdu, zdv, zzdt, zzdq, zdql, zdqi, zdqbs, paprs, text,abortphy,flag_inhib_tend, itap, 0 &
    107107#ifdef ISO
    108         ,zzdxt,zdxtl,zdxti &
     108        ,zzdxt,zdxtl,zdxti,zdxtbs &
    109109#endif
    110110        )
    … …  
    112112    CALL add_phys_tend(zdu, zdv, zdt, zdq, zdql, zdqi, zdqbs, paprs, text,abortphy,flag_inhib_tend, itap, 0 &
    113113#ifdef ISO
    114         ,zdxt,zdxtl,zdxti &
     114        ,zdxt,zdxtl,zdxti,zdxtbs &
    115115#endif
    116116        )
    … …  
    126126                          abortphy,flag_inhib_tend, itap, diag_mode &
    127127#ifdef ISO
    128         ,zdxt,zdxtl,zdxti &
     128        ,zdxt,zdxtl,zdxti,zdxtbs &
    129129#endif
    130130        )
    … …  
    180180#ifdef ISO
    181181REAL zdxt(ntraciso,klon,klev),zdxtl(ntraciso,klon,klev) &
    182         ,zdxti(ntraciso,klon,klev)
     182        ,zdxti(ntraciso,klon,klev), zdxtbs(ntraciso,klon,klev)
    183183integer ixt
    184184#endif
    … …  
    203203#ifdef ISO
    204204REAL, DIMENSION(ntraciso,klon,klev)   :: sav_zdxt
    205 REAL, DIMENSION(ntraciso,klon,klev)   :: sav_xtl_seri, sav_xts_seri, sav_xt_seri
     205REAL, DIMENSION(ntraciso,klon,klev)   :: sav_xtl_seri, sav_xts_seri, sav_xtbs_seri, sav_xt_seri
    206206#endif
    207207!
    … …  
    346346     xtl_seri(:,:,:)=xtl_seri(:,:,:)+zdxtl(:,:,:)
    347347     xts_seri(:,:,:)=xts_seri(:,:,:)+zdxti(:,:,:) 
     348     xtbs_seri(:,:,:)=xtbs_seri(:,:,:)+zdxtbs(:,:,:)
    348349#endif
    349350
    … …  
    661662    xtl_seri(:,:,:) = sav_xtl_seri(:,:,:)
    662663    xts_seri(:,:,:) = sav_xts_seri(:,:,:)
     664    xtbs_seri(:,:,:) = sav_xtbs_seri(:,:,:)
    663665    xt_seri(:,:,:)  = sav_xt_seri(:,:,:)
    664666    zdxt(:,:,:)     = sav_zdxt(:,:,:)

  • LMDZ6/branches/lmdz-snow/libf/phylmdiso/isotopes_routines_mod.F90

    r5748 r5888  
    1890118901subroutine dispatch(klon,klev,qx,q_seri,xt_seri,ql_seri,xtl_seri,qs_seri,xts_seri)
    1890218902
    18903 use infotrac_phy, ONLY: nqtot,nqo,ivap,iliq,isol,iqIsoPha,ntraciso=>ntiso
     18903use infotrac_phy, ONLY: nqtot,nqo,ivap,iliq,isol,ibs,iqIsoPha,ntraciso=>ntiso
    1890418904implicit none
    1890518905
    … …  
    1890718907integer, intent(in) :: klon,klev
    1890818908real,dimension(klon,klev,nqtot), intent(in) ::qx
     18909logical :: ok_bs
    1890918910
    1891018911! outputs

  • LMDZ6/branches/lmdz-snow/libf/phylmdiso/phyaqua_mod.F90

    r5662 r5888  
    388388    ql_ancien   = 0.
    389389    qs_ancien   = 0.
     390    qbs_ancien  = 0.
    390391    u_ancien    = 0.
    391392    v_ancien    = 0.
    … …  
    393394    prlw_ancien = 0.
    394395    prsw_ancien = 0. 
     396    prbsw_ancien= 0.
    395397
    396398    ale_wake    = 0.

  • LMDZ6/branches/lmdz-snow/libf/phylmdiso/phyetat0_mod.F90

    r5662 r5888  
    2929       falb_dir, falb_dif, prw_ancien, prlw_ancien, prsw_ancien, prbsw_ancien, &
    3030       ftsol, pbl_tke, pctsrf, q_ancien, ql_ancien, qs_ancien, qbs_ancien, &
    31        cf_ancien, rvc_ancien, radpas, radsol, rain_fall, ratqs, &
     31       cf_ancien, rvc_ancien, tke_ancien, radpas, radsol, rain_fall, ratqs, &
    3232       rnebcon, rugoro, sig1, snow_fall, bs_fall, solaire_etat0, sollw, sollwdown, &
    3333       solsw, solswfdiff, t_ancien, u_ancien, v_ancien, w01, wake_cstar, wake_deltaq, &
    … …  
    3535       wake_s, awake_s, wake_dens, awake_dens, cv_gen, zgam, zmax0, zmea, zpic, zsig, &
    3636#ifdef ISO
    37        fxtevap, xtsol, xt_ancien, xtl_ancien, xts_ancien, wake_deltaxt, &
    38        xtrain_fall,xtsnow_fall, &
     37       fxtevap, xtsol, xt_ancien, xtl_ancien, xts_ancien, xtbs_ancien, wake_deltaxt, &
     38       xtrain_fall,xtsnow_fall, xtbs_fall, &
    3939#endif
    4040       zstd, zthe, zval, ale_bl, ale_bl_trig, alp_bl, u10m, v10m, treedrg, &
    … …  
    126126    missing_val=missing_val_netcdf
    127127  ENDIF
    128 
     128 
    129129  ! FH1D
    130130  !     real iolat(jjm+1)
    … …  
    376376     bs_fall(:)=0.
    377377  ENDIF
    378 
    379 
    380378!=======================================================================
    381379! Radiation
    … …  
    495493!==================================
    496494!
     495  ! cas specifique de l'advection de TKE
     496  IF (ok_advtke) THEN
     497       ancien_ok=ancien_ok.AND.phyetat0_get(tke_ancien,"TKEANCIEN","TKEANCIEN",0.)
     498  ELSE
     499    tke_ancien(:,:)=0.
     500  ENDIF
     501
     502  IF (ok_advtke) THEN
     503    IF ( (maxval(tke_ancien).EQ.minval(tke_ancien))) THEN
     504       ancien_ok=.false.
     505    ENDIF
     506  ENDIF
     507
    497508  IF (iflag_pbl>1) then
    498509     found=phyetat0_srf(pbl_tke,"TKE","Turb. Kinetic. Energ. ",1.e-8)

  • LMDZ6/branches/lmdz-snow/libf/phylmdiso/phyredem.F90

    r5662 r5888  
    2020                                prw_ancien, prlw_ancien, prsw_ancien, prbsw_ancien,      &
    2121                                ql_ancien, qs_ancien, qbs_ancien, cf_ancien, &
    22                                 rvc_ancien, u_ancien, v_ancien,              &
     22                                rvc_ancien, u_ancien, v_ancien, tke_ancien,  &
    2323                                clwcon, rnebcon, ratqs, pbl_tke,             &
    2424                                wake_delta_pbl_tke, zmax0, f0, sig1, w01,    &
    … …  
    3636!GG
    3737#ifdef ISO
    38   USE phys_state_var_mod, ONLY: xtsol, fxtevap,xtrain_fall, xtsnow_fall,     &
    39                                 xt_ancien, xtl_ancien, xts_ancien,           &
     38  USE phys_state_var_mod, ONLY: xtsol, fxtevap,xtrain_fall, xtsnow_fall, xtbs_fall,    &
     39                                xt_ancien, xtl_ancien, xts_ancien,xtbs_ancien,  &
    4040                                wake_deltaxt                             
    4141#endif
    … …  
    310310
    311311    ! DEB TKE PBL !
     312
     313    IF (ok_advtke) THEN
     314      CALL put_field(pass,"TKEANCIEN", "TKEANCIEN", tke_ancien)
     315    ENDIF
     316
    312317
    313318    IF (iflag_pbl>1) then
    … …  
    433438       CALL put_field(pass, "dT_ns", "delta temperature near surface", dT_ns)
    434439    end if
    435 
     440   
    436441#ifdef ISO
    437442      write(*,*) 'phyredem 342'
    … …  
    571576           call iso_verif_egalite(xts_ancien(iso_eau,i,k),qs_ancien(i,k), &
    572577     &           'phyisoredem 50c')
    573          
     578           call iso_verif_egalite(xtbs_ancien(iso_eau,i,k),qbs_ancien(i,k), &                        |  ---------------------------------------------------------------------------------------------------
     579     &           'phyisoredem 50d')
    574580        enddo !do i=1,klon
    575581      enddo !do k=1,klev
    … …  
    612618      CALL put_field(pass,    "xtsnow_f"//TRIM(outiso), "precipitation solide",iso_tmp)
    613619
     620      iso_tmp(:)=xtbs_fall(ixt,:)                                                                |  ---------------------------------------------------------------------------------------------------
     621      CALL put_field(pass,    "xtbs_f"//TRIM(outiso), "precipitation neige soufflee",iso_tmp)            |  ---------------------------------------------------------------------------------------------------
     622         
    614623      iso_tmp_lonlev(:,:)=xt_ancien(ixt,:,:)
    615624      CALL put_field(pass,    "XTANCIEN"//TRIM(outiso), "QANCIEN",     iso_tmp_lonlev)
    … …  
    620629      iso_tmp_lonlev(:,:)=xts_ancien(ixt,:,:)
    621630      CALL put_field(pass,   "XTSANCIEN"//TRIM(outiso), "QSANCIEN",    iso_tmp_lonlev)
    622 
     631                                                                                                     |  ---------------------------------------------------------------------------------------------------
     632      iso_tmp_lonlev(:,:)=xtbs_ancien(ixt,:,:)                                                      |  ---------------------------------------------------------------------------------------------------
     633      CALL put_field(pass,   "XTBSANCIEN"//TRIM(outiso), "QBSANCIEN",    iso_tmp_lonlev)             |  ---------------------------------------------------------------------------------------------------
     634                               
    623635      iso_tmp_lonlev(:,:)=wake_deltaxt(ixt,:,:)
    624636      CALL put_field(pass,"WAKE_DELTAXT"//TRIM(outiso), "WAKE_DELTAQ", iso_tmp_lonlev)

  • LMDZ6/branches/lmdz-snow/libf/phylmdiso/physiq_mod.F90

    r5872 r5888  
    429429#ifdef ISO
    430430       USE phys_local_var_mod, ONLY: xt_seri,xtl_seri,xts_seri,xtbs_seri, &
    431        d_xt_eva,d_xtl_eva,d_xti_eva,           &
     431       d_xt_eva,d_xtl_eva,d_xti_eva,  &
    432432       d_xt_dyn,d_xtl_dyn,d_xts_dyn,d_xtbs_dyn, &
    433433       d_xt_con, d_xt_wake,                    &
    … …  
    10661066    REAL, dimension(klon,klev):: du0, dv0, dt0, dq0, dql0, dqi0, dqbs0
    10671067#ifdef ISO
    1068     REAL, dimension(ntraciso,klon,klev):: dxt0, dxtl0, dxti0
     1068    REAL, dimension(ntraciso,klon,klev):: dxt0, dxtl0, dxti0, dxtbs0
    10691069#endif
    10701070    REAL, dimension(klon)     :: dsig0, ddens0
    … …  
    25772577      dxtl0(:,:,:)=0.
    25782578      dxti0(:,:,:)=0.
     2579      dxtbs0(:,:,:)=0.
    25792580#endif
    25802581    dsig0(:) = 0.
    … …  
    26792680      write(*,*) 'physiq tmp 1762b: ixt,iqiso_liq=',ixt,iqIsoPha(ixt,iliq)
    26802681      if (nqo.eq.3) then 
    2681         write(*,*) 'physiq tmp 1762c: ixt,iqiso_liq=',ixt,iqIsoPha(ixt,iliq)
     2682        write(*,*) 'physiq tmp 1762c: ixt,iqiso_ice=',ixt,iqIsoPha(ixt,isol)
    26822683      endif !if (nqo.eq.3) then
     2684      if (nqo.ge.4 .and. ok_bs) then
     2685          write(*,*) 'physiq tmp 1762d: ixt,iqiso_bs=',ixt,iqIsoPha(ixt,ibs)
     2686      endif
    26832687#endif
    26842688      if (ixt.gt.niso) write(*,*) 'izone=',tracers(iqIsoPha(ixt,ivap))%iso_iZone
    … …  
    26882692          xt_seri(ixt,i,k)  = qx(i,k,iqIsoPha(ixt,ivap))
    26892693          xtl_seri(ixt,i,k) = qx(i,k,iqIsoPha(ixt,iliq))
    2690           if (nqo.eq.2) then
    2691              xts_seri(ixt,i,k) = 0.
    2692           else if (nqo.eq.3) then
     2694          xts_seri(ixt,i,k) = 0.
     2695          xtbs_seri(ixt,i,k) = 0.
     2696          if (nqo.eq.3) then
    26932697             xts_seri(ixt,i,k) = qx(i,k,iqIsoPha(ixt,isol))
     2698          else if (nqo.ge.4 .and. ok_bs) then
     2699             xtbs_seri(ixt,i,k) = qx(i,k,iqIsoPha(ixt,ibs))
    26942700          endif
     2701
    26952702       enddo !DO i = 1, klon
    26962703      enddo ! DO k = 1, klev
    … …  
    27222729        xtql1(ixt,:)=xtql1(ixt,:)+xts_seri(ixt,:,k)*zmasse(:,k)
    27232730      ENDIF
    2724       IF (ok_bs) THEN
     2731      IF (nqo >= 4 .and. ok_bs) THEN
    27252732        xtql1(ixt,:)=xtql1(ixt,:)+xtbs_seri(ixt,:,k)*zmasse(:,k)
    27262733      ENDIF
    … …  
    28042811              endif !  if (iso_verif_o18_aberrant_nostop
    28052812            endif !if (q_seri(i,k).gt.errmax) then
     2813            if (qbs_seri(i,k).gt.ridicule) then
     2814               if (iso_verif_o18_aberrant_nostop( &
     2815     &              xtbs_seri(iso_HDO,i,k)/qbs_seri(i,k), &
     2816     &              xtbs_seri(iso_O18,i,k)/qbs_seri(i,k), &
     2817     &              'physiq 2099 qs').eq.1) then
     2818                  write(*,*) 'i,k,qbs_seri(i,k)=',i,k,qbs_seri(i,k)
     2819                  !stop
     2820              endif !  if (iso_verif_o18_aberrant_nostop
     2821            endif !if (q_seri(i,k).gt.errmax) then
    28062822          enddo !k=1,klev
    28072823         enddo  !i=1,klon
    … …  
    28132829             call iso_verif_traceur(xtl_seri(1,i,k),'physiq 2620b')
    28142830             call iso_verif_traceur(xts_seri(1,i,k),'physiq 2620c')
     2831             call iso_verif_traceur(xtbs_seri(1,i,k),'physiq 2620d')
    28152832           enddo
    28162833         enddo
    … …  
    31393156     &           xts_seri(iso_eau,i,k),qs_seri(i,k), &
    31403157     &          'reevap 2209b, après réévap totale',errmax,errmaxrel)
     3158              call iso_verif_egalite_choix( &
     3159     &           xtbs_seri(iso_eau,i,k),qbs_seri(i,k), &
     3160     &          'reevap 2209c, après réévap totale',errmax,errmaxrel)
    31413161       endif !if (iso_eau.gt.0) then
    31423162     
    … …  
    31683188               'eva',abortphy,flag_inhib_tend,itap,0 &
    31693189#ifdef ISO
    3170      &    ,d_xt_eva,d_xtl_eva,d_xti_eva &
     3190     &    ,d_xt_eva,d_xtl_eva,d_xti_eva,dxtbs0 &
    31713191#endif     
    31723192     &   )
    … …  
    43474367         'convection',abortphy,flag_inhib_tend,itap,0 &
    43484368#ifdef ISO
    4349      &    ,d_xt_con,dxtl0,dxti0 &
     4369     &    ,d_xt_con,dxtl0,dxti0,dxtbs0 &
    43504370#endif     
    43514371     &   )
    … …  
    45794599            abortphy,flag_inhib_tend,itap,0 &
    45804600#ifdef ISO
    4581      &    ,d_xt_wake,dxtl0,dxti0 &
     4601     &    ,d_xt_wake,dxtl0,dxti0,dxtbs0 &
    45824602#endif     
    45834603     &   )
    … …  
    49174937                             dql0,dqi0,dqbs0,paprs,'thermals', abortphy,flag_inhib_tend,itap,0 &
    49184938#ifdef ISO
    4919      &    ,d_xt_ajs,dxtl0,dxti0 &
     4939     &    ,d_xt_ajs,dxtl0,dxti0,dxtbs0 &
    49204940#endif     
    49214941     &   )
    … …  
    50335053               'ajsb',abortphy,flag_inhib_tend,itap,0 &
    50345054#ifdef ISO
    5035      &    ,d_xt_ajsb,dxtl0,dxti0 &
     5055     &    ,d_xt_ajsb,dxtl0,dxti0,dxtbs0 &
    50365056#endif     
    50375057     &   )
    … …  
    63516371    CALL add_phys_tend(du0,dv0,d_t_swr,dq0,dql0,dqi0,dqbs0,paprs,'SW',abortphy,flag_inhib_tend,itap,0 &
    63526372#ifdef ISO
    6353      &    ,dxt0,dxtl0,dxti0 &
     6373     &    ,dxt0,dxtl0,dxti0,dxtbs0 &
    63546374#endif     
    63556375     &  )
    … …  
    63576377    CALL add_phys_tend(du0,dv0,d_t_lwr,dq0,dql0,dqi0,dqbs0,paprs,'LW',abortphy,flag_inhib_tend,itap,0 &
    63586378#ifdef ISO
    6359      &    ,dxt0,dxtl0,dxti0 &
     6379     &    ,dxt0,dxtl0,dxti0,dxtbs0 &
    63606380#endif     
    63616381     &  )
    … …  
    64446464            abortphy,flag_inhib_tend,itap,0 &
    64456465#ifdef ISO
    6446      &    ,dxt0,dxtl0,dxti0 &
     6466     &    ,dxt0,dxtl0,dxti0,dxtbs0 &
    64476467#endif     
    64486468     &   )
    … …  
    64956515            'lif', abortphy,flag_inhib_tend,itap,0 &
    64966516#ifdef ISO
    6497      &    ,dxt0,dxtl0,dxti0 &
     6517     &    ,dxt0,dxtl0,dxti0,dxtbs0 &
    64986518#endif     
    64996519     &   )
    … …  
    65236543            dqi0, dqbs0, paprs, 'hin', abortphy,flag_inhib_tend,itap,0 &
    65246544#ifdef ISO
    6525      &    ,dxt0,dxtl0,dxti0 &
     6545     &    ,dxt0,dxtl0,dxti0,dxtbs0 &
    65266546#endif     
    65276547     &   )
    … …  
    65466566            paprs, 'front_gwd_rando', abortphy,flag_inhib_tend,itap,0 &
    65476567#ifdef ISO
    6548      &    ,dxt0,dxtl0,dxti0 &
     6568     &    ,dxt0,dxtl0,dxti0,dxtbs0 &
    65496569#endif     
    65506570     &   )
    … …  
    65596579            paprs, 'flott_gwd_rando', abortphy,flag_inhib_tend,itap,0 &
    65606580#ifdef ISO
    6561      &    ,dxt0,dxtl0,dxti0 &
     6581     &    ,dxt0,dxtl0,dxti0,dxtbs0 &
    65626582#endif     
    65636583     &  )
    … …  
    66286648            'q_ch4', abortphy,flag_inhib_tend,itap,0 &
    66296649#ifdef ISO
    6630      &    ,d_xt_ch4_dtime,dxtl0,dxti0 &
     6650     &    ,d_xt_ch4_dtime,dxtl0,dxti0,dxtbs0 &
    66316651#endif     
    66326652     &   )
    … …  
    68136833              endif !  if (iso_verif_o18_aberrant_nostop
    68146834            endif !if (q_seri(i,k).gt.errmax) then
     6835
     6836            if (qbs_seri(i,k).gt.ridicule) then
     6837               if (iso_verif_o18_aberrant_nostop( &
     6838     &              xtbs_seri(iso_HDO,i,k)/qbs_seri(i,k), &
     6839     &              xtbs_seri(iso_O18,i,k)/qbs_seri(i,k), &
     6840     &              'physiq 5937d, juste apres methox, qbs').eq.1) then
     6841                  write(*,*) 'i,k,qbs_seri(i,k)=',i,k,qbs_seri(i,k)
     6842                  !stop
     6843              endif !  if (iso_verif_o18_aberrant_nostop
     6844            endif !if (q_seri(i,k).gt.errmax) then
    68156845          enddo !k=1,klev
    68166846         enddo  !i=1,klon
    … …  
    71327162        qql2(:)=qql2(:)+qs_seri(:,k)*zmasse(:,k)
    71337163      ENDIF
    7134       IF (ok_bs) THEN
     7164      IF (nqo >=4 .and. ok_bs) THEN
    71357165        qql2(:)=qql2(:)+qbs_seri(:,k)*zmasse(:,k)
    71367166      ENDIF
    … …  
    71677197          qs_seri(i,k)=qs_seri(i,k)*corrqql
    71687198        ENDIF
    7169         IF (ok_bs) THEN
     7199        IF (nqo >= 4 .and. ok_bs) THEN
    71707200          qbs_seri(i,k)=qbs_seri(i,k)*corrqql
    71717201        ENDIF
    … …  
    73167346             d_qx(i,k,iq) = ( xts_seri(ixt,i,k) - qx(i,k,iq) ) / phys_tstep
    73177347          endif
     7348          if (nqo.ge.4 .and. ok_bs) then
     7349             iq=iqIsoPha(ixt,ibs)
     7350             d_qx(i,k,iq) = ( xtbs_seri(ixt,i,k) - qx(i,k,iq) ) / phys_tstep
     7351          endif
     7352
    73187353       enddo !DO i = 1, klon
    73197354      enddo ! DO k = 1, klev
Note: See TracChangeset for help on using the changeset viewer.