Ignore:
Timestamp:
Jun 9, 2011, 12:13:33 AM (13 years ago)
Author:
musat
Message:

Ajouts CFMIP2/CMIP5

  • 6eme fichier de sortie "stations" histstn.nc qui necessite 2 fichiers (voir DefLists?): npCFMIP_param.data(_*) contenant le nombre de points (120 pour simulations AMIP, 73 pour aqua) pointlocations.txt(_*) contenant le numero, les coordonnees (lon,lat) et le nom de chaque station
  • flag LOGICAL dans tous les appels histwrite_phy pour pouvoir sortir le fichier histstn.nc

NB: 1) les flags de type phys_ que l'on met dans le physiq.def_L* pour ajouter plus de sorties

necessitent dorenavant 6 valeurs, la 6eme correspondant au fichier histstn.nc

2) par defaut le fichier histstn.nc ne sort pas; pour le sortir ajouter les lignes suivantes

dans physiq.def_L*

### Type de fichier : global (n) ou stations (y)
phys_out_filestations = n n n n n y

  • introduction de 2 jeux de flags pour les taux des GES; taux actuels avec suffixes _act, taux futurs avec "_per" avec 2 appels au rayonnement si taux "_per" different des taux "_act" (utiles pour diags. CFMIP 4CO2)
  • flags "betaCRF" pour calculs CRF pour experiences sensibilite proprietes optiques eau liquide nuageuse avec initialisations par defaut; sinon besoin de fichier beta_crf.data

Ajout flag LOGICAL lCOSP necessaire pour sortir un fichier stations
IM

File:
1 edited

Legend:

Unmodified
Added
Removed

  • LMDZ5/trunk/libf/phylmd/physiq.F

    r1538 r1539  
    4444
    4545
     46!IM stations CFMIP
     47      USE CFMIP_point_locations
    4648      IMPLICIT none
    4749c======================================================================
    … …  
    699701cAA
    700702      REAL coefh(klon,klev)     ! coef d'echange pour phytrac, valable pour 2<=k<=klev
     703      REAL coefm(klon,klev)     ! coef d'echange pour U, V
    701704      REAL u1(klon)             ! vents dans la premiere couche U
    702705      REAL v1(klon)             ! vents dans la premiere couche V
    … …  
    10111014      REAL zx_tmp_fi2d(klon)      ! variable temporaire grille physique
    10121015      REAL zx_tmp_fi3d(klon,klev) ! variable temporaire pour champs 3D
     1016      REAL zx_tmp_fi3d1(klon,klev+1) !variable temporaire pour champs 3D (kelvp1)
    10131017c#ifdef histNMC
    10141018cym   A voir plus tard !!!!
    … …  
    10411045      REAL ue_lay(klon,klev) ! transport zonal de l'energie a chaque niveau vert.
    10421046      REAL uq_lay(klon,klev) ! transport zonal de l'eau a chaque niveau vert.
    1043 c
    1044 cIM 280405 END
    10451047c
    10461048      INTEGER nhori, nvert, nvert1, nvert3
    … …  
    11881190      REAL grain(1), gtsol(1), gt2m(1), gprw(1)
    11891191
     1192cIM stations CFMIP
     1193      INTEGER, SAVE :: nCFMIP
     1194c$OMP THREADPRIVATE(nCFMIP)
     1195      INTEGER, PARAMETER :: npCFMIP=120
     1196      INTEGER, ALLOCATABLE, SAVE :: tabCFMIP(:)
     1197      REAL, ALLOCATABLE, SAVE :: lonCFMIP(:), latCFMIP(:)
     1198c$OMP THREADPRIVATE(tabCFMIP, lonCFMIP, latCFMIP)
     1199      INTEGER, ALLOCATABLE, SAVE :: tabijGCM(:)
     1200      REAL, ALLOCATABLE, SAVE :: lonGCM(:), latGCM(:)
     1201c$OMP THREADPRIVATE(tabijGCM, lonGCM, latGCM)
     1202      INTEGER, ALLOCATABLE, SAVE :: iGCM(:), jGCM(:)
     1203c$OMP THREADPRIVATE(iGCM, jGCM)
     1204      logical, dimension(nfiles)            :: phys_out_filestations
     1205      logical, parameter :: lNMC=.FALSE.
     1206
     1207cIM betaCRF
     1208      REAL, SAVE :: pfree, beta_pbl, beta_free
     1209c$OMP THREADPRIVATE(pfree, beta_pbl, beta_free)
     1210      REAL, SAVE :: lon1_beta,  lon2_beta, lat1_beta, lat2_beta
     1211c$OMP THREADPRIVATE(lon1_beta,  lon2_beta, lat1_beta, lat2_beta)
     1212      LOGICAL, SAVE :: mskocean_beta
     1213c$OMP THREADPRIVATE(mskocean_beta)
     1214      REAL, dimension(klon, klev) :: beta       ! facteur sur cldtaurad et cldemirad pour evaluer les retros liees aux CRF
     1215      REAL, dimension(klon, klev) :: cldtaurad  ! epaisseur optique pour radlwsw,COSP
     1216      REAL, dimension(klon, klev) :: cldemirad  ! emissivite pour radlwsw,COSP
     1217
    11901218cIM for NMC files
    11911219      missing_val=nf90_fill_real
    … …  
    14721500
    14731501c================================================================================
    1474 
     1502cIM stations CFMIP
     1503      nCFMIP=npCFMIP
     1504      OPEN(98,file='npCFMIP_param.data',status='old',
     1505     $          form='formatted',err=999)
     1506      READ(98,*,end=998) nCFMIP
     1507998   CONTINUE
     1508      CLOSE(98)
     1509      CONTINUE
     1510      IF(nCFMIP.GT.npCFMIP) THEN
     1511       print*,'nCFMIP > npCFMIP : augmenter npCFMIP et recompiler'
     1512       CALL abort
     1513      else
     1514       print*,'physiq npCFMIP=',npCFMIP,'nCFMIP=',nCFMIP
     1515      ENDIF
     1516c
     1517      ALLOCATE(tabCFMIP(nCFMIP))
     1518      ALLOCATE(lonCFMIP(nCFMIP), latCFMIP(nCFMIP))
     1519      ALLOCATE(tabijGCM(nCFMIP))
     1520      ALLOCATE(lonGCM(nCFMIP), latGCM(nCFMIP))
     1521      ALLOCATE(iGCM(nCFMIP), jGCM(nCFMIP))
     1522c
     1523c lecture des nCFMIP stations CFMIP, de leur numero
     1524c et des coordonnees geographiques lonCFMIP, latCFMIP
     1525c
     1526         CALL read_CFMIP_point_locations(nCFMIP, tabCFMIP,
     1527     $lonCFMIP, latCFMIP)
     1528c
     1529c identification des
     1530c 1) coordonnees lonGCM, latGCM des points CFMIP dans la grille de LMDZ
     1531c 2) indices points tabijGCM de la grille physique 1d sur klon points
     1532c 3) indices iGCM, jGCM de la grille physique 2d
     1533c
     1534         CALL LMDZ_CFMIP_point_locations(nCFMIP, lonCFMIP, latCFMIP,
     1535     $tabijGCM, lonGCM, latGCM, iGCM, jGCM)
     1536c
     1537999      CONTINUE
    14751538         ENDIF !debut
    1476 
     1539 
    14771540           DO i=1,klon
    14781541             rugoro(i) = f_rugoro * MAX(1.0e-05, zstd(i)*zsig(i)/2.0)
    … …  
    15121575     .                   lmt_pas
    15131576c
    1514 cIM 030306 END
    1515 
    15161577      capemaxcels = 't_max(X)'
    15171578      t2mincels = 't_min(X)'
    … …  
    15301591
    15311592c$OMP MASTER
    1532        call phys_output_open(jjmp1,nlevSTD,clevSTD,nbteta,
    1533      &                        ctetaSTD,dtime,ok_veget,
    1534      &                        type_ocean,iflag_pbl,ok_mensuel,ok_journe,
    1535      &                        ok_hf,ok_instan,ok_LES,ok_ade,ok_aie,
    1536      &                        read_climoz, new_aod, aerosol_couple
     1593       call phys_output_open(rlon,rlat,nCFMIP,tabijGCM,
     1594     &                       iGCM,jGCM,lonGCM,latGCM,
     1595     &                       jjmp1,nlevSTD,clevSTD,
     1596     &                       nbteta, ctetaSTD, dtime,ok_veget,
     1597     &                       type_ocean,iflag_pbl,ok_mensuel,ok_journe,
     1598     &                       ok_hf,ok_instan,ok_LES,ok_ade,ok_aie,
     1599     &                       read_climoz, phys_out_filestations,
     1600     &                       new_aod, aerosol_couple
    15371601     &                        )
    15381602c$OMP END MASTER
    … …  
    15551619#endif
    15561620
    1557 cIM 250308bad guide        ecrit_hf2mth = 30*1/ecrit_hf
    15581621         ecrit_hf2mth = ecrit_mth/ecrit_hf
    15591622
    … …  
    15681631         ecrit_reg = ecrit_reg * un_jour
    15691632         ecrit_tra = ecrit_tra * un_jour
    1570          ecrit_ISCCP = ecrit_ISCCP * un_jour
    15711633         ecrit_LES = ecrit_LES * un_jour
    15721634c
    … …  
    15741636     .   ecrit_hf,ecrit_day,ecrit_mth,ecrit_reg,ecrit_tra,ecrit_ISCCP,
    15751637     .   ecrit_hf2mth
    1576 cIM 030306 END
    1577 
    15781638
    15791639cXXXPB Positionner date0 pour initialisation de ORCHIDEE
    … …  
    16331693      END IF
    16341694C$omp end single
     1695c
     1696cIM betaCRF
     1697      pfree=70000. !Pa
     1698      beta_pbl=1.
     1699      beta_free=1.
     1700      lon1_beta=-180.
     1701      lon2_beta=+180.
     1702      lat1_beta=90.
     1703      lat2_beta=-90.
     1704      mskocean_beta=.FALSE.
     1705
     1706      OPEN(99,file='beta_crf.data',status='old',
     1707     $          form='formatted',err=9999)
     1708      READ(99,*,end=9998) pfree
     1709      READ(99,*,end=9998) beta_pbl
     1710      READ(99,*,end=9998) beta_free
     1711      READ(99,*,end=9998) lon1_beta
     1712      READ(99,*,end=9998) lon2_beta
     1713      READ(99,*,end=9998) lat1_beta
     1714      READ(99,*,end=9998) lat2_beta
     1715      READ(99,*,end=9998) mskocean_beta
     17169998  Continue
     1717      CLOSE(99)
     17189999  Continue
     1719      WRITE(*,*)'pfree=',pfree
     1720      WRITE(*,*)'beta_pbl=',beta_pbl
     1721      WRITE(*,*)'beta_free=',beta_free
     1722      WRITE(*,*)'lon1_beta=',lon1_beta
     1723      WRITE(*,*)'lon2_beta=',lon2_beta
     1724      WRITE(*,*)'lat1_beta=',lat1_beta
     1725      WRITE(*,*)'lat2_beta=',lat2_beta
     1726      WRITE(*,*)'mskocean_beta=',mskocean_beta
    16351727      ENDIF
    16361728!
    … …  
    19492041     s     zxtsol,    zxfluxlat, zt2m,    qsat2m,
    19502042     s     d_t_vdf,   d_q_vdf,   d_u_vdf, d_v_vdf,
    1951      s     coefh,     slab_wfbils,               
     2043     s     coefh,     coefm,     slab_wfbils,               
    19522044     d     qsol,      zq2m,      s_pblh,  s_lcl,
    19532045     d     s_capCL,   s_oliqCL,  s_cteiCL,s_pblT,
    … …  
    19652057! ajout des tendances de la diffusion turbulente
    19662058      CALL add_phys_tend(d_u_vdf,d_v_vdf,d_t_vdf,d_q_vdf,dql0,'vdf')
    1967 
    19682059!-----------------------------------------------------------------------------------------
    19692060
    … …  
    21032194            enddo
    21042195         enddo
    2105 
    21062196     
    21072197cc--   Calcul de l'energie disponible ALE (J/kg) et de la puissance disponible ALP (W/m2)
    … …  
    25242614
    25252615
    2526 
    25272616c  Ajustement sec
    25282617c  ==============
    … …  
    27402829
    27412830
    2742 
    27432831c
    27442832c Appeler le processus de condensation a grande echelle
    … …  
    28732961     &        tausum_aero, tau3d_aero)
    28742962      ELSE
    2875 cIM 170310 BEG
    28762963         tausum_aero(:,:,:) = 0.
    2877 cIM 170310 END
    28782964         tau_aero(:,:,:,:) = 0.
    28792965         piz_aero(:,:,:,:) = 0.
    … …  
    31493235      endif
    31503236c
     3237cIM betaCRF
     3238c
     3239      cldtaurad = cldtau
     3240      cldemirad = cldemi
     3241c
     3242      if(lon1_beta.EQ.-180..AND.lon2_beta.EQ.180..AND.
     3243     $lat1_beta.EQ.90..AND.lat2_beta.EQ.-90.) THEN
     3244c
     3245c global
     3246c
     3247       DO k=1, klev
     3248       DO i=1, klon
     3249        if (pplay(i,k).GE.pfree) THEN
     3250         beta(i,k) = beta_pbl
     3251        else
     3252         beta(i,k) = beta_free
     3253        endif
     3254        if (mskocean_beta) THEN
     3255         beta(i,k) = beta(i,k) * pctsrf(i,is_oce)
     3256        endif
     3257        cldtaurad(i,k) = cldtau(i,k) * beta(i,k)
     3258        cldemirad(i,k) = cldemi(i,k) * beta(i,k)
     3259       ENDDO
     3260       ENDDO
     3261c
     3262      else
     3263c
     3264c regional
     3265c
     3266       DO k=1, klev
     3267       DO i=1,klon
     3268c
     3269        if (rlon(i).ge.lon1_beta.AND.rlon(i).le.lon2_beta.AND.
     3270     $      rlat(i).le.lat1_beta.AND.rlat(i).ge.lat2_beta) THEN
     3271         if (pplay(i,k).GE.pfree) THEN
     3272          beta(i,k) = beta_pbl
     3273         else
     3274          beta(i,k) = beta_free
     3275         endif
     3276         if (mskocean_beta) THEN
     3277          beta(i,k) = beta(i,k) * pctsrf(i,is_oce)
     3278         endif
     3279        cldtaurad(i,k) = cldtau(i,k) * beta(i,k)
     3280        cldemirad(i,k) = cldemi(i,k) * beta(i,k)
     3281        endif
     3282c
     3283       ENDDO
     3284       ENDDO
     3285c
     3286      endif
     3287c
    31513288c Appeler le rayonnement mais calculer tout d'abord l'albedo du sol.
    31523289c
    … …  
    31773314     e        paprs, pplay,zxtsol,albsol1, albsol2, t_seri,q_seri,
    31783315     e        wo(:, :, 1),
    3179      e        cldfra, cldemi, cldtau,
     3316     e        cldfra, cldemirad, cldtaurad,
    31803317     s        heat,heat0,cool,cool0,radsol,albpla,
    31813318     s        topsw,toplw,solsw,sollw,
    … …  
    31953332#endif
    31963333      ELSE
    3197 
     3334c
     3335cIM calcul radiatif pour le cas actuel
     3336c
     3337       RCO2 = RCO2_act
     3338       RCH4 = RCH4_act
     3339       RN2O = RN2O_act
     3340       RCFC11 = RCFC11_act
     3341       RCFC12 = RCFC12_act
     3342c
    31983343         CALL radlwsw
    31993344     e        (dist, rmu0, fract,
    32003345     e        paprs, pplay,zxtsol,albsol1, albsol2,
    32013346     e        t_seri,q_seri,wo,
    3202      e        cldfra, cldemi, cldtau,
     3347     e        cldfra, cldemirad, cldtaurad,
    32033348     e        ok_ade, ok_aie,
    32043349     e        tau_aero, piz_aero, cg_aero,
    … …  
    32183363     o        topswcf_aero, solswcf_aero)
    32193364         
    3220 
     3365c
     3366cIM 2eme calcul radiatif pour le cas perturbe ou au moins un
     3367cIM des taux doit etre different du taux actuel
     3368cIM Par defaut on a les taux perturbes egaux aux taux actuels
     3369c
     3370       if (RCO2_per.NE.RCO2_act.OR.RCH4_per.NE.RCH4_act.OR.
     3371     $RN2O_per.NE.RN2O_act.OR.RCFC11_per.NE.RCFC11_act.OR.
     3372     $RCFC12_per.NE.RCFC12_act) THEN
     3373c
     3374       RCO2 = RCO2_per
     3375       RCH4 = RCH4_per
     3376       RN2O = RN2O_per
     3377       RCFC11 = RCFC11_per
     3378       RCFC12 = RCFC12_per
     3379c
     3380         CALL radlwsw
     3381     e        (dist, rmu0, fract,
     3382     e        paprs, pplay,zxtsol,albsol1, albsol2,
     3383     e        t_seri,q_seri,wo,
     3384     e        cldfra, cldemi, cldtau,
     3385     e        ok_ade, ok_aie,
     3386     e        tau_aero, piz_aero, cg_aero,
     3387     e        cldtaupi,new_aod,
     3388     e        zqsat, flwc, fiwc,
     3389     s        heatp,heat0p,coolp,cool0p,radsolp,albplap,
     3390     s        topswp,toplwp,solswp,sollwp,
     3391     s        sollwdownp,
     3392     s        topsw0p,toplw0p,solsw0p,sollw0p,
     3393     s        lwdn0p, lwdnp, lwup0p, lwupp,
     3394     s        swdn0p, swdnp, swup0p, swupp,
     3395     s        topswad_aerop, solswad_aerop,
     3396     s        topswai_aerop, solswai_aerop,
     3397     o        topswad0_aerop, solswad0_aerop,
     3398     o        topsw_aerop, topsw0_aerop,
     3399     o        solsw_aerop, solsw0_aerop,
     3400     o        topswcf_aerop, solswcf_aerop)
     3401       endif
     3402c
    32213403      ENDIF ! aerosol_couple
    32223404      itaprad = 0
    … …  
    33843566c
    33853567c  ajout des tendances
    3386         CALL add_phys_tend(d_u_hin,d_v_hin,d_t_hin,dq0,dql0,'lif')
     3568        CALL add_phys_tend(d_u_hin,d_v_hin,d_t_hin,dq0,dql0,'hin')
    33873569
    33883570      ENDIF
    … …  
    33973579        call writefield_phy('v_seri',v_seri,llm)
    33983580        call writefield_phy('t_seri',t_seri,llm)
    3399         call writefield_phy('q_seri',q_seri,llm)
     3581        call writefield_phy('q_seri',q_seri,llm)
    34003582      endif
    34013583
    … …  
    34653647     $                   prfl(:,1:klev),psfl(:,1:klev),
    34663648     $                   pmflxr(:,1:klev),pmflxs(:,1:klev),
    3467      $                   mr_ozone,cldtau, cldemi)
     3649     $                   mr_ozone,cldtaurad, cldemirad)
    34683650
    34693651!     L          calipso2D,calipso3D,cfadlidar,parasolrefl,atb,betamol,
    … …  
    36163798c
    36173799cIM initialisation 5eme fichier de sortie
    3618 cIM ajoute 5eme niveau 170310 BEG
    36193800      twriteSTD(:,:,5)=tlevSTD(:,:)
    36203801      qwriteSTD(:,:,5)=qlevSTD(:,:)
    … …  
    36243805      vwriteSTD(:,:,5)=vlevSTD(:,:)
    36253806      wwriteSTD(:,:,5)=wlevSTD(:,:)
     3807c
     3808cIM initialisation 6eme fichier de sortie
     3809      twriteSTD(:,:,6)=tlevSTD(:,:)
     3810      qwriteSTD(:,:,6)=qlevSTD(:,:)
     3811      rhwriteSTD(:,:,6)=rhlevSTD(:,:)
     3812      phiwriteSTD(:,:,6)=philevSTD(:,:)
     3813      uwriteSTD(:,:,6)=ulevSTD(:,:)
     3814      vwriteSTD(:,:,6)=vlevSTD(:,:)
     3815      wwriteSTD(:,:,6)=wlevSTD(:,:)
    36263816cIM for NMC files
    36273817      DO n=1, nlevSTD3
Note: See TracChangeset for help on using the changeset viewer.