!copyright (C) 2001 MSC-RPN COMM %%%RPNPHY%%%
*** S/P PHYEXE1
subroutine phyexe1 (e, d, f, v, 1,29
$ esiz,dsiz,fsiz,vsiz,
$ dt,trnch,kount,task,ni,nk)
*
#include "impnone.cdk"
*
integer esiz,dsiz,fsiz,vsiz,trnch,kount,task,ni,nk
real e(esiz), d(dsiz), f(fsiz), v(vsiz)
real dt
*
*
*Author
* B. Bilodeau Nov 1993 (Routine formerly called "param")
*
*Revisions
* 001 N. Brunet (May 1995) - new surface processes
* 002 M. Gagnon (July 1995) - add reduction parameters to newrad1
* 003 N. Ek (Apr 1995) - Added accumulated surface energy fluxes
* 004 B. Bilodeau (Nov 1995) - Correct bug for MPRECIP
* 005 B. Bilodeau (Jan 1996) - Dynamic allocation of TE
* 006 M. Desgagne and J.-M. Belanger(Nov 1995) -
* Unified physics interface and CLASS
* 007 L. Lefaivre (Nov 95) - Latest version of Mc Farlane (GWDFX95)
* 008 G. Pellerin (Nov 95) - Revised condensation and deep convection
* 009 G. Pellerin (Fev 96) - New options for shallow convection
* 010 B. Dugas (Sep 96) - Correction for stratospheric clouds and
* RADFIX switch added to control usage of radiation fixes
* 011 B. Bilodeau (Sept 1996) - Correct bug in computation of
* tendencies for MPRECIP
* 012 G. Pellerin (Jan 1997) - Execution of whole physics at kount=0.
* Correct extraction of surface fluxes. Change calling sequence
* of vkuocon5.
* 013 M. Desgagne (Apr 1996) - Activate EPLUS (advectke option)
* 014 F. Kong (Dec 1996) - Add new explicit microphysics
* (Ref. Kong & Yau (1996) Atmosphere-Ocean)
* 015 V. Lee (Feb 1996) - Directional roughness length added
* 016 M. Roch (Nov 1997) - Introduce horizontal modulation of sponge
* 017 S. Belair (Spring 1998) - ISBA
* 018 B. Bilodeau (Oct 1998) - Merge phyexe and param4.
* Introduce "entry" bus.
* 019 J. Mailhot (Mar 1999) - Changes for new SURFACE interface
* B. Bilodeau
* 020 S. Belair (January 2000) - Accumulators for ISBA
* 021 B. Bilodeau (Nov 2000) - New comdeck phybus.cdk
* 022 B. Bilodeau (January 2001) - Dynamic memory allocation
* revisited
* 023 B. Dugas (July 2000) - Replace CLIMPHS by CLIMPHS2.
* Add MOYHR cloud average ccnm as well
* as the ttmin, ttmax temperatures.
* 024 J. Mailhot (May 2000) - Changes to add MOISTKE option (ifluvert=3)
* 025 J. Mailhot (Jun 2000) - Correct bug in interface for MIXED-PHASE
* 026 A. Erfani and B. Bilodeau (Oct 2001) - Added the precipitation
* partitioning code developed by A. Methot
* 027 D. Talbot (Oct 2001) - Call to gwd4 (blocking)
* 028 B. Dugas (Nov 2001) - Add the suaf and svaf accumulators
* 029 B. Bilodeau (Mar 2002) - QDIFV tendency = 0 if wet=.false.
* 030 S. Laroche (Apr 2002) - Call simplified physics options as
* suggested by B. Bilodeau and M. Desgagne
* 031 B. Bilodeau (Jun 2002) - Copy level NK of HUMOINS and TMOINS
* into HUCOND and TCOND for Kong-Yau
* 032 A.-M. Leduc, S. Belair and B. Bilodeau (Feb 2002) -
* Add KFC implicit condensate to IWC and LWC
* 033 S.Belair, A-M. Leduc (Nov 2002) - add averaged tendencies for kfc
* add zsqcem
* 034 A-M. Leduc (Nov 2002) - add switch for shallow convection ishlcvt(2)
* 035 B. Bilodeau (Feb 2003) - call to calcdiag and extdiag
* 036 J. Mailhot (Feb 2003) - Changes to the MOISTKE and ADVECTKE options
* 037 S. Belair (Apr 2003) - Changes to the lwc, iwc, and cloud
* fractions (better treatment of shallow
* and convective components
* 038 A. PLante (Sep 2003) - Update doc for key ipcptype (P_cond_pcptype_s)
*
* 039 B. Bilodeau and L. Spacek (Dec 2003) - Move ttmin and ttmax to calcdiag
* 040 B. Bilodeau and L. Spacek (Dec 2003) - Move ttmin and ttmax to calcdiag
* 041 B. Bilodeau (Feb 2004) - Move call to fisimp3 from surface to phyexe1
*
*Object
* this is the main interface subroutine for the
* CMC/RPN unified physics
*
*Arguments
*
* - Input -
* E entry input field
* D dynamics input field
*
* - Input/Output -
* F historic variables for the physics
*
* - Output -
* V physics tendencies and other output fields from the physics
*
* - Input -
* ESIZ dimension of e
* DSIZ dimension of d
* FSIZ dimension of f
* VSIZ dimension of v
* DT timestep (sec.)
* TRNCH slice number
* KOUNT timestep number
* ICPU cpu number executing slice "trnch"
* N horizontal running length
* NK vertical dimension
*
*Notes
* PHYEXE is called by all the models that use the CMC/RPN
* common physics library. It returns tendencies to the
* dynamics.
*
*IMPLICITES
*
#include "phy_macros_f.h"
#include "phybus.cdk"
#include "nbvarsurf.cdk"
#include "dimsurf.cdk"
#include "workspc.cdk"
#include "options.cdk"
#include "consphy.cdk"
*
*MODULES
*
**
*
integer ikk
integer icpu
integer i,j,k,nsups
integer ierget
real cdt1, cdt2, rcdt1
real heurser
*
*
real*8 ppjour,locals,demipa
integer nik,nni,nnik
integer maxadj
parameter (maxadj=20)
integer itadj(maxadj)
*
************************************************************************
* AUTOMATIC ARRAYS
************************************************************************
*
*
AUTOMATIC ( EPLUS0 , REAL , (NI,NK ) )
AUTOMATIC ( HUPLUS0 , REAL , (NI,NK ) )
AUTOMATIC ( QCPLUS0 , REAL , (NI,NK ) )
AUTOMATIC ( QIPLUS0 , REAL , (NI,NK ) )
AUTOMATIC ( QGPLUS0 , REAL , (NI,NK ) )
AUTOMATIC ( QRPLUS0 , REAL , (NI,NK ) )
AUTOMATIC ( TPLUS0 , REAL , (NI,NK ) )
AUTOMATIC ( UPLUS0 , REAL , (NI,NK ) )
AUTOMATIC ( VPLUS0 , REAL , (NI,NK ) )
AUTOMATIC ( WPLUS0 , REAL , (NI,NK ) )
*
AUTOMATIC ( GZMOINS , REAL , (NI,NK ) )
AUTOMATIC ( QCDIFV , REAL , (NI,NK ) )
AUTOMATIC ( QE , REAL , (NI,NK ) )
AUTOMATIC ( SELOC , REAL , (NI,NK ) )
AUTOMATIC ( TRAV1D , REAL , (NI ) )
AUTOMATIC ( TRAV2D , REAL , (NI,NK,4) )
AUTOMATIC ( TVIRT , REAL , (NI,NK ) )
AUTOMATIC ( FICEBL , REAL , (NI,NK ) )
AUTOMATIC ( WORK , REAL , (ESPWORK) )
AUTOMATIC ( QCBLIC , REAL , (NI,NK ) )
*
************************************************************************
*
*
external lin_phyexe1,serxst,vkuocon6,mprecip4
external difuvd9,serget
external gwd5,cldwin
external newrad3, oldrad3
external mzoniv, mzonxst
external difver6,mfotvt,adilwc,climphs2
external ficemxp,shallconv2
*
*
*----------------------------------------------------------------
*
* options de la physique
* ----------------------
*
* iconvec : convection
*
* 0 = 'nil'
* 1 = 'sec'
* 2 = 'manabe'
* 3 = 'oldkuo'
* 4 = 'newkuo'
* 5 = 'fcp'
* 6 = 'kfc'
* 7 = 'kuostd'
* 8 = 'kuosym'
* 9 = 'kuosun'
* 10 = 'ras'
* 11 = 'fcpkuo'
* 12 = 'fcpkuo2'
* 13 = 'bm'
*
* ifluvert : couche limite
*
* 0 = 'nil'
* 1 = 'physimp'
* 2 = 'clef'
* 3 = 'moistke'
*
* igwdrag : gravity wave drag
*
* 0 = 'nil'
* 1 = 'gwd86'
* 2 = 'gwd95'
*
* iradia : radiation
*
* 0 = 'nil'
* 1 = 'oldrad'
* 2 = 'newrad'
*
* ischmsol : svat (transferts energetiques entre le sol,
* la vegetation et l'atmosphere)
*
* 1 = 'fcrest'
* 2 = 'class'
* 3 = 'isba'
*
* ishlcvt : convection restreinte (2 valeurs)
*
* 1) 0 = 'nil'
* 1 = 'geleyn'
* 2 = 'conres'
* 3 = 'shalow'
* 4 = 'shalodqc'
* 5 = 'bmshal'
*
* 2) 0 = 'nil'
* 1 = 'ktrsnt'
* 2 = 'salzen'
*
*
* istcond : condensation stratiforme
*
* 0 = 'nil'
* 1 = 'conds'
* 2 = 'oldsund'
* 3 = 'newsund'
* 4 = 'consun'
* 5 = 'exc' !tremblay
* 6 = 'excr1' !ancien exmo
* 7 = 'excr2'
* 8 = 'excri'
* 9 = 'excrig'
*
* ilongmel : longueur de melange
*
* 0 = 'blac62'
* 1 = 'boujo'
*
* ikfcpcp : conservation de l'eau dans kfcp
*
* 0 = 'ori' ! code original
* 1 = 'conspcpn'
*
* ipcptype : Diagnostique de type de precipitations
*
* 0 = 'nil'
* 1 = 'bourge' ! Methode de Bourgouin etendue
*
*-------------------------------------------------------------
*
*
real eplus
pointer (paeplus, eplus(ni,nk))
*
integer ik
* fonction-formule pour faciliter le calcul des indices
ik(i,k) = (k-1)*ni + i -1
*
*
************************************************************************
* appel a la physique simplifiee si desire *
* ---------------------------------------- *
************************************************************************
if(lin_kph.eq.1) then
*
call lin_phyexe1 (e, d, f, v,
$ esiz,dsiz,fsiz,vsiz,
$ dt,trnch,kount,task,ni,nk)
*
return
endif
*
*
************************************************************************
* preparatifs *
* ----------- *
************************************************************************
*
if (.not.advectke) then
paeplus = loc(f(en))
else
paeplus = loc(d(enplus))
endif
*
icpu = task
*
*
************************************************************************
* garder une copie des champs du temps plus *
* ----------------------------------------- *
************************************************************************
*
do 5 k=1,nk
*vdir nodep
do 5 i = 1,ni
huplus0(i,k) = d(huplus+ik(i,k))
tplus0(i,k) = d( tplus+ik(i,k))
uplus0(i,k) = d( uplus+ik(i,k))
vplus0(i,k) = d( vplus+ik(i,k))
qcplus0(i,k) = d(qcplus+ik(i,k))
qrplus0(i,k) = d(qrplus+ik(i,k))
5 continue
if (advectke) then
do k=1,nk
*vdir nodep
do i = 1,ni
eplus0(i,k) = d(enplus+ik(i,k))
end do
end do
endif
if (diffuw) then
do k=1,nk
*vdir nodep
do i = 1,ni
wplus0(i,k) = d(omegap+ik(i,k))
end do
end do
endif
if(istcond.ge.8) then
do k=1,nk
*vdir nodep
do i=1,ni
qiplus0(i,k) = d(qiplus+ik(i,k))
end do
end do
endif
if(istcond.ge.9) then
do k=1,nk
*vdir nodep
do i=1,ni
qgplus0(i,k) = d(qgplus+ik(i,k))
end do
end do
endif
*
*
*
************************************************************************
* constantes derivees du pas de temps
* ----------------------------------- *
************************************************************************
*
cdt1 = factdt * dt
cdt2 = factdt * dt * facdifv
rcdt1 = 1./cdt1
*
if (climat) then
*
* sommes-nous pres du milieu du jour ?
*
ppjour = 86400 / dble( delt )
locals = kount / ppjour
locals = locals - int( locals )
demipa = 1./(2. * ppjour)
*
if (0.5-demipa .lt. locals
+ .and. locals .le. 0.5+demipa)
+ call climphs2(f,fsiz,kount,ni)
*
end if
*
*
*
************************************************************************
* initialisations *
* --------------- *
************************************************************************
*
* calcul des niveaux intermediaires
call difuvd9(seloc,.true.,d(sigm),ni,nk,nk-1)
*
call inichamp1 (e, esiz, f, fsiz,
$ v, vsiz, d, dsiz,
$ qcplus0, qcdifv,
$ trav2d, seloc, kount, trnch,
$ dt, cdt1, ni, nk)
*
* z0 directionnel
if (z0dir) then
* calcul de z0 avec z1,z2,z3,z4 et umoins,vmoins
call calcz0(f(mg),f(z0),f(z1),f(z2),f(z3),f(z4),
$ d(umoins+ik(1,nk-1)),
$ d(vmoins+ik(1,nk-1)), ni)
endif
*
************************************************************************
* extraction des tendances de la dynamique *
* ---------------------------------------- *
************************************************************************
*
call serget ( 'HEURE' , heurser , 1 , ierget )
*
call sersetm( 'KA', trnch, nk )
call mzoniv ( trnch, nk )
*
do k=1,nk
do i = 1,ni
trav2d(i,k,2) = (d(tplus +ik(i,k)) -
$ d(tmoins+ik(i,k))) * rcdt1
trav2d(i,k,3) = (d(huplus+ik(i,k)) -
$ d(humoins+ik(i,k))) * rcdt1
if(istcond.ge.2) then
trav2d(i,k,4) = (d(qcplus+ik(i,k)) -
$ d(qcmoins+ik(i,k))) * rcdt1
endif
end do
end do
*
call serxst (trav2d(1,1,2),'XT',trnch,ni,0., 1., -1 )
call mzonxst(trav2d(1,1,2),'XT',trnch,ni,heurser,pmoins,-2, icpu)
call serxst (trav2d(1,1,3),'XQ',trnch,ni,0., 1., -1 )
call mzonxst(trav2d(1,1,3),'XQ',trnch,ni,heurser,pmoins,-2, icpu)
*
if(istcond.ge.2) then
call serxst (trav2d(1,1,4),'XL',trnch,ni,0., 1., -1 )
call mzonxst(trav2d(1,1,4),'XL',trnch,ni,heurser,pmoins,-2, icpu)
endif
*
call sersetm( 'KA', trnch, nk-1 )
call mzoniv ( trnch, -(nk-1) )
*
if (istcond.eq.1) then
* trav2d est calcule par inichamp1
call serxst (trav2d ,'LW',trnch,ni,0., 1., -1 )
call mzonxst(trav2d ,'LW',trnch,ni,heurser,pmoins,-2, icpu)
endif
*
*
************************************************************************
* calculs radiatifs *
* ----------------- *
************************************************************************
*
if (iradia.ge.1) then
*
if (iradia.eq.2) then
*
* nouveau scheme de radiation
*
call newrad3 (f, fsiz, v, vsiz, work, espwork,
+ d(tmoins), d(humoins),
+ d(pmoins), d(sigm), delt, kount,
+ trnch, ni, ni, nk-1, icpu, icpu,
+ nk, radnivl(1)-1, radnivl(1), radnivl(2))
*
else if (iradia.eq.1) then
*
* ancien scheme de radiation
*
call oldrad3 (f, fsiz, v, vsiz, work, espwork,
+ d(tmoins), d(humoins),
+ d(pmoins), d(sigm), delt, satuco,
+ radfix, kount, date, kntrad, trnch,
+ ni, ni, nk-1, dbgmem, icpu)
*
endif
*
*
* tendances de la radiation
do i = 1,ni*(nk-1)
v(trad+i-1) = f(ti+i-1) + f(t2+i-1)
end do
*
* tendances nulles au niveau diagnostique
do i = 1,ni
v(trad+ik(i,nk)) = 0.0
end do
*
else if (iradia.eq.0) then
*
* pas de radiation, tendances nulles
*
do i = 1,ni*nk
v(trad+i-1) = 0.0
end do
*
endif
*
*
************************************************************************
* calculs preliminaires a la diffusion verticale *
* ---------------------------------------------- *
************************************************************************
*
* calcul de la temperature virtuelle (tve),
* de l'humidite specifique (qe) et des hauteurs
* geopotentielles (ze) aux niveaux decales
*
do k=1,nk-2
*VDIR NODEP
do i=1,ni
v(tve+ik(i,k)) = ( d(tmoins+ik(i,k)) +
$ d(tmoins+ik(i,k+1)))/2.
qe(i,k ) = ( d(humoins+ik(i,k )) +
$ d(humoins+ik(i,k+1)))/2.
end do
end do
*
do i=1,ni
v(tve+ik(i,nk-1)) = d( tmoins+ik(i,nk-1))
qe( i,nk-1) = d(humoins+ik(i,nk-1))
enddo
*
call mfotvt(v(tve),v(tve),qe,ni,nk-1,ni)
*
do i=1,ni
v(ze+ik(i,nk-1)) = 0.
end do
*
call integ2 ( v(ze), v(tve), -rgasd/grav, seloc,
$ trav2d(1,1,1),trav2d(1,1,2),trav2d(1,1,3),
$ ni, ni, ni, nk-1, .true. )
*
* calcul de la densite de l'air (pour AURAMS)
do k=1,nk
do i=1,ni
v(rhod+ik(i,k)) = d(sigm+ik(i,k))*d(pmoins+ik(i,1))/
$ (rgasd*d(tmoins+ik(i,k))*(1.0+delta*d(humoins+ik(i,k))))
end do
end do
*
*
************************************************************************
* physique simplifiee *
* ------------------- *
************************************************************************
*
if (ifluvert.eq.1) then
*
call fisimp3 ( d,dsiz,f,fsiz,v,vsiz,
$ seloc,trav2d,
$ ni,nk-1,trav1d,kount)
*
*
************************************************************************
* processus de surface *
* -------------------- *
************************************************************************
*
else if (ifluvert.ge.2) then
*
call surface ( d, dsiz, f, fsiz, v, vsiz,
$ work, espwork, seloc, trnch,
$ kount, dt, ni, ni, nk, icpu )
*
*
************************************************************************
*
* energie cinetique turbulente, operateurs de diffusion *
* ----------------------------------------------------- *
* *
* et hauteur de la couche limite stable ou instable *
* ------------------------------------------------- *
* *
************************************************************************
*
*
call turbul ( d, dsiz, f, fsiz, v, vsiz,
$ work, espwork,
$ eplus, qe, seloc,
$ kount, trnch, ni, ni, nk-1, icpu )
*
*
* calcul des tendances de TKE
* ---------------------------
*
if (advectke) then
do k = 1,nk-1
do i = 1,ni
v(enphytd+ik(i,k)) = (eplus(i,k)-eplus0(i,k))*rcdt1
end do
end do
do i = 1,ni
v(enphytd+ik(i,nk)) = 0.0
end do
endif
*
*
endif
*
*
************************************************************************
* gravity wave drag *
* ----------------- *
************************************************************************
*
if(igwdrag.eq.1 .or. igwdrag.eq.2 ) then
*
* calcul de la temperature virtuelle au temps moins
call mfotvt(tvirt,d(tplus),d(huplus),ni,nk,ni)
*
CALL GWD4 (F, FSIZ, D(UPLUS), D(VPLUS), TVIRT, D(PPLUS), D(SIGM),
+ V(UGWD), V(VGWD), CDT1, KOUNT, TRNCH, NI, NI, NK-1,
+ ICPU )
*
* note : les tendances dues au gravity wave drag sont
* appliquees dans les sous-programmes gwdflx2 et
* gwdfx95
*
* tendances dues au gravity wave drag mises a zero
* au niveau diagnostique
do i=1,ni
v(ugwd+ik(i,nk)) = 0.
v(vgwd+ik(i,nk)) = 0.
end do
*
endif
*
* application des tendances de la radiation
* -----------------------------------------
*
if (iradia.ge.1) then
do k = 1,nk-1
do i = 1,ni
d(tplus+ik(i,k)) = d(tplus+ik(i,k)) + cdt1*v(trad+ik(i,k))
end do
end do
endif
*
*
************************************************************************
* diffusion verticale *
* ------------------- *
************************************************************************
*
if (ifluvert.ge.2.or.(ifluvert.eq.1.and.dmom)) then
*
* calcul des tendances de la diffusion au niveau diagnostique
* (dont les series temporelles sont extraites dans difver6)
*
*VDIR NODEP
do i=1,ni
qcdifv(i,nk) = 0.0
*
* tendances d'humidite specifique nulles si le modele est sec
if (wet) then
v( qdifv+ik(i,nk)) = (f(qdiag+i-1)-d(huplus+ik(i,nk)))*rcdt1
endif
*
v( tdifv+ik(i,nk)) = (f(tdiag+i-1)-d( tplus+ik(i,nk)))*rcdt1
v( udifv+ik(i,nk)) = (f(udiag+i-1)-d( uplus+ik(i,nk)))*rcdt1
v( vdifv+ik(i,nk)) = (f(vdiag+i-1)-d( vplus+ik(i,nk)))*rcdt1
end do
*
else if (ifluvert.eq.1.and..not.dmom) then
*
do i=1,ni
qcdifv(i,nk) = 0.0
v( qdifv+ik(i,nk)) = 0.0
v( tdifv+ik(i,nk)) = 0.0
v( udifv+ik(i,nk)) = 0.0
v( vdifv+ik(i,nk)) = 0.0
end do
if (diffuw) then
do i=1,ni
v(wdifv+ik(i,nk)) = 0.0
end do
endif
*
endif
*
call difver6 (d, dsiz, f, fsiz, v, vsiz,
$ work, espwork, qcdifv, seloc,
$ cdt1, kount, trnch, ni, nk-1,
$ icpu )
*
*
* application des tendances de la diffusion
* -----------------------------------------
*
if (ifluvert. ge. 1) then
do k = 1,nk-1
*VDIR NODEP
do i = 1,ni
d(huplus+ik(i,k)) = d(huplus+ik(i,k)) +
$ cdt1 * v(qdifv+ik(i,k))
d( tplus+ik(i,k)) = d( tplus+ik(i,k)) +
$ cdt1 * v(tdifv+ik(i,k))
d( uplus+ik(i,k)) = d( uplus+ik(i,k)) +
$ cdt1 * v(udifv+ik(i,k))
d( vplus+ik(i,k)) = d( vplus+ik(i,k)) +
$ cdt1 * v(vdifv+ik(i,k))
d(qcplus+ik(i,k)) = d(qcplus+ik(i,k)) +
$ cdt1 * qcdifv(i,k)
end do
end do
if (diffuw) then
do k = 1,nk-1
do i = 1,ni
d(omegap+ik(i,k))= d(omegap+ik(i,k)) +
$ cdt1 * v(wdifv+ik(i,k))
end do
end do
endif
endif
*
if (ifluvert. eq. 3) then
* BL ice fraction for later use (in cloud water section)
call ficemxp (ficebl, trav2d(1,1,1), trav2d(1,1,2),
$ d(tplus), ni, ni, nk-1)
endif
*
************************************************************************
* processus de convection/condensation *
* ------------------------------------ *
************************************************************************
*
*
*
* calcul de la temperature virtuelle au temps moins
call mfotvt(tvirt,d(tmoins),d(humoins),ni,nk,ni)
*
* calcul du geopotentiel au temps moins
do i=1,ni
gzmoins(i,nk) = 0.0
end do
call integ2 ( gzmoins, tvirt, -rgasd, d(sigm),
$ trav2d(1,1,1), trav2d(1,1,2), trav2d(1,1,3),
$ ni, ni, ni, nk, .true. )
*
*
*-----------------------------------------------------
* shallow convection calculation: ouput v(tshal) et v(hushal)
*
*
if (ishlcvt(2).gt.0) then
CALL shallconv3( d, dsiz, f, fsiz, v, vsiz, kount, trnch,
1 cdt1, ni, nk )
endif
*
****************************************************************************
do k = 1,nk-1
*VDIR NODEP
do i = 1,ni
d(huplus+ik(i,k)) = d(huplus+ik(i,k)) +
$ cdt1 * v(hushal+ik(i,k))
d( tplus+ik(i,k)) = d( tplus+ik(i,k)) +
$ cdt1 * v(tshal+ik(i,k))
end do
end do
*----------------------------------------------------------
if ( iconvec.eq.1 .or. iconvec.ge.3 .or. istcond.ge.1 ) then
*
* transvidage
*VDIR NODEP
do i=1,ni*nk
v(hucond+i-1) = d(humoins+i-1)
v( tcond+i-1) = d( tmoins+i-1)
end do
*
*
call vkuocon6 (d, dsiz, f, fsiz, v, vsiz,
$ work, espwork, gzmoins, seloc,
$ dt, ni, ni, nk-1,
$ kount, trnch, icpu)
*
if ( iconvec.eq.4 .or. istcond.eq.3 .or. istcond.eq.4 ) then
*
* transvider l'avant-dernier niveau dans le niveau diagnostique
*vdir nodep
do i=1,ni
*
v( rnflx+ik(i,nk)) = v( rnflx+ik(i,nk-1))
v(snoflx+ik(i,nk)) = v(snoflx+ik(i,nk-1))
*
end do
*
endif
*
*
else if ( iconvec .eq. 2 ) then
*
* convection selon schema de manabe
* ---------------------------------
*
*VDIR NODEP
do i=1,ni*(nk-1)
v(hucond+i-1) = d(huplus+i-1)
v( tcond+i-1) = d( tplus+i-1)
end do
*
call mprecip4 (v(tcond), v(hucond), f(tls),
+ d(omegap), d(pplus), v(kcl),
+ satuco, d(sigm),
+ cdt1, itadj, maxadj, nsups,
+ ni, nk, nk-1, ni)
*
* calcul des tendances de la convection
*VDIR NODEP
do i=1,ni*(nk-1)
v(hucond+i-1) = (v(hucond+i-1) -
+ d(huplus+i-1)) * rcdt1
v( tcond+i-1) = (v( tcond+i-1) -
+ d( tplus+i-1)) * rcdt1
end do
*
endif
*
*
*
* tendances de la convection/condensation nulles au niveau diagnostique
* --------------------------------------------------------
*
do i=1,ni
*
v( tcond+ik(i,nk)) = 0.0
v(hucond+ik(i,nk)) = 0.0
v(qccond+ik(i,nk)) = 0.0
f(fice +ik(i,nk)) = 0.0
v(tshal +ik(i,nk)) = 0.0
v(hushal+ik(i,nk)) = 0.0
*
end do
*
*
*
* application des tendances de la convection/condensation
* ------------------------------------------
*
if (iconvec.gt.0 .or. istcond.gt.0) then
*
*VDIR NODEP
do i=0,ni*(nk-1)-1
d(huplus+i) = d(huplus+i) + cdt1 * v(hucond+i)
d( tplus+i) = d( tplus+i) + cdt1 * v( tcond+i)
d(qcplus+i) = d(qcplus+i) + cdt1 * v(qccond+i)
end do
*
*
*
* add shallow convection tendencies to convection/condensation tendencies
* ------------------------------------------
*
do i=0,ni*(nk-1) -1
v(hucond+i) = v(hucond+i) + v(hushal+i)
v(tcond+i) = v(tcond+i) + v(tshal +i)
end do
*
*
*
* combine diffusion tendency for cloud water to condensation tendency
* -------------------------------------------------------------------
*
do k=1,nk-1
do i=1,ni
*
v(qccond+ik(i,k)) = v(qccond+ik(i,k)) + qcdifv(i,k)
*
* qccond must be greater or equal to zero
if (kount.eq.0) then
v(qccond+ik(i,k)) = max( v(qccond+ik(i,k)), 0. )
else if (kount.gt.0) then
v(qccond+ik(i,k)) = max( v(qccond+ik(i,k)),
+ -max(0.0 ,qcplus0(i,k)) /cdt1 )
endif
*
end do
end do
*
*
*
* eau nuageuse
* ------------
*
if (istcond.ge.2) then
*
* transvidage de qcplus dans une variable permanente qui
* sera utilisee dans la radiation au pas de temps suivant
do i=1,ni*(nk-1)
f(lwc+i-1) = d(qcplus+i-1)
enddo
*
*
*
* If we have the CONSUN scheme, the cloud
* water from MoisTKE has priority over the
* cloud water from the grid-scale scheme
* In this case, everything goes into the
* LWC (for liquid cloud water) variable
*
if (istcond.eq.4) then
*VDIR NODEP
do i=1,ni*(nk-1)
if (ifluvert.eq.3. .and. f(qcbl+i-1).gt.d(qcplus+i-1)) then
f(lwc+i-1) = f(qcbl+i-1)
v(ccs+i-1) = f(fn+i-1)
endif
end do
*
endif
*
*
*
* Same for MixPhase ...
* BUT the partition of cloud water is also done
* using fice and ficebl
*
if (istcond.eq.5) then
*VDIR NODEP
do i=1,ni*(nk-1)
f(lwc+i-1) = d(qcplus+i-1) * (1.0 - f(fice+i-1) )
f(iwc+i-1) = d(qcplus+i-1) * f(fice+i-1)
*
if (ifluvert.eq.3. .and. f(qcbl+i-1).gt.d(qcplus+i-1)) then
f(lwc+i-1) = f(qcbl+i-1) * (1.0 - ficebl(i,1) )
f(iwc+i-1) = f(qcbl+i-1) * ficebl(i,1)
v(ccs+i-1) = f(fn+i-1)
endif
end do
*
endif
*
*
*
* For the Kong and Yau scheme, MoisTKE has
* priority over the grid-scale condensation
* scheme
*
if (istcond.ge.6) then ! schemas "EXCR1"
* "EXCR2"
* "EXCRI"
* "EXCRG"
*VDIR NODEP
do i=1,ni*(nk-1)
* transvidage de qcplus (eau nuageuse) et qrplus (pluie)
* dans une variable permanente qui sera utilisee dans la
* radiation au pas de temps suivant
f(lwc+i-1) = d(qcplus+i-1) + d(qrplus+i-1)
end do
endif
*
if (istcond.eq.8) then ! schema "EXCRI"
do i=1,ni*(nk-1)
* transvidage de qiplus (glace)
f(iwc+i-1) = d(qiplus+i-1)
end do
endif
*
if (istcond.eq.9) then ! schema "EXCRG"
do i=1,ni*(nk-1)
* transvidage de qiplus (glace) et
* qgplus (neige roulee ou "graupel")
f(iwc+i-1) = d(qiplus+i-1) + d(qgplus+i-1)
*
if (ifluvert.eq.3. .and. f(qcbl+i-1).gt.
+ f(lwc+i-1)+f(iwc+i-1)) then
*
f(lwc+i-1) = f(qcbl+i-1) * (1.0 - ficebl(i,1) )
f(iwc+i-1) = f(qcbl+i-1) * ficebl(i,1)
v(ccs+i-1) = f(fn+i-1)
*
endif
*
end do
endif
*
if(kount.eq.1) then
do i=1,ni*(nk-1)
* necessaire jusqu'a ce que les zonaux a t=0
* soient corriges
d(qcplus+i-1)=qcplus0(i,1)
enddo
endif
*
endif
*
*
*
* Add the cloud water (liquid and solid) coming from
* shallow and deep cumulus clouds (only for the
* Kuo Transient and Kain-Fritsch schemes).
* Note that no conditions are used for these
* calculations ... rliqkfc, ricekfc, and qkt
* are zero is these schemes are not used.
* Also note that rliqkfc and ricekfc are IN-CLOUD values.
*
* If the MixPhase or Kong-Yau schemes are used, it is
* better to partition cloud water into liquid and
* solid phases (possible for the Kain-Fritsch scheme)
* For the other grid-scale schemes (e.g., Sundqvist)
* all the cloud water is put in LWC (and will be
* partition later in the radiation subroutines)
*
if (istcond.GE.5) then
*
do i=1,ni*(nk-1)
f(lwc+i-1) = f(lwc+i-1) + f(rliqkfc+i-1)*f(cck+i-1)
1 + v(qktl+i-1) *v(ckt+i-1)
f(iwc+i-1) = f(iwc+i-1) + f(ricekfc+i-1)*f(cck+i-1)
1 + v(qkti+i-1) *v(ckt+i-1)
end do
*
else
*
do i=1,ni*(nk-1)
f(lwc+i-1) = f(lwc+i-1)
1 + ( f(rliqkfc+i-1)+f(ricekfc+i-1) ) * f(cck+i-1)
1 + ( v(qktl+i-1) +v(qkti+i-1) ) * v(ckt+i-1)
end do
*
end if
*
*
*
*
* Combine explicit and implicit clouds.
* using the random overlap approximation
* CCS is for grid-scale condensation
* (but can also include the PBL
* clouds of MoisTKE)
* CCK is for deep convection clouds
* CKT is for the shallow convection clouds
*
do i=1,ni*(nk-1)
F(CCN+I-1) = MIN( MAX(
1 1. - (1.-V(CCS+I-1))*(1.-F(CCK+I-1))*(1.-V(CKT+I-1)) ,
1 0. ) ,
1 1.)
end do
*
endif
*
*
************************************************************************
* Calcul de moyennes et d'accumulateurs *
* ------------------------------------- *
************************************************************************
*
call calcdiag (d,f,v,dsiz,fsiz,vsiz,dt,trnch,kount,ni,nk)
*
*
************************************************************************
* extraction de diagnostics *
* ------------------------- *
************************************************************************
*
call extdiag (d,f,v,dsiz,fsiz,vsiz,trnch,icpu,ni,nk)
*
*
************************************************************************
* Remettre les champs du temps plus a leur valeur initiale *
* -------------------------------------------------------- *
************************************************************************
*
do 6 k=1,nk
*VDIR NODEP
do 6 i = 1,ni
d(huplus+ik(i,k)) = huplus0(i,k)
d( tplus+ik(i,k)) = tplus0(i,k)
d( uplus+ik(i,k)) = uplus0(i,k)
d( vplus+ik(i,k)) = vplus0(i,k)
d(qcplus+ik(i,k)) = qcplus0(i,k)
d(qrplus+ik(i,k)) = qrplus0(i,k)
6 continue
*
if(istcond.ge.8) then
do k=1,nk
*VDIR NODEP
do i=1,ni
d(qiplus+ik(i,k)) = qiplus0(i,k)
end do
end do
endif
if(istcond.ge.9) then
do k=1,nk
do i=1,ni
d(qgplus+ik(i,k)) = qgplus0(i,k)
end do
end do
endif
if (advectke) then
do k=1,nk
do i = 1,ni
d(enplus+ik(i,k)) = eplus0(i,k)
end do
end do
endif
if (diffuw) then
do k=1,nk
do i = 1,ni
d(omegap+ik(i,k)) = wplus0(i,k)
end do
end do
endif
*
return
end