source: trunk/LMDZ.VENUS/libf/phyvenus/lwi.F @ 1572

Last change on this file since 1572 was 1530, checked in by emillour, 9 years ago

Venus and Titan GCMs:
Updates in the physics to keep up with updates in LMDZ5 (up to
LMDZ5 trunk, rev 2350) concerning dynamics/physics separation:

  • Adapted makelmdz and makelmdz_fcm script to stop if trying to compile 1d model or newstart or start2archive in parallel.
  • got rid of references to "dimensions.h" in physics. Within physics packages, use nbp_lon (=iim), nbp_lat (=jjmp1) and nbp_lev (=llm) from module mod_grid_phy_lmdz (in phy_common) instead. Only partially done for Titan, because of many hard-coded commons; a necessary first step will be to clean these up (using modules).

EM
File size: 4.4 KB
Line 
1      subroutine lwi(nl,netrad,deltapsi,deltap,temp,coolrate)
2
3      use dimphy
4      use cpdet_mod, only: cpdet
5      implicit none
6
7
8#include "YOMCST.h"
9#include "timerad.h"
10     
11CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
12C
13C                             -   lwi    -   
14C
15C     PURPOSE:       Schema semi - implicite
16C
17C
18CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
19
20
21c************************************************************************
22c
23c        0.    Declarations
24c              ------------
25c
26c-------------------------------------------------------------------------
27c        0.1   Arguments
28c              ---------
29c Inputs
30c.......
31
32      integer nl
33      real    deltapsi(0:nl+1,0:nl+1) ! Dpsi/DT = sum(nu)[ksi DB/DT]
34      real    netrad(0:nl)     ! radiative budget (W/m**2)
35      real    deltap(nl) !epaisseur de la couche en pression (Pa)
36! ADAPTATION GCM POUR CP(T)
37      real    temp(nl)   !temperature de la couche
38     
39c Outputs
40c........
41
42      real    coolrate(nl)
43     
44c-------------------------------------------------------------------------
45c        0.2   local arrays
46c              ------------
47c
48      real    di(nl)
49     .      , spec(nl)
50     .      , hi(nl)
51     .      , bi(nl)
52
53      real    ci(nl)
54     .      , ai(nl)
55
56      real    deltat
57
58      real    semi,semit,denom
59
60      integer i,jl,band
61
62c************************************************************************
63c
64c        1.    Initialisations
65c              ---------------
66c
67c-----------------------------------------------------------------------
68        deltat = dtimerad
69c       print*,'SEMI = ',semi, '(expl:0  semi-implicite:0.5  impl:1)'
70c        semi  = 0.5
71c       semi = 0.5
72        semi = 1.
73        semit = semi * deltat
74
75c       print*,'dtimerad,deltat,semit:',dtimerad,deltat,semit
76
77c************************************************************************
78c
79c        2.   
80c              ---------------
81c
82c-------------------------------------------------------------------------
83c        2.1   Calcul des di
84c              -------------
85c
86
87      do i = 1 , nl-1
88        spec(i) =
89     .      ( deltapsi(i,nl+1)
90     .      + deltapsi(i,i+1)
91     .      + deltapsi(i,i-1) )
92     
93        di(i) =  1. + semit * (RG/deltap(i)/cpdet(temp(i))) * spec(i)
94c     print*,i,' di(i)=',di(i)
95      enddo
96         
97
98c couche nl
99c ------------
100c      , on enleve i,i+1 sinon on a 2 fois le cooling2space
101 
102        spec(nl) =
103     .      ( deltapsi(nl,nl+1)
104     .      + deltapsi(nl,nl-1) )
105     
106        di(nl) =  1. + semit * (RG/deltap(nl)/cpdet(temp(nl)))*spec(nl)
107
108c-------------------------------------------------------------------------
109c        2.2   Calcul des hi
110c              -------------
111c
112
113      do i = 1 , nl-1
114     
115        spec(i) = deltapsi(i+1,i)
116
117        hi(i) = - semit * (RG/deltap(i)/cpdet(temp(i))) * spec(i)
118
119      enddo
120                       
121c     print*,'hi(i)',hi(i)
122
123c-------------------------------------------------------------------------
124c        2.3   Calcul des bi
125c              -------------
126c
127
128
129      do i = 2 , nl
130     
131        spec(i) = deltapsi(i-1,i)
132       
133        bi(i) = - semit * (RG/deltap(i)/cpdet(temp(i))) * spec(i)
134       
135      enddo
136
137c     print*,'bi(i)',bi(i)
138
139c couche 1
140c --------
141c  tant qu'on a pas un calcul propre de deltab(0) qui tienne compte de
142c    la discontinuite de temperature au sol , on met  b1 = 0
143
144     
145        bi(1) = 0
146
147c-------------------------------------------------------------------------
148c        2.4   
149c              -------------
150c
151
152c couche nl
153c ------------
154
155c     ci(nl) = ((RG/CP) * netrad(nl) / deltap(nl))
156c    .                   / di(nl)
157      ci(nl) = netrad(nl) *RG/cpdet(temp(nl)) / deltap(nl)
158
159      ai(nl) = - bi(nl) / di(nl)
160
161      do i = nl-1 , 1 , -1
162        denom = di(i) + hi(i) * ai(i+1)
163
164c       ci(i) = (  (RG/CP) * netrad(i) / deltap(i)
165c    .               - hi(i) * ci(i+1)  )  / denom
166        ci(i) = netrad(i) *RG/cpdet(temp(i)) / deltap(i)
167 
168        ai(i) = -bi(i) / denom
169      enddo
170
171c-------------------------------------------------------------------------
172c        2.5   
173c              -------------
174c
175
176      coolrate(1) = ci(1)
177
178      do i = 2 , nl
179           coolrate(i) = ci(i) + ai(i) * coolrate(i-1)
180c       print*,i,' coolrate(i)=',coolrate(i)
181      enddo
182
183c-------------------------------------------------------------------------
184
185      RETURN
186      END
Note: See TracBrowser for help on using the repository browser.