source: LMDZ6/branches/Portage_acc/libf/phylmd/thermcell_qsat.F90 @ 4584

Last change on this file since 4584 was 4435, checked in by Laurent Fairhead, 19 months ago

Routines adapted to GPU. Speed-up is some 200 w.r.t. 1 CPU as measured by SYSTEM_CLOCK calls
around the call to the routine, CPU takes .14 (seconds?), GPU around 6E-004 (s?) at a resolution of 144x142x95.
Comparison of the netcdf result files output with the replay method yields no difference between the
two runs
  • Property copyright set to
    Name of program: LMDZ
    Creation date: 1984
    Version: LMDZ5
    License: CeCILL version 2
    Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
    See the license file in the root directory
File size: 3.1 KB
Line 
1subroutine thermcell_qsat(klon,active,pplev,ztemp,zqta,zqsat)
2  USE thermcell_ini_mod, ONLY : RLvCp,RETV, RCPD, RLVTT, RTT
3  implicit none
4
5!!#include "YOMCST.h"
6#include "YOETHF.h"
7#include "FCTTRE.h"
8
9
10!====================================================================
11! DECLARATIONS
12!====================================================================
13
14! Arguments
15INTEGER, INTENT(IN) :: klon
16!REAL zpspsk(klon),pplev(klon)
17!REAL ztemp(klon),zqta(klon),zqsat(klon)
18REAL, INTENT(IN) :: pplev(klon)
19REAL, INTENT(IN) :: ztemp(klon),zqta(klon)
20REAL, INTENT(OUT) :: zqsat(klon)
21LOGICAL, INTENT(IN) :: active(klon)
22
23! Variables locales
24INTEGER ig,iter
25REAL Tbef(klon),DT(klon)
26REAL tdelta,qsatbef,zcor,qlbef,zdelta,zcvm5,dqsat,num,denom,dqsat_dT
27logical Zsat
28!REAL RLvCp
29
30REAL, SAVE :: DDT0=.01
31!$OMP THREADPRIVATE(DDT0)
32
33     !$acc data create (dt, tbef) &
34     !$acc present(pplev, ztemp, zqta, active) &
35     !$acc present(zqsat) &
36     !$acc &
37
38!====================================================================
39! INITIALISATIONS
40!====================================================================
41
42!RLvCp = RLVTT/RCPD
43!$acc kernels default(none)
44DT(:)=0.
45!!LF !$acc end kernels
46
47!====================================================================
48! Routine a vectoriser en copiant active dans converge et en mettant
49! la boucle sur les iterations a l'exterieur est en mettant
50! converge= false des que la convergence est atteinte.
51!====================================================================
52
53!!LF   !$acc kernels default(none)
54do ig=1,klon
55   if (active(ig)) then
56               Tbef(ig)=ztemp(ig)
57               zdelta=MAX(0.,SIGN(1.,RTT-Tbef(ig)))
58               qsatbef= R2ES * FOEEW(Tbef(ig),zdelta)/pplev(ig)
59               qsatbef=MIN(0.5,qsatbef)
60               zcor=1./(1.-retv*qsatbef)
61               qsatbef=qsatbef*zcor
62               qlbef=max(0.,zqta(ig)-qsatbef)
63               DT(ig) = 0.5*RLvCp*qlbef
64               zqsat(ig)=qsatbef
65     endif
66enddo
67!!LF   !$acc end kernels
68
69! Traitement du cas ou il y a condensation mais faible
70! On ne condense pas mais on dit que le qsat est le qta
71!!LF    !$acc kernels default(none)
72do ig=1,klon
73   if (active(ig)) then
74      if (0.<abs(DT(ig)).and.abs(DT(ig))<=DDT0) then
75         zqsat(ig)=zqta(ig)
76      endif
77   endif
78enddo
79!!LF   !$acc end kernels
80
81!!LF   !$acc kernels default(none)
82!!   !$acc parallel loop private(tbef,dt,zqsat)
83do iter=1,10
84    do ig=1,klon
85       if (abs(DT(ig)).gt.DDT0) then
86          Tbef(ig)=Tbef(ig)+DT(ig)
87          zdelta=MAX(0.,SIGN(1.,RTT-Tbef(ig)))
88          qsatbef= R2ES * FOEEW(Tbef(ig),zdelta)/pplev(ig)
89          qsatbef=MIN(0.5,qsatbef)
90          zcor=1./(1.-retv*qsatbef)
91          qsatbef=qsatbef*zcor
92          qlbef=zqta(ig)-qsatbef
93          zdelta=MAX(0.,SIGN(1.,RTT-Tbef(ig)))
94          zcvm5=R5LES*(1.-zdelta) + R5IES*zdelta
95          zcor=1./(1.-retv*qsatbef)
96          dqsat_dT=FOEDE(Tbef(ig),zdelta,zcvm5,qsatbef,zcor)
97          num=-Tbef(ig)+ztemp(ig)+RLvCp*qlbef
98          denom=1.+RLvCp*dqsat_dT
99          zqsat(ig) = qsatbef
100          DT(ig)=num/denom
101       endif
102    enddo
103enddo
104   !$acc end kernels
105
106
107   !$acc end data
108
109return
110end
Note: See TracBrowser for help on using the repository browser.