Ignore:
Timestamp:
Jan 11, 2013, 10:19:19 AM (12 years ago)
Author:
Laurent Fairhead
Message:

Version testing basée sur la r1706

Testing release based on r1706

Location:
LMDZ5/branches/testing
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • LMDZ5/branches/testing

  • LMDZ5/branches/testing/libf/phylmd/physiq.F

    r1669 r1707  
    178178      save iflag_ratqs
    179179c$OMP THREADPRIVATE(iflag_ratqs)
    180       real facteur,zfratqs1,zfratqs2
     180      real facteur
    181181
    182182      REAL zz,znum,zden
    … …  
    257257c variables a une pression donnee
    258258c
    259       real rlevSTD(nlevSTD)
    260       DATA rlevSTD/100000., 92500., 85000., 70000.,
    261      .60000., 50000., 40000., 30000., 25000., 20000.,
    262      .15000., 10000., 7000., 5000., 3000., 2000., 1000./
    263       SAVE rlevstd
    264 c$OMP THREADPRIVATE(rlevstd)
    265       CHARACTER*4 clevSTD(nlevSTD)
    266       DATA clevSTD/'1000','925 ','850 ','700 ','600 ',
    267      .'500 ','400 ','300 ','250 ','200 ','150 ','100 ',
    268      .'70  ','50  ','30  ','20  ','10  '/
    269       SAVE clevSTD
    270 c$OMP THREADPRIVATE(clevSTD)
     259#include "declare_STDlev.h"
    271260c
    272261      CHARACTER*4 bb2
    273262      CHARACTER*2 bb3
    274 
    275       real twriteSTD(klon,nlevSTD,nfiles)
    276       real qwriteSTD(klon,nlevSTD,nfiles)
    277       real rhwriteSTD(klon,nlevSTD,nfiles)
    278       real phiwriteSTD(klon,nlevSTD,nfiles)
    279       real uwriteSTD(klon,nlevSTD,nfiles)
    280       real vwriteSTD(klon,nlevSTD,nfiles)
    281       real wwriteSTD(klon,nlevSTD,nfiles)
    282 cIM for NMC files
    283       REAL geo500(klon)
    284       real :: rlevSTD3(nlevSTD3)
    285       DATA rlevSTD3/85000., 50000., 25000./
    286       SAVE rlevSTD3
    287 c$OMP THREADPRIVATE(rlevSTD3)
    288       real :: rlevSTD8(nlevSTD8)
    289       DATA rlevSTD8/100000., 85000., 70000., 50000., 25000., 10000.,
    290      $     5000., 1000./
    291       SAVE rlevSTD8
    292 c$OMP THREADPRIVATE(rlevSTD8)
    293       real twriteSTD3(klon,nlevSTD3)
    294       real qwriteSTD3(klon,nlevSTD3)
    295       real rhwriteSTD3(klon,nlevSTD3)
    296       real phiwriteSTD3(klon,nlevSTD3)
    297       real uwriteSTD3(klon,nlevSTD3)
    298       real vwriteSTD3(klon,nlevSTD3)
    299       real wwriteSTD3(klon,nlevSTD3)
    300 c
    301       real tnondefSTD8(klon,nlevSTD8)
    302       real twriteSTD8(klon,nlevSTD8)
    303       real qwriteSTD8(klon,nlevSTD8)
    304       real rhwriteSTD8(klon,nlevSTD8)
    305       real phiwriteSTD8(klon,nlevSTD8)
    306       real uwriteSTD8(klon,nlevSTD8)
    307       real vwriteSTD8(klon,nlevSTD8)
    308       real wwriteSTD8(klon,nlevSTD8)
    309 c
    310 c plevSTD3 END
    311 c
    312 c nout : niveau de output des variables a une pression donnee
    313       logical oknondef(klon,nlevSTD,nout)
    314 c
    315 c les produits uvSTD, vqSTD, .., T2STD sont calcules
    316 c a partir des valeurs instantannees toutes les 6 h
    317 c qui sont moyennees sur le mois
    318263c
    319264#include "radopt.h"
    … …  
    958903      REAL snow_lsc(klon)
    959904c
    960       REAL ratqss(klon,klev),ratqsc(klon,klev)
     905      REAL ratqsc(klon,klev)
    961906      real ratqsbas,ratqshaut,tau_ratqs
    962907      save ratqsbas,ratqshaut,tau_ratqs
    … …  
    1050995      REAL zx_tmp_fi3d(klon,klev) ! variable temporaire pour champs 3D
    1051996      REAL zx_tmp_fi3d1(klon,klev+1) !variable temporaire pour champs 3D (kelvp1)
    1052 c#ifdef histNMC
    1053 cym   A voir plus tard !!!!
    1054 cym      REAL zx_tmp_NC(iim,jjmp1,nlevSTD)
    1055       REAL zx_tmp_fiNC(klon,nlevSTD)
    1056 c#endif
    1057997      REAL(KIND=8) zx_tmp2_fi3d(klon,klev) ! variable temporaire pour champs 3D
    1058998      REAL zx_tmp_2d(iim,jjmp1), zx_tmp_3d(iim,jjmp1,klev)
    1059999      REAL zx_lon(iim,jjmp1), zx_lat(iim,jjmp1)
    1060 cIM for NMC files
    1061       REAL missing_val
    1062       REAL, SAVE :: freq_moyNMC(nout)
    1063 c$OMP THREADPRIVATE(freq_moyNMC)
    10641000c
    10651001      INTEGER nid_day, nid_mth, nid_ins, nid_mthnmc, nid_daynmc
    … …  
    11371073      REAL q2m(klon,nbsrf)  ! humidite a 2m
    11381074
    1139 cIM: t2m, q2m, u10m, v10m et t2mincels, t2maxcels
     1075cIM: t2m, q2m, ustar, u10m, v10m et t2mincels, t2maxcels
    11401076      REAL zt2m(klon), zq2m(klon)             !temp., hum. 2m moyenne s/ 1 maille
    1141       REAL zu10m(klon), zv10m(klon)           !vents a 10m moyennes s/1 maille
     1077      REAL zustar(klon),zu10m(klon), zv10m(klon)  ! u* et vents a 10m moyennes s/1 maille
    11421078      CHARACTER*40 t2mincels, t2maxcels       !t2m min., t2m max
    11431079      CHARACTER*40 tinst, tave, typeval
    … …  
    12551191      integer iostat
    12561192
    1257 cIM for NMC files
    1258       missing_val=nf90_fill_real
    12591193c======================================================================
    12601194! Gestion calendrier : mise a jour du module phys_cal_mod
    … …  
    13261260      call phys_output_var_init
    13271261      print*, '================================================='
    1328 cIM for NMC files
    1329 cIM freq_moyNMC = frequences auxquelles on moyenne les champs accumules
    1330 cIM               sur les niveaux de pression standard du NMC
    1331       DO n=1, nout
    1332        freq_moyNMC(n)=freq_outNMC(n)/freq_calNMC(n)
    1333       ENDDO
    1334 c
    1335 cIM beg
     1262c
    13361263          dnwd0=0.0
    13371264          ftd=0.0
    … …  
    13811308         lalim_conv(:)=1
    13821309cRC
     1310         ustar(:,:)=0.
    13831311         u10m(:,:)=0.
    13841312         v10m(:,:)=0.
    … …  
    17681696!
    17691697      CALL change_srf_frac(itap, dtime, days_elapsed+1,
    1770      *     pctsrf, falb1, falb2, ftsol, u10m, v10m, pbl_tke)
     1698     *     pctsrf, falb1, falb2, ftsol, ustar, u10m, v10m, pbl_tke)
    17711699
    17721700
    … …  
    20782006     e     t_seri,    q_seri,    u_seri,  v_seri,   
    20792007     e     pplay,     paprs,     pctsrf,           
    2080      +     ftsol,     falb1,     falb2,   u10m,   v10m,
     2008     +     ftsol,     falb1,     falb2,   ustar, u10m,   v10m,
    20812009     s     sollwdown, cdragh,    cdragm,  u1,    v1,
    20822010     s     albsol1,   albsol2,   sens,    evap, 
    … …  
    20872015     d     s_capCL,   s_oliqCL,  s_cteiCL,s_pblT,
    20882016     d     s_therm,   s_trmb1,   s_trmb2, s_trmb3,
    2089      d     zxrugs,    zu10m,     zv10m,   fder,
     2017     d     zxrugs,    zustar, zu10m,     zv10m,   fder,
    20902018     d     zxqsurf,   rh2m,      zxfluxu, zxfluxv,
    20912019     d     frugs,     agesno,    fsollw,  fsolsw,
    … …  
    28162744
    28172745c-------------------------------------------------------------------------
    2818 c  Caclul des ratqs
    2819 c-------------------------------------------------------------------------
    2820 
    2821 c      print*,'calcul des ratqs'
    2822 c   ratqs convectifs a l'ancienne en fonction de q(z=0)-q / q
    2823 c   ----------------
    2824 c   on ecrase le tableau ratqsc calcule par clouds_gno
    2825       if (iflag_cldcon.eq.1) then
    2826          do k=1,klev
    2827          do i=1,klon
    2828             if(ptconv(i,k)) then
    2829               ratqsc(i,k)=ratqsbas
    2830      s        +fact_cldcon*(q_seri(i,1)-q_seri(i,k))/q_seri(i,k)
    2831             else
    2832                ratqsc(i,k)=0.
    2833             endif
    2834          enddo
    2835          enddo
    2836 
    2837 c-----------------------------------------------------------------------
    2838 c  par nversion de la fonction log normale
    2839 c-----------------------------------------------------------------------
    2840       else if (iflag_cldcon.eq.4) then
    2841          ptconvth(:,:)=.false.
    2842          ratqsc(:,:)=0.
    2843          if(prt_level.ge.9) print*,'avant clouds_gno thermique'
    2844          call clouds_gno
    2845      s   (klon,klev,q_seri,zqsat,clwcon0th,ptconvth,ratqsc,rnebcon0th)
    2846          if(prt_level.ge.9) print*,' CLOUDS_GNO OK'
    2847        
    2848        endif
    2849 
    2850 c   ratqs stables
    2851 c   -------------
    2852 
    2853       if (iflag_ratqs.eq.0) then
    2854 
    2855 ! Le cas iflag_ratqs=0 correspond a la version IPCC 2005 du modele.
    2856          do k=1,klev
    2857             do i=1, klon
    2858                ratqss(i,k)=ratqsbas+(ratqshaut-ratqsbas)*
    2859      s         min((paprs(i,1)-pplay(i,k))/(paprs(i,1)-30000.),1.)
    2860             enddo
    2861          enddo
    2862 
    2863 ! Pour iflag_ratqs=1 ou 2, le ratqs est constant au dessus de
    2864 ! 300 hPa (ratqshaut), varie lineariement en fonction de la pression
    2865 ! entre 600 et 300 hPa et est soit constant (ratqsbas) pour iflag_ratqs=1
    2866 ! soit lineaire (entre 0 a la surface et ratqsbas) pour iflag_ratqs=2
    2867 ! Il s'agit de differents tests dans la phase de reglage du modele
    2868 ! avec thermiques.
    2869 
    2870       else if (iflag_ratqs.eq.1) then
    2871 
    2872          do k=1,klev
    2873             do i=1, klon
    2874                if (pplay(i,k).ge.60000.) then
    2875                   ratqss(i,k)=ratqsbas
    2876                else if ((pplay(i,k).ge.30000.).and.
    2877      s            (pplay(i,k).lt.60000.)) then
    2878                   ratqss(i,k)=ratqsbas+(ratqshaut-ratqsbas)*
    2879      s            (60000.-pplay(i,k))/(60000.-30000.)
    2880                else
    2881                   ratqss(i,k)=ratqshaut
    2882                endif
    2883             enddo
    2884          enddo
    2885 
    2886       else if (iflag_ratqs.eq.2) then
    2887 
    2888          do k=1,klev
    2889             do i=1, klon
    2890                if (pplay(i,k).ge.60000.) then
    2891                   ratqss(i,k)=ratqsbas
    2892      s            *(paprs(i,1)-pplay(i,k))/(paprs(i,1)-60000.)
    2893                else if ((pplay(i,k).ge.30000.).and.
    2894      s             (pplay(i,k).lt.60000.)) then
    2895                     ratqss(i,k)=ratqsbas+(ratqshaut-ratqsbas)*
    2896      s              (60000.-pplay(i,k))/(60000.-30000.)
    2897                else
    2898                     ratqss(i,k)=ratqshaut
    2899                endif
    2900             enddo
    2901          enddo
    2902 
    2903       else if (iflag_ratqs==3) then
    2904          do k=1,klev
    2905            ratqss(:,k)=ratqsbas+(ratqshaut-ratqsbas)
    2906      s     *min( ((paprs(:,1)-pplay(:,k))/70000.)**2 , 1. )
    2907          enddo
    2908 
    2909       else if (iflag_ratqs==4) then
    2910          do k=1,klev
    2911            ratqss(:,k)=ratqsbas+0.5*(ratqshaut-ratqsbas)
    2912      s     *( tanh( (50000.-pplay(:,k))/20000.) + 1.)
    2913          enddo
    2914 
    2915       endif
    2916 
    2917 
    2918 
    2919 
    2920 c  ratqs final
    2921 c  -----------
    2922 
    2923       if (iflag_cldcon.eq.1 .or.iflag_cldcon.eq.2
    2924      s    .or.iflag_cldcon.eq.4) then
    2925 
    2926 ! On ajoute une constante au ratqsc*2 pour tenir compte de
    2927 ! fluctuations turbulentes de petite echelle
    2928 
    2929          do k=1,klev
    2930             do i=1,klon
    2931                if ((fm_therm(i,k).gt.1.e-10)) then
    2932                   ratqsc(i,k)=sqrt(ratqsc(i,k)**2+0.05**2)
    2933                endif
    2934             enddo
    2935          enddo
    2936 
    2937 !   les ratqs sont une combinaison de ratqss et ratqsc
    2938        if(prt_level.ge.9)
    2939      $       write(lunout,*)'PHYLMD NOUVEAU TAU_RATQS ',tau_ratqs
    2940 
    2941          if (tau_ratqs>1.e-10) then
    2942             facteur=exp(-pdtphys/tau_ratqs)
    2943          else
    2944             facteur=0.
    2945          endif
    2946          ratqs(:,:)=ratqsc(:,:)*(1.-facteur)+ratqs(:,:)*facteur
    2947 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    2948 ! FH 22/09/2009
    2949 ! La ligne ci-dessous faisait osciller le modele et donnait une solution
    2950 ! assymptotique bidon et dépendant fortement du pas de temps.
    2951 !        ratqs(:,:)=sqrt(ratqs(:,:)**2+ratqss(:,:)**2)
    2952 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    2953          ratqs(:,:)=max(ratqs(:,:),ratqss(:,:))
    2954       else if (iflag_cldcon<=6) then
    2955 !   on ne prend que le ratqs stable pour fisrtilp
    2956          ratqs(:,:)=ratqss(:,:)
    2957       else
    2958           zfratqs1=exp(-pdtphys/10800.)
    2959           zfratqs2=exp(-pdtphys/10800.)
    2960 !         print*,'RAPPEL RATQS 1 ',zfratqs1,zfratqs2
    2961 !    s    ,ratqss(1,14),ratqs(1,14),ratqsc(1,14)
    2962           do k=1,klev
    2963              do i=1,klon
    2964                 if (ratqsc(i,k).gt.1.e-10) then
    2965                    ratqs(i,k)=ratqs(i,k)*zfratqs2
    2966      s             +(iflag_cldcon/100.)*ratqsc(i,k)*(1.-zfratqs2)
    2967                 endif
    2968                 ratqs(i,k)=min(ratqs(i,k)*zfratqs1
    2969      s          +ratqss(i,k)*(1.-zfratqs1),0.5)
    2970              enddo
    2971           enddo
    2972       endif
     2746! Computation of ratqs, the width (normalized) of the subrid scale
     2747! water distribution
     2748      CALL  calcratqs(klon,klev,prt_level,lunout,       
     2749     s     iflag_ratqs,iflag_con,iflag_cldcon,pdtphys,
     2750     s     ratqsbas,ratqshaut,tau_ratqs,fact_cldcon, 
     2751     s     ptconv,ptconvth,clwcon0th, rnebcon0th,   
     2752     s     paprs,pplay,q_seri,zqsat,fm_therm,
     2753     s     ratqs,ratqsc)
    29732754
    29742755
    … …  
    38433624     I     cdragh,   coefh,     fm_therm, entr_therm,
    38443625     I     u1,       v1,        ftsol,    pctsrf,
     3626     I     ustar,     u10m,      v10m,
    38453627     I     rlat,     frac_impa, frac_nucl,rlon,
    38463628     I     presnivs, pphis,     pphi,     albsol1,
    … …  
    39333715c
    39343716#include "calcul_STDlev.h"
    3935       twriteSTD(:,:,1)=tsumSTD(:,:,1)
    3936       qwriteSTD(:,:,1)=qsumSTD(:,:,1)
    3937       rhwriteSTD(:,:,1)=rhsumSTD(:,:,1)
    3938       phiwriteSTD(:,:,1)=phisumSTD(:,:,1)
    3939       uwriteSTD(:,:,1)=usumSTD(:,:,1)
    3940       vwriteSTD(:,:,1)=vsumSTD(:,:,1)
    3941       wwriteSTD(:,:,1)=wsumSTD(:,:,1)
    3942 
    3943       twriteSTD(:,:,2)=tsumSTD(:,:,2)
    3944       qwriteSTD(:,:,2)=qsumSTD(:,:,2)
    3945       rhwriteSTD(:,:,2)=rhsumSTD(:,:,2)
    3946       phiwriteSTD(:,:,2)=phisumSTD(:,:,2)
    3947       uwriteSTD(:,:,2)=usumSTD(:,:,2)
    3948       vwriteSTD(:,:,2)=vsumSTD(:,:,2)
    3949       wwriteSTD(:,:,2)=wsumSTD(:,:,2)
    3950 
    3951       twriteSTD(:,:,3)=tlevSTD(:,:)
    3952       qwriteSTD(:,:,3)=qlevSTD(:,:)
    3953       rhwriteSTD(:,:,3)=rhlevSTD(:,:)
    3954       phiwriteSTD(:,:,3)=philevSTD(:,:)
    3955       uwriteSTD(:,:,3)=ulevSTD(:,:)
    3956       vwriteSTD(:,:,3)=vlevSTD(:,:)
    3957       wwriteSTD(:,:,3)=wlevSTD(:,:)
    3958 
    3959       twriteSTD(:,:,4)=tlevSTD(:,:)
    3960       qwriteSTD(:,:,4)=qlevSTD(:,:)
    3961       rhwriteSTD(:,:,4)=rhlevSTD(:,:)
    3962       phiwriteSTD(:,:,4)=philevSTD(:,:)
    3963       uwriteSTD(:,:,4)=ulevSTD(:,:)
    3964       vwriteSTD(:,:,4)=vlevSTD(:,:)
    3965       wwriteSTD(:,:,4)=wlevSTD(:,:)
    3966 c
    3967 cIM initialisation 5eme fichier de sortie
    3968       twriteSTD(:,:,5)=tlevSTD(:,:)
    3969       qwriteSTD(:,:,5)=qlevSTD(:,:)
    3970       rhwriteSTD(:,:,5)=rhlevSTD(:,:)
    3971       phiwriteSTD(:,:,5)=philevSTD(:,:)
    3972       uwriteSTD(:,:,5)=ulevSTD(:,:)
    3973       vwriteSTD(:,:,5)=vlevSTD(:,:)
    3974       wwriteSTD(:,:,5)=wlevSTD(:,:)
    3975 c
    3976 cIM initialisation 6eme fichier de sortie
    3977       twriteSTD(:,:,6)=tlevSTD(:,:)
    3978       qwriteSTD(:,:,6)=qlevSTD(:,:)
    3979       rhwriteSTD(:,:,6)=rhlevSTD(:,:)
    3980       phiwriteSTD(:,:,6)=philevSTD(:,:)
    3981       uwriteSTD(:,:,6)=ulevSTD(:,:)
    3982       vwriteSTD(:,:,6)=vlevSTD(:,:)
    3983       wwriteSTD(:,:,6)=wlevSTD(:,:)
    3984 cIM for NMC files
    3985       DO n=1, nlevSTD3
    3986        DO k=1, nlevSTD
    3987         if(rlevSTD3(n).EQ.rlevSTD(k)) THEN
    3988          twriteSTD3(:,n)=tlevSTD(:,k)
    3989          qwriteSTD3(:,n)=qlevSTD(:,k)
    3990          rhwriteSTD3(:,n)=rhlevSTD(:,k)
    3991          phiwriteSTD3(:,n)=philevSTD(:,k)
    3992          uwriteSTD3(:,n)=ulevSTD(:,k)
    3993          vwriteSTD3(:,n)=vlevSTD(:,k)
    3994          wwriteSTD3(:,n)=wlevSTD(:,k)
    3995         endif !rlevSTD3(n).EQ.rlevSTD(k)
    3996        ENDDO
    3997       ENDDO
    3998 c
    3999       DO n=1, nlevSTD8
    4000        DO k=1, nlevSTD
    4001         if(rlevSTD8(n).EQ.rlevSTD(k)) THEN
    4002          tnondefSTD8(:,n)=tnondef(:,k,2)
    4003          twriteSTD8(:,n)=tsumSTD(:,k,2)
    4004          qwriteSTD8(:,n)=qsumSTD(:,k,2)
    4005          rhwriteSTD8(:,n)=rhsumSTD(:,k,2)
    4006          phiwriteSTD8(:,n)=phisumSTD(:,k,2)
    4007          uwriteSTD8(:,n)=usumSTD(:,k,2)
    4008          vwriteSTD8(:,n)=vsumSTD(:,k,2)
    4009          wwriteSTD8(:,n)=wsumSTD(:,k,2)
    4010         endif !rlevSTD8(n).EQ.rlevSTD(k)
    4011        ENDDO
    4012       ENDDO
    40133717c
    40143718c slp sea level pressure
Note: See TracChangeset for help on using the changeset viewer.