Ignore:
Timestamp:
Oct 22, 2010, 6:18:27 PM (14 years ago)
Author:
jghattas
Message:
  • Added variables written to file phystokenc.nc by option offline.
  • initphysto and phystokenc rewritten in F90
  • ener.h modified to be compatible with F77 and F90 syntax
File:
1 moved

Legend:

Unmodified
Added
Removed

  • LMDZ5/branches/LMDZ5V1.0-dev/libf/phylmd/phystokenc.F90

    r1436 r1447  
    1 !
    2 c
    3 c
    4       SUBROUTINE phystokenc (
    5      I                   nlon,nlev,pdtphys,rlon,rlat,
    6      I                   pt,pmfu, pmfd, pen_u, pde_u, pen_d, pde_d,
    7      I                   pfm_therm,pentr_therm,
    8      I                   cdragh, pcoefh,yu1,yv1,ftsol,pctsrf,
    9      I                   frac_impa,frac_nucl,
    10      I                   pphis,paire,dtime,itap)
    11       USE ioipsl
    12       USE dimphy
    13       USE infotrac, ONLY : nqtot
    14       USE iophy
    15       USE control_mod
    16 
    17       IMPLICIT none
    18 
    19 c======================================================================
    20 c Auteur(s) FH
    21 c Objet: Moniteur general des tendances traceurs
    22 c
    23 
    24 c======================================================================
    25 #include "dimensions.h"
    26 #include "tracstoke.h"
    27 #include "indicesol.h"
    28 c======================================================================
    29 
    30 c Arguments:
    31 c
    32 c   EN ENTREE:
    33 c   ==========
    34 c
    35 c   divers:
    36 c   -------
    37 c
    38       integer nlon ! nombre de points horizontaux
    39       integer nlev ! nombre de couches verticales
    40       real pdtphys ! pas d'integration pour la physique (seconde)
    41 c
    42       integer physid, itap
    43       save physid
    44 c$OMP THREADPRIVATE(physid)
    45       integer ndex2d(iim*(jjm+1)),ndex3d(iim*(jjm+1)*klev)
    46 
    47 c   convection:
    48 c   -----------
    49 c
    50       REAL pmfu(klon,klev)  ! flux de masse dans le panache montant
    51       REAL pmfd(klon,klev)  ! flux de masse dans le panache descendant
    52       REAL pen_u(klon,klev) ! flux entraine dans le panache montant
    53       REAL pde_u(klon,klev) ! flux detraine dans le panache montant
    54       REAL pen_d(klon,klev) ! flux entraine dans le panache descendant
    55       REAL pde_d(klon,klev) ! flux detraine dans le panache descendant
    56       real pt(klon,klev)
    57       REAL,allocatable,save :: t(:,:)
    58 c$OMP THREADPRIVATE(t)
    59 c
    60       REAL rlon(klon), rlat(klon), dtime
    61       REAL zx_tmp_3d(iim,jjm+1,klev),zx_tmp_2d(iim,jjm+1)
    62 
    63 c   Couche limite:
    64 c   --------------
    65 c
    66       REAL cdragh(klon)          ! cdrag
    67       REAL pcoefh(klon,klev)     ! coeff melange CL
    68       REAL pcoefh_buf(klon,klev) ! coeff melange CL + cdrag
    69       REAL yv1(klon)
    70       REAL yu1(klon),pphis(klon),paire(klon)
    71 
    72 c   Les Thermiques : (Abderr 25 11 02)
    73 c   ---------------
    74       REAL pfm_therm(klon,klev+1)
    75       real fm_therm1(klon,klev)
    76       REAL pentr_therm(klon,klev)
    77    
    78       REAL,allocatable,save :: entr_therm(:,:)
    79       REAL,allocatable,save :: fm_therm(:,:)
    80 c$OMP THREADPRIVATE(entr_therm)
    81 c$OMP THREADPRIVATE(fm_therm)
    82 c
    83 c   Lessivage:
    84 c   ----------
    85 c
    86       REAL frac_impa(klon,klev)
    87       REAL frac_nucl(klon,klev)
    88 c
    89 c Arguments necessaires pour les sources et puits de traceur
    90 C
    91       real ftsol(klon,nbsrf)  ! Temperature du sol (surf)(Kelvin)
    92       real pctsrf(klon,nbsrf) ! Pourcentage de sol f(nature du sol)
    93 c======================================================================
    94 c
    95       INTEGER i, k
    96 c
    97       REAL,allocatable,save :: mfu(:,:)  ! flux de masse dans le panache montant
    98       REAL,allocatable,save :: mfd(:,:)  ! flux de masse dans le panache descendant
    99       REAL,allocatable,save :: en_u(:,:) ! flux entraine dans le panache montant
    100       REAL,allocatable,save :: de_u(:,:) ! flux detraine dans le panache montant
    101       REAL,allocatable,save :: en_d(:,:) ! flux entraine dans le panache descendant
    102       REAL,allocatable,save :: de_d(:,:) ! flux detraine dans le panache descendant
    103       REAL,allocatable,save :: coefh(:,:) ! flux detraine dans le panache descendant
    104 
    105       REAL,allocatable,save :: pyu1(:)
    106       REAL,allocatable,save :: pyv1(:)
    107       REAL,allocatable,save :: pftsol(:,:)
    108       REAL,allocatable,save :: ppsrf(:,:)
    109 c$OMP THREADPRIVATE(mfu,mfd,en_u,de_u,en_d,de_d,coefh)
    110 c$OMP THREADPRIVATE(pyu1,pyv1,pftsol,ppsrf)
    111       real pftsol1(klon),pftsol2(klon),pftsol3(klon),pftsol4(klon)
    112       real ppsrf1(klon),ppsrf2(klon),ppsrf3(klon),ppsrf4(klon)
    113 
    114       REAL dtcum
    115 
    116       integer iadvtr,irec
    117       real zmin,zmax
    118       logical ok_sync
    119  
    120       save dtcum
    121       save iadvtr,irec
    122 c$OMP THREADPRIVATE(dtcum,iadvtr,irec)
    123       data iadvtr,irec/0,1/
    124       logical,save :: first=.true.
    125 c$OMP THREADPRIVATE(first)
    126 c
    127 c   Couche limite:
    128 c======================================================================
    129 
    130 c Dans le meme vecteur on recombine le drag et les coeff d'echange
    131       pcoefh_buf(:,1)      = cdragh(:)
    132       pcoefh_buf(:,2:klev) = pcoefh(:,2:klev)
    133 
    134       ok_sync = .true.
    135         print*,'Dans phystokenc.F'
    136       print*,'iadvtr= ',iadvtr
    137       print*,'istphy= ',istphy
    138       print*,'istdyn= ',istdyn
    139 
    140       if (first) then
    141      
    142         allocate( t(klon,klev))
    143         allocate( mfu(klon,klev)) 
    144         allocate( mfd(klon,klev)) 
    145         allocate( en_u(klon,klev))
    146         allocate( de_u(klon,klev))
    147         allocate( en_d(klon,klev))
    148         allocate( de_d(klon,klev))
    149         allocate( coefh(klon,klev))
    150         allocate( entr_therm(klon,klev))
    151         allocate( fm_therm(klon,klev))
    152         allocate( pyu1(klon))
    153         allocate( pyv1(klon))
    154         allocate( pftsol(klon,nbsrf))
    155         allocate( ppsrf(klon,nbsrf))
    156  
    157         first=.false.
    158       endif
    159      
    160       IF (iadvtr.eq.0) THEN
    161        
    162         CALL initphysto('phystoke',
    163      . rlon,rlat,dtime, dtime*istphy,dtime*istphy,nqtot,physid)
    164        
    165         write(*,*) 'apres initphysto ds phystokenc'
    166 
    167        
    168       ENDIF
    169 c
    170       ndex2d = 0
    171       ndex3d = 0
    172       i=itap
    173 cym      CALL gr_fi_ecrit(1,klon,iim,jjm+1,pphis,zx_tmp_2d)
    174       CALL histwrite_phy(physid,"phis",i,pphis)
    175 c
    176       i=itap
    177 cym      CALL gr_fi_ecrit(1,klon,iim,jjm+1,paire,zx_tmp_2d)
    178       CALL histwrite_phy(physid,"aire",i,paire)
    179 
    180       iadvtr=iadvtr+1
    181 c
    182       if (mod(iadvtr,istphy).eq.1.or.istphy.eq.1) then
    183         print*,'reinitialisation des champs cumules
    184      s          a iadvtr=',iadvtr
    185          do k=1,klev
    186             do i=1,klon
    187                mfu(i,k)=0.
    188                mfd(i,k)=0.
    189                en_u(i,k)=0.
    190                de_u(i,k)=0.
    191                en_d(i,k)=0.
    192                de_d(i,k)=0.
    193                coefh(i,k)=0.
    194                 t(i,k)=0.
    195                 fm_therm(i,k)=0.
    196                entr_therm(i,k)=0.
    197             enddo
    198          enddo
    199          do i=1,klon
    200             pyv1(i)=0.
    201             pyu1(i)=0.
    202          end do
    203          do k=1,nbsrf
    204              do i=1,klon
    205                pftsol(i,k)=0.
    206                ppsrf(i,k)=0.
    207             enddo
    208          enddo
    209 
    210          dtcum=0.
    211       endif
    212 
    213       do k=1,klev
    214          do i=1,klon
    215             mfu(i,k)=mfu(i,k)+pmfu(i,k)*pdtphys
    216             mfd(i,k)=mfd(i,k)+pmfd(i,k)*pdtphys
    217             en_u(i,k)=en_u(i,k)+pen_u(i,k)*pdtphys
    218             de_u(i,k)=de_u(i,k)+pde_u(i,k)*pdtphys
    219             en_d(i,k)=en_d(i,k)+pen_d(i,k)*pdtphys
    220             de_d(i,k)=de_d(i,k)+pde_d(i,k)*pdtphys
    221             coefh(i,k)=coefh(i,k)+pcoefh_buf(i,k)*pdtphys
    222                 t(i,k)=t(i,k)+pt(i,k)*pdtphys
    223        fm_therm(i,k)=fm_therm(i,k)+pfm_therm(i,k)*pdtphys
    224        entr_therm(i,k)=entr_therm(i,k)+pentr_therm(i,k)*pdtphys
    225          enddo
    226       enddo
    227          do i=1,klon
    228             pyv1(i)=pyv1(i)+yv1(i)*pdtphys
    229             pyu1(i)=pyu1(i)+yu1(i)*pdtphys
    230          end do
    231          do k=1,nbsrf
    232              do i=1,klon
    233                pftsol(i,k)=pftsol(i,k)+ftsol(i,k)*pdtphys
    234                ppsrf(i,k)=ppsrf(i,k)+pctsrf(i,k)*pdtphys
    235             enddo
    236          enddo
    237 
    238       dtcum=dtcum+pdtphys
    239 
    240       IF(mod(iadvtr,istphy).eq.0) THEN
    241 c
    242 c   normalisation par le temps cumule
    243          do k=1,klev
    244             do i=1,klon
    245                mfu(i,k)=mfu(i,k)/dtcum
    246                mfd(i,k)=mfd(i,k)/dtcum
    247                en_u(i,k)=en_u(i,k)/dtcum
    248                de_u(i,k)=de_u(i,k)/dtcum
    249                en_d(i,k)=en_d(i,k)/dtcum
    250                de_d(i,k)=de_d(i,k)/dtcum
    251                coefh(i,k)=coefh(i,k)/dtcum
    252 c Unitel a enlever
    253               t(i,k)=t(i,k)/dtcum       
    254                fm_therm(i,k)=fm_therm(i,k)/dtcum
    255                entr_therm(i,k)=entr_therm(i,k)/dtcum
    256             enddo
    257          enddo
    258          do i=1,klon
    259             pyv1(i)=pyv1(i)/dtcum
    260             pyu1(i)=pyu1(i)/dtcum
    261          end do
    262          do k=1,nbsrf
    263              do i=1,klon
    264                pftsol(i,k)=pftsol(i,k)/dtcum
    265                pftsol1(i) = pftsol(i,1)
    266                pftsol2(i) = pftsol(i,2)
    267                pftsol3(i) = pftsol(i,3)
    268                pftsol4(i) = pftsol(i,4)
    269 
    270                ppsrf(i,k)=ppsrf(i,k)/dtcum
    271                ppsrf1(i) = ppsrf(i,1)
    272                ppsrf2(i) = ppsrf(i,2)
    273                ppsrf3(i) = ppsrf(i,3)
    274                ppsrf4(i) = ppsrf(i,4)
    275 
    276             enddo
    277          enddo
    278 c
    279 c   ecriture des champs
    280 c
    281          irec=irec+1
    282 
    283 ccccc
    284 cym         CALL gr_fi_ecrit(klev,klon,iim,jjm+1, t, zx_tmp_3d)
    285          CALL histwrite_phy(physid,"t",itap,t)
    286 
    287 cym         CALL gr_fi_ecrit(klev,klon,iim,jjm+1, mfu, zx_tmp_3d)
    288       CALL histwrite_phy(physid,"mfu",itap,mfu)
    289 cym     CALL gr_fi_ecrit(klev,klon,iim,jjm+1, mfd, zx_tmp_3d)
    290       CALL histwrite_phy(physid,"mfd",itap,mfd)
    291 cym        CALL gr_fi_ecrit(klev,klon,iim,jjm+1, en_u, zx_tmp_3d)
    292       CALL histwrite_phy(physid,"en_u",itap,en_u)
    293 cym        CALL gr_fi_ecrit(klev,klon,iim,jjm+1, de_u, zx_tmp_3d)
    294       CALL histwrite_phy(physid,"de_u",itap,de_u)
    295 cym        CALL gr_fi_ecrit(klev,klon,iim,jjm+1, en_d, zx_tmp_3d)
    296       CALL histwrite_phy(physid,"en_d",itap,en_d)
    297 cym        CALL gr_fi_ecrit(klev,klon,iim,jjm+1, de_d, zx_tmp_3d)       
    298       CALL histwrite_phy(physid,"de_d",itap,de_d)
    299 cym        CALL gr_fi_ecrit(klev,klon,iim,jjm+1, coefh, zx_tmp_3d)         
    300       CALL histwrite_phy(physid,"coefh",itap,coefh)     
    301 
    302 c ajou...
    303         do k=1,klev
    304            do i=1,klon
    305          fm_therm1(i,k)=fm_therm(i,k)   
    306            enddo
    307         enddo
    308 
    309 cym      CALL gr_fi_ecrit(klev,klon,iim,jjm+1, fm_therm1, zx_tmp_3d)
    310       CALL histwrite_phy(physid,"fm_th",itap,fm_therm1)
    311 c
    312 cym      CALL gr_fi_ecrit(klev,klon,iim,jjm+1, entr_therm, zx_tmp_3d)
    313       CALL histwrite_phy(physid,"en_th",itap,entr_therm)
    314 cccc
    315 cym       CALL gr_fi_ecrit(klev,klon,iim,jjm+1,frac_impa,zx_tmp_3d)
    316         CALL histwrite_phy(physid,"frac_impa",itap,frac_impa)
    317 
    318 cym        CALL gr_fi_ecrit(klev,klon,iim,jjm+1,frac_nucl,zx_tmp_3d)
    319         CALL histwrite_phy(physid,"frac_nucl",itap,frac_nucl)
    320  
    321 cym        CALL gr_fi_ecrit(1, klon,iim,jjm+1, pyu1,zx_tmp_2d)
    322       CALL histwrite_phy(physid,"pyu1",itap,pyu1)
    323        
    324 cym     CALL gr_fi_ecrit(1, klon,iim,jjm+1, pyv1,zx_tmp_2d)
    325       CALL histwrite_phy(physid,"pyv1",itap,pyv1)
    326        
    327 cym     CALL gr_fi_ecrit(1,klon,iim,jjm+1, pftsol1, zx_tmp_2d)
    328       CALL histwrite_phy(physid,"ftsol1",itap,pftsol1)
    329 cym         CALL gr_fi_ecrit(1,klon,iim,jjm+1, pftsol2, zx_tmp_2d)
    330       CALL histwrite_phy(physid,"ftsol2",itap,pftsol2)
    331 cym          CALL gr_fi_ecrit(1,klon,iim,jjm+1, pftsol3, zx_tmp_2d)
    332       CALL histwrite_phy(physid,"ftsol3",itap,pftsol3)
    333 cym         CALL gr_fi_ecrit(1,klon,iim,jjm+1, pftsol4, zx_tmp_2d)
    334       CALL histwrite_phy(physid,"ftsol4",itap,pftsol4)
    335 
    336 cym        CALL gr_fi_ecrit(1,klon,iim,jjm+1, ppsrf1, zx_tmp_2d)
    337       CALL histwrite_phy(physid,"psrf1",itap,ppsrf1)
    338 cym        CALL gr_fi_ecrit(1,klon,iim,jjm+1, ppsrf2, zx_tmp_2d)
    339       CALL histwrite_phy(physid,"psrf2",itap,ppsrf2)
    340 cym        CALL gr_fi_ecrit(1,klon,iim,jjm+1, ppsrf3, zx_tmp_2d)
    341       CALL histwrite_phy(physid,"psrf3",itap,ppsrf3)
    342 cym        CALL gr_fi_ecrit(1,klon,iim,jjm+1, ppsrf4, zx_tmp_2d)
    343       CALL histwrite_phy(physid,"psrf4",itap,ppsrf4)
    344 
    345 c$OMP MASTER
    346       if (ok_sync) call histsync(physid)
    347 c$OMP END MASTER
    348 c     if (ok_sync) call histsync
    349        
    350 c
    351 cAA Test sur la valeur des coefficients de lessivage
    352 c
    353          zmin=1e33
    354          zmax=-1e33
    355          do k=1,klev
    356             do i=1,klon
    357                   zmax=max(zmax,frac_nucl(i,k))
    358                   zmin=min(zmin,frac_nucl(i,k))
    359             enddo
    360          enddo
    361          Print*,'------ coefs de lessivage (min et max) --------'
    362          Print*,'facteur de nucleation ',zmin,zmax
    363          zmin=1e33
    364          zmax=-1e33
    365          do k=1,klev
    366             do i=1,klon
    367                   zmax=max(zmax,frac_impa(i,k))
    368                   zmin=min(zmin,frac_impa(i,k))
    369             enddo
    370          enddo
    371          Print*,'facteur d impaction ',zmin,zmax
    372 
    373       ENDIF
    374 
    375 c   reinitialisation des champs cumules
    376         go to 768
    377       if (mod(iadvtr,istphy).eq.1) then
    378          do k=1,klev
    379             do i=1,klon
    380                mfu(i,k)=0.
    381                mfd(i,k)=0.
    382                en_u(i,k)=0.
    383                de_u(i,k)=0.
    384                en_d(i,k)=0.
    385                de_d(i,k)=0.
    386                coefh(i,k)=0.
    387                t(i,k)=0.
    388                fm_therm(i,k)=0.
    389                entr_therm(i,k)=0.
    390             enddo
    391          enddo
    392          do i=1,klon
    393             pyv1(i)=0.
    394             pyu1(i)=0.
    395          end do
    396          do k=1,nbsrf
    397              do i=1,klon
    398                pftsol(i,k)=0.
    399                ppsrf(i,k)=0.
    400             enddo
    401          enddo
    402 
    403          dtcum=0.
    404       endif
    405 
    406       do k=1,klev
    407          do i=1,klon
    408             mfu(i,k)=mfu(i,k)+pmfu(i,k)*pdtphys
    409             mfd(i,k)=mfd(i,k)+pmfd(i,k)*pdtphys
    410             en_u(i,k)=en_u(i,k)+pen_u(i,k)*pdtphys
    411             de_u(i,k)=de_u(i,k)+pde_u(i,k)*pdtphys
    412             en_d(i,k)=en_d(i,k)+pen_d(i,k)*pdtphys
    413             de_d(i,k)=de_d(i,k)+pde_d(i,k)*pdtphys
    414             coefh(i,k)=coefh(i,k)+pcoefh_buf(i,k)*pdtphys
    415                 t(i,k)=t(i,k)+pt(i,k)*pdtphys
    416        fm_therm(i,k)=fm_therm(i,k)+pfm_therm(i,k)*pdtphys
    417        entr_therm(i,k)=entr_therm(i,k)+pentr_therm(i,k)*pdtphys
    418          enddo
    419       enddo
    420          do i=1,klon
    421             pyv1(i)=pyv1(i)+yv1(i)*pdtphys
    422             pyu1(i)=pyu1(i)+yu1(i)*pdtphys
    423          end do
    424          do k=1,nbsrf
    425              do i=1,klon
    426                pftsol(i,k)=pftsol(i,k)+ftsol(i,k)*pdtphys
    427                ppsrf(i,k)=ppsrf(i,k)+pctsrf(i,k)*pdtphys
    428             enddo
    429          enddo
    430 
    431       dtcum=dtcum+pdtphys
    432 768   continue
    433 
    434       RETURN
    435       END
     1SUBROUTINE phystokenc (nlon,nlev,pdtphys,rlon,rlat, &
     2     pt,pmfu, pmfd, pen_u, pde_u, pen_d, pde_d, &
     3     pfm_therm,pentr_therm, &
     4     cdragh, pcoefh,yu1,yv1,ftsol,pctsrf, &
     5     frac_impa,frac_nucl, &
     6     pphis,paire,dtime,itap, &
     7     psh, pda, pphi, pmp, pupwd, pdnwd)
     8 
     9  USE ioipsl
     10  USE dimphy
     11  USE infotrac, ONLY : nqtot
     12  USE iophy
     13  USE control_mod
     14 
     15  IMPLICIT NONE
     16 
     17!======================================================================
     18! Auteur(s) FH
     19! Objet: Ecriture des variables pour transport offline
     20!
     21!======================================================================
     22  INCLUDE "dimensions.h"
     23  INCLUDE "tracstoke.h"
     24  INCLUDE "indicesol.h"
     25  INCLUDE "iniprint.h"
     26!======================================================================
     27
     28! Arguments:
     29!
     30  REAL,DIMENSION(klon,klev), INTENT(IN)     :: psh   ! humidite specifique
     31  REAL,DIMENSION(klon,klev), INTENT(IN)     :: pda
     32  REAL,DIMENSION(klon,klev,klev), INTENT(IN):: pphi
     33  REAL,DIMENSION(klon,klev), INTENT(IN)     :: pmp
     34  REAL,DIMENSION(klon,klev), INTENT(IN)     :: pupwd ! saturated updraft mass flux
     35  REAL,DIMENSION(klon,klev), INTENT(IN)     :: pdnwd ! saturated downdraft mass flux
     36
     37!   EN ENTREE:
     38!   ==========
     39!
     40!   divers:
     41!   -------
     42!
     43  INTEGER nlon ! nombre de points horizontaux
     44  INTEGER nlev ! nombre de couches verticales
     45  REAL pdtphys ! pas d'integration pour la physique (seconde)
     46  INTEGER itap
     47  INTEGER, SAVE :: physid
     48!$OMP THREADPRIVATE(physid)
     49
     50!   convection:
     51!   -----------
     52!
     53  REAL pmfu(klon,klev)  ! flux de masse dans le panache montant
     54  REAL pmfd(klon,klev)  ! flux de masse dans le panache descendant
     55  REAL pen_u(klon,klev) ! flux entraine dans le panache montant
     56  REAL pde_u(klon,klev) ! flux detraine dans le panache montant
     57  REAL pen_d(klon,klev) ! flux entraine dans le panache descendant
     58  REAL pde_d(klon,klev) ! flux detraine dans le panache descendant
     59  REAL pt(klon,klev)
     60  REAL,ALLOCATABLE,SAVE :: t(:,:)
     61!$OMP THREADPRIVATE(t)
     62!
     63  REAL rlon(klon), rlat(klon), dtime
     64  REAL zx_tmp_3d(iim,jjm+1,klev),zx_tmp_2d(iim,jjm+1)
     65
     66!   Couche limite:
     67!   --------------
     68!
     69  REAL cdragh(klon)          ! cdrag
     70  REAL pcoefh(klon,klev)     ! coeff melange CL
     71  REAL pcoefh_buf(klon,klev) ! coeff melange CL + cdrag
     72  REAL yv1(klon)
     73  REAL yu1(klon),pphis(klon),paire(klon)
     74
     75!   Les Thermiques : (Abderr 25 11 02)
     76!   ---------------
     77  REAL, INTENT(IN) ::  pfm_therm(klon,klev+1)
     78  REAL pentr_therm(klon,klev)
     79 
     80  REAL,ALLOCATABLE,SAVE :: entr_therm(:,:)
     81  REAL,ALLOCATABLE,SAVE :: fm_therm(:,:)
     82!$OMP THREADPRIVATE(entr_therm)
     83!$OMP THREADPRIVATE(fm_therm)
     84!
     85!   Lessivage:
     86!   ----------
     87!
     88  REAL frac_impa(klon,klev)
     89  REAL frac_nucl(klon,klev)
     90!
     91! Arguments necessaires pour les sources et puits de traceur
     92!
     93  REAL ftsol(klon,nbsrf)  ! Temperature du sol (surf)(Kelvin)
     94  REAL pctsrf(klon,nbsrf) ! Pourcentage de sol f(nature du sol)
     95!======================================================================
     96!
     97  INTEGER i, k, kk
     98  REAL,ALLOCATABLE,SAVE :: mfu(:,:)  ! flux de masse dans le panache montant
     99  REAL,ALLOCATABLE,SAVE :: mfd(:,:)  ! flux de masse dans le panache descendant
     100  REAL,ALLOCATABLE,SAVE :: en_u(:,:) ! flux entraine dans le panache montant
     101  REAL,ALLOCATABLE,SAVE :: de_u(:,:) ! flux detraine dans le panache montant
     102  REAL,ALLOCATABLE,SAVE :: en_d(:,:) ! flux entraine dans le panache descendant
     103  REAL,ALLOCATABLE,SAVE :: de_d(:,:) ! flux detraine dans le panache descendant
     104  REAL,ALLOCATABLE,SAVE :: coefh(:,:) ! flux detraine dans le panache descendant
     105 
     106  REAL,ALLOCATABLE,SAVE :: pyu1(:)
     107  REAL,ALLOCATABLE,SAVE :: pyv1(:)
     108  REAL,ALLOCATABLE,SAVE :: pftsol(:,:)
     109  REAL,ALLOCATABLE,SAVE :: ppsrf(:,:)
     110!$OMP THREADPRIVATE(mfu,mfd,en_u,de_u,en_d,de_d,coefh)
     111!$OMP THREADPRIVATE(pyu1,pyv1,pftsol,ppsrf)
     112
     113
     114  REAL,DIMENSION(:,:), ALLOCATABLE,SAVE     :: sh 
     115  REAL,DIMENSION(:,:), ALLOCATABLE,SAVE     :: da
     116  REAL,DIMENSION(:,:,:), ALLOCATABLE,SAVE   :: phi
     117  REAL,DIMENSION(:,:), ALLOCATABLE,SAVE     :: mp
     118  REAL,DIMENSION(:,:), ALLOCATABLE,SAVE     :: upwd
     119  REAL,DIMENSION(:,:), ALLOCATABLE,SAVE     :: dnwd
     120 
     121  REAL, SAVE :: dtcum
     122  INTEGER, SAVE:: iadvtr=0
     123!$OMP THREADPRIVATE(dtcum,iadvtr)
     124  REAL zmin,zmax
     125  LOGICAL ok_sync
     126  CHARACTER(len=12) :: nvar
     127!
     128!======================================================================
     129
     130  iadvtr=iadvtr+1
     131
     132! Dans le meme vecteur on recombine le drag et les coeff d'echange
     133  pcoefh_buf(:,1)      = cdragh(:)
     134  pcoefh_buf(:,2:klev) = pcoefh(:,2:klev)
     135 
     136  ok_sync = .TRUE.
     137
     138! Initialization done only once
     139!======================================================================
     140  IF (iadvtr==1) THEN
     141     ALLOCATE( t(klon,klev))
     142     ALLOCATE( mfu(klon,klev)) 
     143     ALLOCATE( mfd(klon,klev)) 
     144     ALLOCATE( en_u(klon,klev))
     145     ALLOCATE( de_u(klon,klev))
     146     ALLOCATE( en_d(klon,klev))
     147     ALLOCATE( de_d(klon,klev))
     148     ALLOCATE( coefh(klon,klev))
     149     ALLOCATE( entr_therm(klon,klev))
     150     ALLOCATE( fm_therm(klon,klev))
     151     ALLOCATE( pyu1(klon))
     152     ALLOCATE( pyv1(klon))
     153     ALLOCATE( pftsol(klon,nbsrf))
     154     ALLOCATE( ppsrf(klon,nbsrf))
     155     
     156     ALLOCATE(sh(klon,klev))
     157     ALLOCATE(da(klon,klev))
     158     ALLOCATE(phi(klon,klev,klev))
     159     ALLOCATE(mp(klon,klev))
     160     ALLOCATE(upwd(klon,klev))
     161     ALLOCATE(dnwd(klon,klev))
     162
     163     CALL initphysto('phystoke', dtime, dtime*istphy,dtime*istphy,physid)
     164     
     165     ! Write field phis and aire only once
     166     CALL histwrite_phy(physid,"phis",itap,pphis)
     167     CALL histwrite_phy(physid,"aire",itap,paire)
     168     CALL histwrite_phy(physid,"longitudes",itap,rlon)
     169     CALL histwrite_phy(physid,"latitudes",itap,rlat)
     170
     171  END IF
     172 
     173 
     174! Set to zero cumulating fields
     175!======================================================================
     176  IF (MOD(iadvtr,istphy)==1.OR.istphy==1) THEN
     177     WRITE(lunout,*)'reinitialisation des champs cumules a iadvtr=',iadvtr
     178     mfu(:,:)=0.
     179     mfd(:,:)=0.
     180     en_u(:,:)=0.
     181     de_u(:,:)=0.
     182     en_d(:,:)=0.
     183     de_d(:,:)=0.
     184     coefh(:,:)=0.
     185     t(:,:)=0.
     186     fm_therm(:,:)=0.
     187     entr_therm(:,:)=0.
     188     pyv1(:)=0.
     189     pyu1(:)=0.
     190     pftsol(:,:)=0.
     191     ppsrf(:,:)=0.
     192     sh(:,:)=0.
     193     da(:,:)=0.
     194     phi(:,:,:)=0.
     195     mp(:,:)=0.
     196     upwd(:,:)=0.
     197     dnwd(:,:)=0.
     198     dtcum=0.
     199  ENDIF
     200 
     201
     202! Cumulate fields at each time step
     203!======================================================================
     204  DO k=1,klev
     205     DO i=1,klon
     206        mfu(i,k)=mfu(i,k)+pmfu(i,k)*pdtphys
     207        mfd(i,k)=mfd(i,k)+pmfd(i,k)*pdtphys
     208        en_u(i,k)=en_u(i,k)+pen_u(i,k)*pdtphys
     209        de_u(i,k)=de_u(i,k)+pde_u(i,k)*pdtphys
     210        en_d(i,k)=en_d(i,k)+pen_d(i,k)*pdtphys
     211        de_d(i,k)=de_d(i,k)+pde_d(i,k)*pdtphys
     212        coefh(i,k)=coefh(i,k)+pcoefh_buf(i,k)*pdtphys
     213        t(i,k)=t(i,k)+pt(i,k)*pdtphys
     214        fm_therm(i,k)=fm_therm(i,k)+pfm_therm(i,k)*pdtphys
     215        entr_therm(i,k)=entr_therm(i,k)+pentr_therm(i,k)*pdtphys
     216        sh(i,k) = sh(i,k) + psh(i,k)*pdtphys
     217        da(i,k) = da(i,k) + pda(i,k)*pdtphys
     218        mp(i,k) = mp(i,k) + pmp(i,k)*pdtphys
     219        upwd(i,k) = upwd(i,k) + pupwd(i,k)*pdtphys
     220        dnwd(i,k) = dnwd(i,k) + pdnwd(i,k)*pdtphys
     221     ENDDO
     222  ENDDO
     223
     224  DO kk=1,klev
     225     DO k=1,klev
     226        DO i=1,klon
     227           phi(i,k,kk) = phi(i,k,kk) + pphi(i,k,kk)*pdtphys
     228        END DO
     229     END DO
     230  END DO
     231
     232  DO i=1,klon
     233     pyv1(i)=pyv1(i)+yv1(i)*pdtphys
     234     pyu1(i)=pyu1(i)+yu1(i)*pdtphys
     235  END DO
     236  DO k=1,nbsrf
     237     DO i=1,klon
     238        pftsol(i,k)=pftsol(i,k)+ftsol(i,k)*pdtphys
     239        ppsrf(i,k)=ppsrf(i,k)+pctsrf(i,k)*pdtphys
     240     ENDDO
     241  ENDDO
     242 
     243! Add time step to cumulated time
     244  dtcum=dtcum+pdtphys
     245 
     246
     247! Write fields to file, if it is time to do so
     248!======================================================================
     249  IF(MOD(iadvtr,istphy)==0) THEN
     250
     251     ! normalize with time period
     252     DO k=1,klev
     253        DO i=1,klon
     254           mfu(i,k)=mfu(i,k)/dtcum
     255           mfd(i,k)=mfd(i,k)/dtcum
     256           en_u(i,k)=en_u(i,k)/dtcum
     257           de_u(i,k)=de_u(i,k)/dtcum
     258           en_d(i,k)=en_d(i,k)/dtcum
     259           de_d(i,k)=de_d(i,k)/dtcum
     260           coefh(i,k)=coefh(i,k)/dtcum
     261           t(i,k)=t(i,k)/dtcum 
     262           fm_therm(i,k)=fm_therm(i,k)/dtcum
     263           entr_therm(i,k)=entr_therm(i,k)/dtcum
     264           sh(i,k)=sh(i,k)/dtcum
     265           da(i,k)=da(i,k)/dtcum
     266           mp(i,k)=mp(i,k)/dtcum
     267           upwd(i,k)=upwd(i,k)/dtcum
     268           dnwd(i,k)=dnwd(i,k)/dtcum
     269        ENDDO
     270     ENDDO
     271     DO kk=1,klev
     272        DO k=1,klev
     273           DO i=1,klon
     274              phi(i,k,kk) = phi(i,k,kk)/dtcum
     275           END DO
     276        END DO
     277     END DO
     278     DO i=1,klon
     279        pyv1(i)=pyv1(i)/dtcum
     280        pyu1(i)=pyu1(i)/dtcum
     281     END DO
     282     DO k=1,nbsrf
     283        DO i=1,klon
     284           pftsol(i,k)=pftsol(i,k)/dtcum
     285           ppsrf(i,k)=ppsrf(i,k)/dtcum
     286        ENDDO
     287     ENDDO
     288
     289     ! write fields
     290     CALL histwrite_phy(physid,"t",itap,t)
     291     CALL histwrite_phy(physid,"mfu",itap,mfu)
     292     CALL histwrite_phy(physid,"mfd",itap,mfd)
     293     CALL histwrite_phy(physid,"en_u",itap,en_u)
     294     CALL histwrite_phy(physid,"de_u",itap,de_u)
     295     CALL histwrite_phy(physid,"en_d",itap,en_d)
     296     CALL histwrite_phy(physid,"de_d",itap,de_d)
     297     CALL histwrite_phy(physid,"coefh",itap,coefh)     
     298     CALL histwrite_phy(physid,"fm_th",itap,fm_therm)
     299     CALL histwrite_phy(physid,"en_th",itap,entr_therm)
     300     CALL histwrite_phy(physid,"frac_impa",itap,frac_impa)
     301     CALL histwrite_phy(physid,"frac_nucl",itap,frac_nucl)
     302     CALL histwrite_phy(physid,"pyu1",itap,pyu1)
     303     CALL histwrite_phy(physid,"pyv1",itap,pyv1)
     304     CALL histwrite_phy(physid,"ftsol1",itap,pftsol(:,1))
     305     CALL histwrite_phy(physid,"ftsol2",itap,pftsol(:,2))
     306     CALL histwrite_phy(physid,"ftsol3",itap,pftsol(:,3))
     307     CALL histwrite_phy(physid,"ftsol4",itap,pftsol(:,4))
     308     CALL histwrite_phy(physid,"psrf1",itap,ppsrf(:,1))
     309     CALL histwrite_phy(physid,"psrf2",itap,ppsrf(:,2))
     310     CALL histwrite_phy(physid,"psrf3",itap,ppsrf(:,3))
     311     CALL histwrite_phy(physid,"psrf4",itap,ppsrf(:,4))
     312     CALL histwrite_phy(physid,"sh",itap,sh)
     313     CALL histwrite_phy(physid,"da",itap,da)
     314     CALL histwrite_phy(physid,"mp",itap,mp)
     315     CALL histwrite_phy(physid,"upwd",itap,upwd)
     316     CALL histwrite_phy(physid,"dnwd",itap,dnwd)
     317
     318
     319! phi
     320     DO k=1,klev
     321        IF (k<10) THEN
     322           WRITE(nvar,'(i1)') k
     323        ELSE IF (k<100) THEN
     324           WRITE(nvar,'(i2)') k
     325        ELSE
     326           WRITE(nvar,'(i3)') k
     327        END IF
     328        nvar='phi_lev'//trim(nvar)
     329       
     330        CALL histwrite_phy(physid,nvar,itap,phi(:,:,k))
     331     END DO
     332     
     333     ! Syncronize file
     334!$OMP MASTER
     335     IF (ok_sync) CALL histsync(physid)
     336!$OMP END MASTER
     337     
     338     
     339     ! Calculate min and max values for some fields (coefficients de lessivage)
     340     zmin=1e33
     341     zmax=-1e33
     342     DO k=1,klev
     343        DO i=1,klon
     344           zmax=MAX(zmax,frac_nucl(i,k))
     345           zmin=MIN(zmin,frac_nucl(i,k))
     346        ENDDO
     347     ENDDO
     348     WRITE(lunout,*)'------ coefs de lessivage (min et max) --------'
     349     WRITE(lunout,*)'facteur de nucleation ',zmin,zmax
     350     zmin=1e33
     351     zmax=-1e33
     352     DO k=1,klev
     353        DO i=1,klon
     354           zmax=MAX(zmax,frac_impa(i,k))
     355           zmin=MIN(zmin,frac_impa(i,k))
     356        ENDDO
     357     ENDDO
     358     WRITE(lunout,*)'facteur d impaction ',zmin,zmax
     359     
     360  ENDIF ! IF(MOD(iadvtr,istphy)==0)
     361
     362END SUBROUTINE phystokenc
Note: See TracChangeset for help on using the changeset viewer.