/[MITgcm]/MITgcm/model/src/ini_p_ground.F
ViewVC logotype

Diff of /MITgcm/model/src/ini_p_ground.F

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph | View Patch Patch

revision 1.3 by jmc, Wed Nov 6 03:45:46 2002 UTC revision 1.4 by jmc, Tue Dec 10 02:55:47 2002 UTC
# Line 2  C $Header$ Line 2  C $Header$
2  C $Name$  C $Name$
3    
4  #include "CPP_OPTIONS.h"  #include "CPP_OPTIONS.h"
5    #undef CHECK_ANALYLIC_THETA
6    
7  CBOP  CBOP
8  C     !ROUTINE: INI_P_GROUND  C     !ROUTINE: INI_P_GROUND
9  C     !INTERFACE:  C     !INTERFACE:
10        SUBROUTINE INI_P_GROUND(        SUBROUTINE INI_P_GROUND(selectMode,
11       I                        Hfld,       &                        Hfld, Pfld,
      O                        Pfld,  
12       I                        myThid )       I                        myThid )
13  C     !DESCRIPTION: \bv  C     !DESCRIPTION: \bv
14  C     *==========================================================*  C     *==========================================================*
# Line 27  C     == Global variables == Line 27  C     == Global variables ==
27  #include "GRID.h"  #include "GRID.h"
28  #include "EEPARAMS.h"  #include "EEPARAMS.h"
29  #include "PARAMS.h"  #include "PARAMS.h"
30    #include "SURFACE.h"
31    
32  C     !INPUT/OUTPUT PARAMETERS:  C     !INPUT/OUTPUT PARAMETERS:
33  C     == Routine arguments ==  C     == Routine arguments ==
34  C     Hfld (input)  :: Ground elevation [m]  C     selectMode :: 0 = find Pfld from Hfld using Theta_Ref profile
35  C     Pfld (output) :: reference Pressure at the ground [Pa]  C                   1 = find Pfld from Hfld using Analytic fct Theta(Lat,P)
36    C                  -1 = compute Hfld from Pfld using Analytic  Theta(Lat,P)
37    C     Hfld (input/outp) :: Ground elevation [m]
38    C     Pfld (outp/input) :: reference Pressure at the ground [Pa]
39          INTEGER selectMode
40        _RS Hfld(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)        _RS Hfld(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
41        _RS Pfld(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)        _RS Pfld(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
42        INTEGER myThid        INTEGER myThid
43    
 c #ifdef INCLUDE_PHIHYD_CALCULATION_CODE  
   
44  C     !LOCAL VARIABLES:  C     !LOCAL VARIABLES:
45  C     == Local variables ==  C     == Local variables ==
46    C     msgBuf :: Informational/error meesage buffer
47    C-
48    C       For an accurate definition of the reference surface pressure,
49    C       define a High vertical resolution (in P):
50    C     nLevHvR :: Number of P-level used for High vertical Resolution (HvR)
51    C     plowHvR :: Lower bound of the High vertical Resolution
52    C     dpHvR   :: normalized pressure increment (HvR)
53    C     pLevHvR :: normalized P-level of the High vertical Resolution
54    C     pMidHvR :: normalized mid point level (HvR)
55    C     thetaHvR :: potential temperature at mid point level (HvR)
56    C     PiHvR  :: Exner function at P-level
57    C     dPiHvR :: Exner function difference between 2 P-levels
58    C     recip_kappa :: 1/kappa = Cp/R
59    C     PiLoc, zLoc, dzLoc, yLatLoc, phiLoc :: hold on temporary values
60    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
61          CHARACTER*(MAX_LEN_MBUF) msgBuf
62        INTEGER bi,bj,i,j,K, Ks        INTEGER bi,bj,i,j,K, Ks
63        _RL ddPI, Po_surf        _RL ddPI, Po_surf
64        _RL phiRef(2*Nr+1), rHalf(2*Nr+1)        _RL phiRef(2*Nr+1), rHalf(2*Nr+1)
65    
66          INTEGER nLevHvR
67          PARAMETER ( nLevHvR = 60 )
68          _RL plowHvR, dpHvR, pLevHvR(nLevHvR+1), pMidHvR(nLevHvR)
69          _RL thetaHvR(nLevHvR), PiHvR(nLevHvR+1), dPiHvR(nLevHvR)
70          _RL recip_kappa, PiLoc, zLoc, dzLoc, yLatLoc, phiLoc
71          INTEGER kLev
72    #ifdef CHECK_ANALYLIC_THETA
73          _RL tmpVar(nLevAn,61)
74    #endif
75  CEOP  CEOP
76    
77          IF ( selectFindRoSurf.LT.0 .OR. selectFindRoSurf.GT.1 ) THEN
78            WRITE(msgBuf,'(A,I2,A)')
79         &   'INI_P_GROUND: selectFindRoSurf =', selectFindRoSurf,
80         &        ' <== bad value !'
81            CALL PRINT_ERROR( msgBuf , myThid)
82            STOP 'INI_P_GROUND'
83          ENDIF
84    
85        DO K=1,Nr        DO K=1,Nr
86          rHalf(2*K-1) = rF(K)          rHalf(2*K-1) = rF(K)
87          rHalf(2*K)   = rC(K)          rHalf(2*K)   = rC(K)
88        ENDDO        ENDDO
89         rHalf(2*Nr+1) = rF(Nr+1)         rHalf(2*Nr+1) = rF(Nr+1)
90    
91          IF (selectMode .LE. 0) THEN
92  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
93  C- Compute Reference Geopotential at Half levels :  C- Compute Reference Geopotential at Half levels :
94  C    Tracer level: phiRef(2K)  ;  Interface_W level: phiRef(2K+1)  C    Tracer level: phiRef(2K)  ;  Interface_W level: phiRef(2K+1)
95    
96        phiRef(1) = 0.        phiRef(1) = 0.
97    
98        IF (Integr_GeoPot.EQ.1) THEN        IF (integr_GeoPot.EQ.1) THEN
99  C- Finite Volume Form, linear by half level :  C- Finite Volume Form, linear by half level :
100          DO K=1,2*Nr          DO K=1,2*Nr
101            Ks = (K+1)/2            Ks = (K+1)/2
102            ddPI=atm_cp*( ((rHalf( K )/atm_po)**atm_kappa)            ddPI=atm_Cp*( ((rHalf( K )/atm_Po)**atm_kappa)
103       &                 -((rHalf(K+1)/atm_po)**atm_kappa) )       &                 -((rHalf(K+1)/atm_Po)**atm_kappa) )
104            phiRef(K+1) = phiRef(K)+ddPI*tRef(Ks)            phiRef(K+1) = phiRef(K)+ddPI*tRef(Ks)
105          ENDDO          ENDDO
106  C------  C------
107        ELSE        ELSE
108  C- Finite Difference Form, linear between Tracer level :  C- Finite Difference Form, linear between Tracer level :
109  C   works with Integr_GeoPot = 0, 2 or 3  C   works with integr_GeoPot = 0, 2 or 3
110          K = 1          K = 1
111            ddPI=atm_cp*( ((rF(K)/atm_po)**atm_kappa)            ddPI=atm_Cp*( ((rF(K)/atm_Po)**atm_kappa)
112       &                 -((rC(K)/atm_po)**atm_kappa) )       &                 -((rC(K)/atm_Po)**atm_kappa) )
113            phiRef(2*K)   = phiRef(1) + ddPI*tRef(K)            phiRef(2*K)   = phiRef(1) + ddPI*tRef(K)
114          DO K=1,Nr-1          DO K=1,Nr-1
115            ddPI=atm_cp*( ((rC( K )/atm_po)**atm_kappa)            ddPI=atm_Cp*( ((rC( K )/atm_Po)**atm_kappa)
116       &                 -((rC(K+1)/atm_po)**atm_kappa) )       &                 -((rC(K+1)/atm_Po)**atm_kappa) )
117            phiRef(2*K+1) = phiRef(2*K) + ddPI*0.5*tRef(K)            phiRef(2*K+1) = phiRef(2*K) + ddPI*0.5*tRef(K)
118            phiRef(2*K+2) = phiRef(2*K)            phiRef(2*K+2) = phiRef(2*K)
119       &                  + ddPI*0.5*(tRef(K)+tRef(K+1))       &                  + ddPI*0.5*(tRef(K)+tRef(K+1))
120          ENDDO          ENDDO
121          K = Nr          K = Nr
122            ddPI=atm_cp*( ((rC( K )/atm_po)**atm_kappa)            ddPI=atm_Cp*( ((rC( K )/atm_Po)**atm_kappa)
123       &                 -((rF(K+1)/atm_po)**atm_kappa) )       &                 -((rF(K+1)/atm_Po)**atm_kappa) )
124            phiRef(2*K+1) = phiRef(2*K) + ddPI*tRef(K)            phiRef(2*K+1) = phiRef(2*K) + ddPI*tRef(K)
125  C------  C------
126         ENDIF         ENDIF
# Line 105  C- Write to check : Line 143  C- Write to check :
143        ENDDO        ENDDO
144        _END_MASTER( myThid )        _END_MASTER( myThid )
145    
146    C-- endif selectMode =< 0
147          ENDIF
148    
149          IF (selectMode .EQ. 0) THEN
150  C- Find Po_surf : Linear between 2 half levels :  C- Find Po_surf : Linear between 2 half levels :
151        DO bj = myByLo(myThid), myByHi(myThid)        DO bj = myByLo(myThid), myByHi(myThid)
152         DO bi = myBxLo(myThid), myBxHi(myThid)         DO bi = myBxLo(myThid), myBxHi(myThid)
# Line 118  C- Find Po_surf : Linear between 2 half Line 160  C- Find Po_surf : Linear between 2 half
160       &       (Hfld(i,j,bi,bj)-phiRef(Ks))/(phiRef(Ks+1)-phiRef(Ks))       &       (Hfld(i,j,bi,bj)-phiRef(Ks))/(phiRef(Ks+1)-phiRef(Ks))
161    
162  c          IF (Hfld(i,j,bi,bj).LT.phiRef(1)) Po_surf= rHalf(1)  c          IF (Hfld(i,j,bi,bj).LT.phiRef(1)) Po_surf= rHalf(1)
163  c          IF (Hfld(i,j,bi,bj).GT.phiRef(2*Nr+1)) Po_surf= rHalf(2*Nr+1)  c          IF (Hfld(i,j,bi,bj).GT.phiRef(2*Nr+1)) Po_surf=rHalf(2*Nr+1)
164             Pfld(i,j,bi,bj) = Po_surf             Pfld(i,j,bi,bj) = Po_surf
165           ENDDO           ENDDO
166          ENDDO          ENDDO
# Line 126  c          IF (Hfld(i,j,bi,bj).GT.phiRef Line 168  c          IF (Hfld(i,j,bi,bj).GT.phiRef
168        ENDDO        ENDDO
169    
170  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
171    C-- endif selectMode=0
172          ENDIF
173    
174          IF (abs(selectMode) .EQ. 1) THEN
175    C-- define high resolution Pressure discretization:
176    
177          recip_kappa = 1. _d 0 / atm_kappa
178          plowHvR = 0.4 _d 0
179          dpHvR = nLevHvR
180          dpHvR = (1. - plowHvR) / dpHvR
181            pLevHvR(1)= Ro_SeaLevel/atm_Po
182            PiHvR(1) = atm_Cp*(pLevHvR(1)**atm_kappa)
183          DO k=1,nLevHvR
184            pLevHvR(k+1)= pLevHvR(1) - float(k)*dpHvR
185            PiHvR(k+1)  = atm_Cp*(pLevHvR(k+1)**atm_kappa)
186            pMidHvR(k)= (pLevHvR(k)+pLevHvR(k+1))*0.5 _d 0
187            dPiHvR(k) = PiHvR(k) - PiHvR(k+1)
188          ENDDO
189    
190    #ifdef CHECK_ANALYLIC_THETA
191          _BEGIN_MASTER( mythid )
192          DO j=1,61
193            yLatLoc =-90.+(j-1)*3.
194            CALL ANALYLIC_THETA( yLatLoc , pMidHvR,
195         &                       tmpVar(1,j), nLevHvR, mythid)
196          ENDDO  
197          OPEN(88,FILE='analytic_theta',
198         &      STATUS='unknown',FORM='unformatted')
199          WRITE(88) tmpVar
200          CLOSE(88)
201          _END_MASTER( mythid )
202          STOP 'CHECK_ANALYLIC_THETA'
203    #endif /* CHECK_ANALYLIC_THETA */
204    
205    C-- endif |selectMode|=1
206          ENDIF
207    
208          IF (selectMode .EQ. 1) THEN
209    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
210    C- Find Po_surf such as g*Hfld = Phi[Po_surf,theta(yLat,p)]:
211    
212          DO bj = myByLo(myThid), myByHi(myThid)
213           DO bi = myBxLo(myThid), myBxHi(myThid)
214    C- start bi,bj loop:
215            DO j=1,sNy
216             DO i=1,sNx
217              IF ( Hfld(i,j,bi,bj) .LE. 0. _d 0) THEN
218               Pfld(i,j,bi,bj) = Ro_SeaLevel
219              ELSE
220               yLatLoc  = yC(i,j,bi,bj)
221               CALL ANALYLIC_THETA( yLatLoc , pMidHvR,
222         &                       thetaHvR, nLevHvR, mythid)
223               zLoc = 0.
224               DO k=1,nLevHvR
225                IF (zLoc.GE.0.) THEN
226    C-    compute  phi/g corresponding to next p_level:
227                 dzLoc = dPiHvR(k)*thetaHvR(k)*recip_gravity
228                 IF ( Hfld(i,j,bi,bj) .LE. zLoc+dzLoc ) THEN
229    C-    compute the normalized surf. Pressure Po_surf
230                   PiLoc = PiHvR(k)
231         &               - gravity*(Hfld(i,j,bi,bj)-zLoc)/thetaHvR(k)
232                   Po_surf = (PiLoc/atm_Cp)**recip_kappa
233    C- use linear interpolation:
234    c              Po_surf = pLevHvR(k)
235    c    &                 - dpHvR*(Hfld(i,j,bi,bj)-zLoc)/dzLoc
236                   zLoc = -1.
237                 ELSE
238                   zLoc = zLoc + dzLoc
239                 ENDIF
240                ENDIF
241               ENDDO
242               IF (zLoc.GE.0.) THEN
243                 WRITE(msgBuf,'(2A)')
244         &        'INI_P_GROUND: FAIL in trying to find Pfld:',
245         &        ' selectMode,i,j,bi,bj=',selectMode,i,j,bi,bj
246                 CALL PRINT_ERROR( msgBuf , myThid)
247                 WRITE(msgBuf,'(A,F7.1,2A,F6.4,A,F8.0)')
248         &        'INI_P_GROUND: Hfld=', Hfld(i,j,bi,bj), ' exceeds',
249         &        ' Zloc(lowest P=', pLevHvR(1+nLevHvR),' )=',zLoc
250                 CALL PRINT_ERROR( msgBuf , myThid)
251                 STOP 'ABNORMAL END: S/R INI_P_GROUND'  
252               ELSE
253                 Pfld(i,j,bi,bj) = Po_surf*atm_Po
254               ENDIF
255              ENDIF
256             ENDDO
257            ENDDO
258    C- end bi,bj loop.
259           ENDDO
260          ENDDO
261    
262    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
263    C-- endif selectMode=1
264          ENDIF
265    
266          IF (selectMode .EQ. -1) THEN
267    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
268    C  goal: evaluate phi0surf (used for computing geopotential'= Phi - PhiRef):
269    C   phi0surf = g*Ztopo-PhiRef(Ro_surf) if no truncation was applied to Ro_surf;
270    C  but because of hFacMin, surf.ref.pressure (=Ro_surf) is truncated and
271    C   phi0surf = Phi(Theta-Analytic,P=Ro_surf) - PhiRef(P=Ro_surf)
272    C-----
273    
274          DO bj = myByLo(myThid), myByHi(myThid)
275           DO bi = myBxLo(myThid), myBxHi(myThid)
276    C- start bi,bj loop:
277    
278    C--   Compute Hfld from g*Hfld = Phi[Po_surf,theta(yLat,p)]:
279            DO j=1,sNy
280             DO i=1,sNx
281              IF ( Pfld(i,j,bi,bj) .GE. Ro_SeaLevel) THEN
282               Hfld(i,j,bi,bj) = 0.
283              ELSE
284               Po_surf = Pfld(i,j,bi,bj)/atm_Po
285               kLev = 1 + INT( (pLevHvR(1)-Po_surf)/dpHvR )
286               yLatLoc  = yC(i,j,bi,bj)
287               CALL ANALYLIC_THETA( yLatLoc , pMidHvR,
288         &                       thetaHvR, kLev, mythid)
289    C-    compute height at level pLev(kLev) just below Pfld:
290               zLoc = 0.
291               DO k=1,kLev-1
292                 dzLoc = dPiHvR(k)*thetaHvR(k)*recip_gravity
293                 zLoc = zLoc + dzLoc
294               ENDDO
295               dzLoc = ( PiHvR(kLev)-atm_Cp*(Po_surf**atm_kappa) )
296         &         * thetaHvR(kLev)*recip_gravity
297               Hfld(i,j,bi,bj) = zLoc + dzLoc
298              ENDIF
299             ENDDO
300            ENDDO
301    
302    C--  compute phi0surf (stored in Hfld):
303            DO j=1,sNy
304             DO i=1,sNx
305    C-   compute phiLoc = PhiRef(Ro_surf):
306              ks = ksurfC(i,j,bi,bj)
307              IF (ks.LE.Nr) THEN
308               IF ( Pfld(i,j,bi,bj).GE.rC(ks) ) THEN
309                phiLoc = phiRef(2*ks)
310         &       + (phiRef(2*ks-1)-phiRef(2*ks))
311         &        *(Pfld(i,j,bi,bj)-rC(ks))/(rHalf(2*ks-1)-rHalf(2*ks))
312               ELSE
313                phiLoc = phiRef(2*ks)
314         &       + (phiRef(2*ks+1)-phiRef(2*ks))
315         &        *(Pfld(i,j,bi,bj)-rC(ks))/(rHalf(2*ks+1)-rHalf(2*ks))
316               ENDIF
317               Hfld(i,j,bi,bj) = gravity*(Hfld(i,j,bi,bj) - phiLoc)
318              ELSE
319               Hfld(i,j,bi,bj) = 0.
320              ENDIF
321             ENDDO
322            ENDDO
323    C- end bi,bj loop.
324           ENDDO
325          ENDDO
326    
327    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
328    C-- endif selectMode=-1
329          ENDIF
330    
331          RETURN
332          END
333    
334    CBOP
335    C     !SUBROUTINE: ANALYLIC_THETA
336    C     !INTERFACE:
337          SUBROUTINE ANALYLIC_THETA( yLat, pNlev,
338         O                           thetaLev,
339         I                           kSize,myThid)
340    
341    C     !DESCRIPTION: \bv
342    C     *==========================================================*
343    C     | SUBROUTINE ANALYLIC_THETA                                    
344    C     | o Conpute analyticaly the potential temperature Theta
345    C     |   as a function of Latitude (yLat) and normalized
346    C     |   pressure pNlev.
347    C     |   approximatively match the N-S symetric, zonal-mean and
348    C     |   annual-mean NCEP climatology in the troposphere.
349    C     *==========================================================*
350    C     \ev
351    
352    C     !USES:
353          IMPLICIT NONE
354    
355    C     == Global variables ==
356    #include "SIZE.h"
357    #include "EEPARAMS.h"
358    #include "PARAMS.h"
359    
360    C     !INPUT/OUTPUT PARAMETERS:
361    C     == Routine arguments ==
362    C     yLat   :: latitude (degre)
363    C     pNlev  :: normalized pressure levels
364    C     kSize  :: dimension of pNlev & ANALYLIC_THETA
365    C     myThid :: Thread number for this instance of the routine
366          INTEGER kSize
367          _RL  yLat
368          _RL  pNlev  (kSize)
369          _RL  thetaLev(kSize)
370          INTEGER myThid
371    
372    C     !LOCAL VARIABLES:
373    C     == Local variables ==
374          INTEGER k
375          _RL  yyA, yyB, yyC, yyAd, yyBd, yyCd
376          _RL  cAtmp, cBtmp, ttdC
377          _RL  ppN0, ppN1, ppN2, ppN3a, ppN3b, ppN4
378          _RL  ttp1, ttp2, ttp3, ttp4, ttp5
379          _RL  yAtmp, yBtmp, yCtmp, yDtmp
380          _RL  ttp2y, ttp4y, a1tmp
381          _RL  ppl, ppm, pph, ppr
382    
383    CEOP
384    
385          DATA yyA ,    yyB ,     yyC ,     yyAd ,   yyBd ,   yyCd
386         &  / 45. _d 0, 65. _d 0, 65. _d 0, .9 _d 0, .9 _d 0, 10. _d 0 /
387          DATA  cAtmp ,   cBtmp ,   ttdC
388         &   /  2.6 _d 0, 1.5 _d 0, 3.3 _d 0 /
389          DATA  ppN0  ,   ppN1  ,  ppN2  ,  ppN3a ,  ppN3b ,  ppN4
390         &   / .1 _d 0, .19 _d 0, .3 _d 0, .9 _d 0, .7 _d 0, .925 _d 0 /
391          DATA ttp1 ,     ttp2 ,     ttp3 ,     ttp4 ,     ttp5
392         &   / 350. _d 0, 342. _d 0, 307. _d 0, 301. _d 0, 257. _d 0 /
393    
394    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
395    
396           yAtmp = ABS(yLat) - yyA
397           yAtmp = yyA + MIN(0. _d 0,yAtmp/yyAd) + MAX(yAtmp, 0. _d 0)
398           yAtmp = COS( deg2rad*MAX(yAtmp, 0. _d 0) )
399           yBtmp = ABS(yLat) - yyB
400           yBtmp = yyB + yBtmp/yyBd
401           yBtmp = COS( deg2rad*MAX( 0. _d 0, MIN(yBtmp,90. _d 0) ) )
402           yCtmp = ABS(yLat) - yyC
403           yCtmp = MAX( 0. _d 0, 1. _d 0 - (yCtmp/yyCd)**2 )
404           yDtmp = ppN3a +(ppN3b - ppN3a)*yCtmp
405           ttp2y = ttp3 + (ttp2-ttp3)*yAtmp**cAtmp
406           ttp4y = ttp5 + (ttp4-ttp5)*yBtmp**cBtmp
407           a1tmp = (ttp1-ttp2y)*ppN1*ppN2/(ppN2-ppN1)
408          DO k=1,kSize
409           ppl = MIN(pNlev(k),ppN1)
410           ppm = MIN(MAX(pNlev(k),ppN1),ppN2)
411           pph = MAX(pNlev(k),ppN2)
412           ppr =( ppN0 + ABS(ppl-ppN0) - ppN1 )/(ppN2-ppN1)
413           thetaLev(k) =
414         &       ( (1. _d 0 -ppr)*ttp1*ppN1**atm_kappa
415         &        + ppr*ttp2y*ppN2**atm_kappa
416         &       )*ppl**(-atm_kappa)
417         &     + a1tmp*(1. _d 0 /ppm - 1. _d 0/ppN1)
418         &     + (ttp4y-ttp2y)*(pph-ppN2)/(ppN4-ppN2)
419         &     + (ttdC+yCtmp)*MAX(0. _d 0,pNlev(k)-yDtmp)/(1-yDtmp)
420          ENDDO
421    
422  c #endif/* INCLUDE_PHIHYD_CALCULATION_CODE */  C---------------------------------------------------  
423    C matlab script, input: pN, yp=abs(yLat)
424    C pN0=.1; pN1=.19 ; pN2=.3; pN4=.925;
425    C pm=min(max(pN,pN1),pN2); pp=max(pN,pN2);
426    C pl=min(pN,pN1); pr=(pN0+abs(pl-pN0)-pN1)/(pN2-pN1);
427    C  
428    C  yA=yp-45; yA=45+min(0,yA/.9)+max(0,yA); yA=max(0,yA); cosyA=cos(yA*rad) ;
429    C  yB=yp-65; yB=65+yB/.9; yB=min(max(0,yB),90); cosyB=cos(yB*rad) ;
430    C  tp1=350*ones(nyA,1);
431    C  tp2=307+(342-307)*cosyA.^2.6;
432    C  tp4=257+(301-257)*cosyB.^1.5;
433    C  a1=(tp1-tp2)*pN1*pN2/(pN2-pN1);
434    C  pF=max(0,1.-((yp-65)/10).^2); pT=.9+(0.7-.9)*pF;
435    C  
436    C  tA0=((1-pr(k))*tp1(j)*pN1^kappa+pr(k)*tp2(j)*pN2^kappa)*pl(k)^-kappa;
437    C  tA1=a1(j)*(1./pm(k)-1./pN1);
438    C  tA2=(tp4(j)-tp2(j))*(pp(k)-pN2)/(pN4-pN2);
439    C  tA3=(3.3+pF(j))*max(0,pN(k)-pT(j))/(1-pT(j));
440    C  tAn(j,k)=tA0+tA1+tA2+tA3;
441    C---------------------------------------------------  
442    
443        RETURN        RETURN
444        END        END
Legend:
Removed from v.1.3  
changed lines
  Added in v.1.4
  ViewVC Help
Powered by ViewVC 1.1.22