/[MITgcm]/MITgcm/verification/fizhi-gridalt-hs/code/fizhi_init_vars.F
ViewVC logotype

Diff of /MITgcm/verification/fizhi-gridalt-hs/code/fizhi_init_vars.F

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

revision 1.1 by molod, Fri Jun 4 16:20:00 2004 UTC revision 1.4 by jmc, Tue Mar 27 15:49:37 2012 UTC
# Line 1  Line 1 
1  C $Header$  C $Header$
2  C $Name$  C $Name$
3    
4         subroutine fizhi_init_vars (myThid)  #include "FIZHI_OPTIONS.h"
5           SUBROUTINE FIZHI_INIT_VARS (myThid)
6  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
7  c  Routine to initialise the fizhi state.  c  Routine to initialise the fizhi state.
8  c    c
9  c  Input: myThid       - Process number calling this routine  c  Input: myThid       - Process number calling this routine
10  c  c
11  c  Notes:  c  Notes:
12  c   1) For a Cold Start -  c   1) For a Cold Start -
13  c      This routine takes the initial condition on the dynamics grid  c      This routine takes the initial condition on the dynamics grid
14  c      and interpolates to the physics grid to initialize the state  c      and interpolates to the physics grid to initialize the state
15  c      variables that are on both grids. It initializes the variables  c      variables that are on both grids. It initializes the variables
# Line 19  c   3) The velocity component physics fi Line 20  c   3) The velocity component physics fi
20  c  c
21  c Calls: dyn2phys (x4)  c Calls: dyn2phys (x4)
22  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
23         implicit none         IMPLICIT NONE
 #include "CPP_OPTIONS.h"  
24  #include "SIZE.h"  #include "SIZE.h"
25  #include "fizhi_SIZE.h"  #include "fizhi_SIZE.h"
26  #include "land_SIZE.h"  #include "fizhi_land_SIZE.h"
27  #include "GRID.h"  #include "GRID.h"
28  #include "DYNVARS.h"  #include "DYNVARS.h"
29  #include "gridalt_mapping.h"  #include "gridalt_mapping.h"
30  #include "fizhi_coms.h"  #include "fizhi_coms.h"
31  #include "land_coms.h"  #include "fizhi_land_coms.h"
32    #include "fizhi_earth_coms.h"
33  #include "EEPARAMS.h"  #include "EEPARAMS.h"
34  #include "SURFACE.h"  #include "SURFACE.h"
35  #include "PARAMS.h"  #include "PARAMS.h"
36    #include "chronos.h"
37    
38         integer myThid         INTEGER myThid
39    
40  c pe on dynamics and physics grid refers to bottom edge  c pe on dynamics and physics grid refers to bottom edge
41         _RL pephy(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nrphys+1,nSx,nSy)         _RL pephy(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nrphys+1,nSx,nSy)
# Line 41  c pe on dynamics and physics grid refers Line 43  c pe on dynamics and physics grid refers
43         _RL windphy(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nrphys,nSx,nSy)         _RL windphy(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nrphys,nSx,nSy)
44         _RL udyntemp(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)         _RL udyntemp(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
45         _RL vdyntemp(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)         _RL vdyntemp(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
46           _RL tempphy(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nrphys,nSx,nSy)
47    
48         integer i, j, L, bi, bj, Lbotij         INTEGER i, j, L, bi, bj, Lbotij
49         integer im1, im2, jm1, jm2, idim1, idim2, jdim1, jdim2         INTEGER im1, im2, jm1, jm2, idim1, idim2, jdim1, jdim2
50           LOGICAL alarm
51           EXTERNAL alarm
52    
53         im1 = 1-OLx         im1 = 1-OLx
54         im2 = sNx+OLx         im2 = sNx+OLx
# Line 54  c pe on dynamics and physics grid refers Line 59  c pe on dynamics and physics grid refers
59         jdim1 = 1         jdim1 = 1
60         jdim2 = sNy         jdim2 = sNy
61    
62        IF ( startTime.EQ.0. .AND. nIter0.EQ.0 ) THEN  c   First Check to see if we can start a fizhi experiment at current time
63        print *,' In fizhi_init_vars: Cold start '  c    All Fizhi alarms must be on for the first time step of a segment
64    
65          if( .not.alarm('moist') .or. .not.alarm('turb')   .or.
66         &    .not.alarm('radsw') .or. .not.alarm('radlw') ) then
67           print *,' Cant Start Fizhi experiment at ',nymd,' ',nhms
68           stop
69          endif
70    
71    C Deal Here with Variables that are on a Fizhi Pickup or need Initialization
72    
73          IF ( startTime.EQ.baseTime .AND. nIter0.EQ.0 ) THEN
74          print *,' In fizhi_init_vars: Beginning of New Experiment '
75    
76         do bj = myByLo(myThid), myByHi(myThid)         do bj = myByLo(myThid), myByHi(myThid)
77         do bi = myBxLo(myThid), myBxHi(myThid)         do bi = myBxLo(myThid), myBxHi(myThid)
# Line 70  C Build pressures on dynamics grid Line 86  C Build pressures on dynamics grid
86          enddo          enddo
87          do j = 1,sNy          do j = 1,sNy
88          do i = 1,sNx          do i = 1,sNx
89           Lbotij = ksurfC(i,j,bi,bj)           Lbotij = kSurfC(i,j,bi,bj)
90           if(Lbotij.ne.0.)           if(Lbotij.ne.0.)
91       .    pedyn(i,j,Lbotij,bi,bj) = Ro_surf(i,j,bi,bj) + etaH(i,j,bi,bj)       &    pedyn(i,j,Lbotij,bi,bj) = Ro_surf(i,j,bi,bj) + etaH(i,j,bi,bj)
92          enddo          enddo
93          enddo          enddo
94          do j = 1,sNy          do j = 1,sNy
95          do i = 1,sNx          do i = 1,sNx
96           Lbotij = ksurfC(i,j,bi,bj)           Lbotij = kSurfC(i,j,bi,bj)
97           do L = Lbotij+1,Nr+1           do L = Lbotij+1,Nr+1
98            pedyn(i,j,L,bi,bj) = pedyn(i,j,L-1,bi,bj) -            pedyn(i,j,L,bi,bj) = pedyn(i,j,L-1,bi,bj) -
99       .                        drF(L-1)*hfacC(i,j,L-1,bi,bj)       &                        drF(L-1)*hfacC(i,j,L-1,bi,bj)
100           enddo           enddo
101  c Do not use a zero field as the top edge pressure for interpolation  c Do not use a zero field as the top edge pressure for interpolation
102           if(pedyn(i,j,Nr+1,bi,bj).lt.1.e-5)           if(pedyn(i,j,Nr+1,bi,bj).lt.1.e-5)
103       .                               pedyn(i,j,Nr+1,bi,bj) = 1.e-5       &                               pedyn(i,j,Nr+1,bi,bj) = 1.e-5
104          enddo          enddo
105          enddo          enddo
106  C Build pressures on physics grid  C Build pressures on physics grid
# Line 96  C Build pressures on physics grid Line 112  C Build pressures on physics grid
112           enddo           enddo
113  c Do not use a zero field as the top edge pressure for interpolation  c Do not use a zero field as the top edge pressure for interpolation
114           if(pephy(i,j,Nrphys+1,bi,bj).lt.1.e-5)           if(pephy(i,j,Nrphys+1,bi,bj).lt.1.e-5)
115       .                               pephy(i,j,Nrphys+1,bi,bj) = 1.e-5       &                               pephy(i,j,Nrphys+1,bi,bj) = 1.e-5
116          enddo          enddo
117          enddo          enddo
118  c  c
# Line 107  c   do units and get u = .025*ln(dP*10), Line 123  c   do units and get u = .025*ln(dP*10),
123          do j = 1,sNy          do j = 1,sNy
124          do i = 1,sNx          do i = 1,sNx
125           windphy(i,j,L,bi,bj) = 0.025 *           windphy(i,j,L,bi,bj) = 0.025 *
126       .             log((pephy(i,j,1,bi,bj)-pephy(i,j,L+1,bi,bj))*10.)       &             log((pephy(i,j,1,bi,bj)-pephy(i,j,L+1,bi,bj))*10.)
127          enddo          enddo
128          enddo          enddo
129          enddo          enddo
130    
131         enddo         enddo
132         enddo         enddo
133                                                                                      
134  c Create initial fields on phys. grid - Move Dynamics u and v to A-Grid  c Create initial fields on phys. grid - Move Dynamics u and v to A-Grid
135         call CtoA(myThid,uvel,vvel,maskW,maskS,im1,im2,jm1,jm2,Nr,         call CtoA(myThid,uvel,vvel,maskW,maskS,im1,im2,jm1,jm2,Nr,
136       .                     Nsx,Nsy,1,sNx,1,sNy,udyntemp,vdyntemp)       &                     nSx,nSy,1,sNx,1,sNy,udyntemp,vdyntemp)
137    
138         do bj = myByLo(myThid), myByHi(myThid)         do bj = myByLo(myThid), myByHi(myThid)
139         do bi = myBxLo(myThid), myBxHi(myThid)         do bi = myBxLo(myThid), myBxHi(myThid)
140    
141  c Create initial fields on phys. grid - interpolate from dyn. grid  c Create initial fields on phys. grid - interpolate from dyn. grid
142          call dyn2phys(udyntemp,pedyn,im1,im2,jm1,jm2,Nr,Nsx,Nsy,          call dyn2phys(udyntemp,pedyn,im1,im2,jm1,jm2,Nr,nSx,nSy,
143       .    1,sNx,1,sNy,bi,bj,windphy,pephy,ksurfC,Nrphys,nlperdyn,1,uphy)       & 1,sNx,1,sNy,bi,bj,windphy,pephy,kSurfC,Nrphys,nlperdyn,1,tempphy)
144          call dyn2phys(vdyntemp,pedyn,im1,im2,jm1,jm2,Nr,Nsx,Nsy,  c   Note: Interpolation gives bottom-up arrays (level 1 is bottom),
145       .    1,sNx,1,sNy,bi,bj,windphy,pephy,ksurfC,Nrphys,nlperdyn,1,vphy)  c         Physics works top-down. so -> need to flip arrays
146          call dyn2phys(theta,pedyn,im1,im2,jm1,jm2,Nr,Nsx,Nsy,          do L = 1,Nrphys
147       .   1,sNx,1,sNy,bi,bj,windphy,pephy,ksurfC,Nrphys,nlperdyn,0,thphy)          do j = 1,sNy
148          call dyn2phys(salt,pedyn,im1,im2,jm1,jm2,Nr,Nsx,Nsy,          do i = 1,sNx
149       .    1,sNx,1,sNy,bi,bj,windphy,pephy,ksurfC,Nrphys,nlperdyn,0,sphy)           uphy(i,j,Nrphys+1-L,bi,bj) = tempphy(i,j,L,bi,bj)
   
 c Now initialize tke, xlmt, khmt, xxmt, yymt, ctmt, zetamt,  
 c Now initialize land state too - tcanopy, etc...  
        do L = 1,Nrphys  
         do i = 1,nchp  
          tke(i,L,bi,bj) = 0.  
          xlmt(i,L,bi,bj) = 0.  
          khmt(i,L,bi,bj) = 0.  
150          enddo          enddo
151         enddo          enddo
152  c Now initialize land state too - tcanopy, etc... ZERO FOR NOW,          enddo
153            call dyn2phys(vdyntemp,pedyn,im1,im2,jm1,jm2,Nr,nSx,nSy,
154         & 1,sNx,1,sNy,bi,bj,windphy,pephy,kSurfC,Nrphys,nlperdyn,1,tempphy)
155            do L = 1,Nrphys
156            do j = 1,sNy
157            do i = 1,sNx
158             vphy(i,j,Nrphys+1-L,bi,bj) = tempphy(i,j,L,bi,bj)
159            enddo
160            enddo
161            enddo
162            call dyn2phys(theta,pedyn,im1,im2,jm1,jm2,Nr,nSx,nSy,
163         & 1,sNx,1,sNy,bi,bj,windphy,pephy,kSurfC,Nrphys,nlperdyn,0,tempphy)
164            do L = 1,Nrphys
165            do j = 1,sNy
166            do i = 1,sNx
167             thphy(i,j,Nrphys+1-L,bi,bj) = tempphy(i,j,L,bi,bj)
168            enddo
169            enddo
170            enddo
171    
172            call dyn2phys(salt,pedyn,im1,im2,jm1,jm2,Nr,nSx,nSy,
173         & 1,sNx,1,sNy,bi,bj,windphy,pephy,kSurfC,Nrphys,nlperdyn,0,tempphy)
174            do L = 1,Nrphys
175            do j = 1,sNy
176            do i = 1,sNx
177             sphy(i,j,Nrphys+1-L,bi,bj) = tempphy(i,j,L,bi,bj)
178            enddo
179            enddo
180            enddo
181    
182    c Zero out fizhi tendency arrays on the fizhi grid
183            do L = 1,Nrphys
184            do j = 1,sNy
185            do i = 1,sNx
186             duphy(i,j,L,bi,bj) = 0.
187             dvphy(i,j,L,bi,bj) = 0.
188             dthphy(i,j,L,bi,bj) = 0.
189             dsphy(i,j,L,bi,bj) = 0.
190            enddo
191            enddo
192            enddo
193    
194    c Zero out fizhi tendency arrays on the dynamics grid
195            do L = 1,Nr
196            do j = jm1,jm2
197            do i = im1,im2
198             guphy(i,j,L,bi,bj) = 0.
199             gvphy(i,j,L,bi,bj) = 0.
200             gthphy(i,j,L,bi,bj) = 0.
201             gsphy(i,j,L,bi,bj) = 0.
202            enddo
203            enddo
204            enddo
205    
206    c Initialize vegetation tile tke, xlmt, khmt, xxmt, yymt, ctmt, zetamt,
207            if( (nhms.eq.nhms0) .and. (nymd.eq.nymd0) ) then
208             print *,' Cold Start: Zero out Turb second moments '
209             do i = 1,nchp
210              ctmt(i,bi,bj) = 0.
211              xxmt(i,bi,bj) = 0.
212              yymt(i,bi,bj) = 0.
213              zetamt(i,bi,bj) = 0.
214             enddo
215             do L = 1,Nrphys
216             do i = 1,nchp
217              tke(i,L,bi,bj) = 0.
218              xlmt(i,L,bi,bj) = 0.
219              khmt(i,L,bi,bj) = 0.
220             enddo
221             enddo
222            else
223             print *,' Need initial Values for TKE - dont have them! '
224             stop
225            endif
226            turbStart(bi,bj) = .TRUE.
227    
228    c Now initialize vegetation tile land state too - tcanopy, etc...
229    c       call fizhi_init_vegsurftiles( nymd,nhms, 'D', myThid )
230    c Now initialize land state too - tcanopy, etc... SET FOR NOW,
231  c                                              READ CLIM FOR REAL  c                                              READ CLIM FOR REAL
232         do i = 1,nchp          do i = 1,nchp
233          tcanopy(i) = 0.          tcanopy(i,bi,bj) = 283.
234          tdeep(i) = 0.          tdeep(i,bi,bj) = 282.5
235          ecanopy(i) = 0.          ecanopy(i,bi,bj) = 2.e-2
236          swetshal(i) = 0.          swetshal(i,bi,bj) = 0.6
237          swetroot(i) = 0.          swetroot(i,bi,bj) = 0.5
238          swetdeep(i) = 0.          swetdeep(i,bi,bj) = 0.5
239          capac(i) = 0.          capac(i,bi,bj) = 0.
240          snodep(i) = 0.          snodep(i,bi,bj) = 0.
241         enddo          enddo
242    
243    c Now initialize fizhi arrays that will be on a pickup
244            print *,' Initialize fizhi arrays that will be on pickup '
245            imstturblw(bi,bj) = 0
246            imstturbsw(bi,bj) = 0
247            iras(bi,bj) = 0
248            nlwcld(bi,bj) = 0
249            nlwlz(bi,bj) = 0
250            nswcld(bi,bj) = 0
251            nswlz(bi,bj) = 0
252            do L = 1,Nrphys
253            do j = 1,sNy
254            do i = 1,sNx
255             swlz(i,j,L,bi,bj) = 0.
256             lwlz(i,j,L,bi,bj) = 0.
257             qliqavesw(i,j,L,bi,bj) = 0.
258             qliqavelw(i,j,L,bi,bj) = 0.
259             fccavesw(i,j,L,bi,bj) = 0.
260             fccavelw(i,j,L,bi,bj) = 0.
261             cldtot_sw(i,j,L,bi,bj) = 0.
262             cldras_sw(i,j,L,bi,bj) = 0.
263             cldlsp_sw(i,j,L,bi,bj) = 0.
264             cldtot_lw(i,j,L,bi,bj) = 0.
265             cldras_lw(i,j,L,bi,bj) = 0.
266             cldlsp_lw(i,j,L,bi,bj) = 0.
267            enddo
268            enddo
269            enddo
270            do j = 1,sNy
271            do i = 1,sNx
272             rainlsp(i,j,bi,bj) = 0.
273             raincon(i,j,bi,bj) = 0.
274             snowfall(i,j,bi,bj) = 0.
275            enddo
276            enddo
277    
278         enddo         enddo
279         enddo         enddo
280    
281        ELSE        ELSE
282        print *,' In fizhi_init_vars: Read from restart '        print *,' In fizhi_init_vars: Read from restart '
283                                                                                      
284  C--   Read fizhi package state variables from pickup file  C--   Read fizhi package state variables from pickup file
285          call fizhi_read_pickup( nIter0, myThid )  
286                                                                                             call fizhi_read_pickup( nIter0, myThid )
287           CALL FIZHI_READ_VEGTILES( nIter0, 'D', myThid )
288           do bj = myByLo(myThid), myByHi(myThid)
289           do bi = myBxLo(myThid), myBxHi(myThid)
290             turbStart(bi,bj) = .FALSE.
291           enddo
292           enddo
293    
294        ENDIF        ENDIF
295    
296         return        RETURN
297         end        END
Legend:
Removed from v.1.1  
changed lines
  Added in v.1.4
  ViewVC Help
Powered by ViewVC 1.1.22