/[MITgcm]/MITgcm_contrib/gael/pkg/smooth/smooth_ensinv.F
ViewVC logotype

Annotation of /MITgcm_contrib/gael/pkg/smooth/smooth_ensinv.F

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


Revision 1.2 - (hide annotations) (download)
Fri Oct 16 03:36:34 2009 UTC (15 years, 9 months ago) by gforget
Branch: MAIN
CVS Tags: HEAD
Changes since 1.1: +8 -8 lines
bring pkg/smooth up to date

1 gforget 1.1 #include "CPP_OPTIONS.h"
2    
3     subroutine smooth_ensinv (mythid)
4    
5    
6     IMPLICIT NONE
7     #include "SIZE.h"
8     #include "EEPARAMS.h"
9     #include "GRID.h"
10     #include "PARAMS.h"
11     c#include "tamc.h"
12     #include "smooth.h"
13     #ifdef ALLOW_PROFILES
14     #include "profiles.h"
15     #include "netcdf.inc"
16     #endif
17     #include "DYNVARS.h"
18     c#include "adcommon.h"
19     #include "optim.h"
20     #include "ecco_cost.h"
21    
22     integer myThid
23    
24     #ifdef ALLOW_SMOOTH_ENSINV
25    
26     double precision adsalt(1-olx:snx+olx,1-oly:sny+oly,nr,nsx,nsy)
27     double precision adtheta(1-olx:snx+olx,1-oly:sny+oly,nr,nsx,nsy)
28     common /addynvars_r/ adtheta, adsalt
29    
30     character*( 80) fnamegeneric
31     integer i,j,k,bi,bj
32     integer itlo,ithi
33     integer jtlo,jthi
34     _RL port_rand, port_rand_norm
35     integer iloop,iii,ii,nbRand
36     integer nbRandPrev
37     integer choice_bckgd
38    
39     _RL mytime
40     integer myiter
41     integer num_file,num_var,prof_num
42     integer iG,jG,err,fid
43     _RL tmp_lon
44     _RL prof_traj1D(NLEVELMAX), prof_traj1D_mean(NLEVELMAX)
45     _RL prof_data1D(NLEVELMAX), prof_weights1D(NLEVELMAX)
46    
47     integer nbt_in
48    
49     jtlo = mybylo(mythid)
50     jthi = mybyhi(mythid)
51     itlo = mybxlo(mythid)
52     ithi = mybxhi(mythid)
53    
54    
55     write(fnamegeneric(1:80),'(1a)')
56     & 'wvar_T_4_50levels.bin'
57     call mdsreadfield(fnamegeneric,32,'RL',nR,
58     & wthetalev,1,mythid)
59 gforget 1.2 _EXCH_XYZ_RL( wthetalev, mythid )
60 gforget 1.1 write(fnamegeneric(1:80),'(1a)')
61     & 'wvar_S_4_50levels.bin'
62     call mdsreadfield(fnamegeneric,32,'RL',nR,
63     & wsaltlev,1,mythid)
64 gforget 1.2 _EXCH_XYZ_RL( wsaltlev, mythid )
65 gforget 1.1
66     C initialize the random number generator:
67     theta(1,1,1,1,1)=port_rand(1)
68     nbRand=200
69     nbRandPrev=119
70     choice_bckgd=0
71    
72     C load the background fields:
73     if (choice_bckgd.GT.0) then
74     write(fnamegeneric(1:80),'(1a)')
75     & 'wc01_3Densinv.bckgd'
76     call mdsreadfield(fnamegeneric,32,'RL',nR,wtheta2,
77     & 1,mythid)
78     call mdsreadfield(fnamegeneric,32,'RL',nR,wsalt2,
79     & 2,mythid)
80     endif
81    
82     C main loop:
83     DO ii=1,nbRand
84     print*,"ii",ii
85    
86     C define a perturbation:
87     DO bj=jtlo,jthi
88     DO bi=itlo,ithi
89     DO k=1,Nr
90     DO j=1-OLy,sNy+OLy
91     DO i=1-OLx,sNx+OLx
92     theta(i,j,k,bi,bj)=0.
93     salt(i,j,k,bi,bj)=0.
94     if (maskC(i,j,k,bi,bj).NE.0) then
95     theta(i,j,k,bi,bj)=port_rand_norm()
96     salt(i,j,k,bi,bj)=port_rand_norm()
97     if (choice_bckgd.EQ.1) then
98     theta(i,j,k,bi,bj)=theta(i,j,k,bi,bj)+wtheta2(i,j,k,bi,bj)
99     salt(i,j,k,bi,bj)=salt(i,j,k,bi,bj)+wsalt2(i,j,k,bi,bj)
100     endif
101     endif
102     ENDDO
103     ENDDO
104     ENDDO
105     ENDDO
106     ENDDO
107 gforget 1.2 _EXCH_XYZ_RL ( theta, myThid )
108     _EXCH_XYZ_RL ( salt, myThid )
109 gforget 1.1
110     if (ii.GT.nbRandPrev) then
111    
112     write(fnamegeneric(1:80),'(1a,i3.3)')
113     & 'wc01_3Densinv.',ii
114    
115     C the smoothing itself:
116     call smooth_correl3D(theta,1,mythid)
117     call smooth_correl3D(salt,1,mythid)
118    
119     C scale the variance
120     DO bj=jtlo,jthi
121     DO bi=itlo,ithi
122     DO k=1,Nr
123     DO j=1-OLy,sNy+OLy
124     DO i=1-OLx,sNx+OLx
125     if (maskC(i,j,k,bi,bj).NE.0) then
126     theta(i,j,k,bi,bj)=theta(i,j,k,bi,bj)
127     & * wthetaLev(i,j,k,bi,bj) * maskc(i,j,k,bi,bj)
128     salt(i,j,k,bi,bj)=salt(i,j,k,bi,bj)
129     & * wsaltLev(i,j,k,bi,bj) * maskc(i,j,k,bi,bj)
130     if (choice_bckgd.EQ.2) then
131     theta(i,j,k,bi,bj)=theta(i,j,k,bi,bj)+wtheta2(i,j,k,bi,bj)
132     salt(i,j,k,bi,bj)=salt(i,j,k,bi,bj)+wsalt2(i,j,k,bi,bj)
133     endif
134     endif
135     ENDDO
136     ENDDO
137     ENDDO
138     ENDDO
139     ENDDO
140 gforget 1.2 _EXCH_XYZ_RL ( theta, myThid )
141     _EXCH_XYZ_RL ( salt, myThid )
142 gforget 1.1
143     c do the profiles_inloop operations
144     do iloop=1,nTimeSteps
145     mytime = startTime + float(iloop-1)*deltaTclock
146     call profiles_inloop( mytime, mythid )
147     enddo
148    
149     c compute the cost and write to ad*equi.bin
150     DO bj=jtlo,jthi
151     DO bi=itlo,ithi
152     do num_file=1,NFILESPROFMAX
153     fid=fiddata(num_file,bi,bj)
154     do prof_num=1,NOBSGLOB
155     if (prof_num.LE.ProfNo(num_file,bi,bj)) then
156    
157     do num_var=1,NVARMAX
158    
159     do K=1,NLEVELMAX
160     prof_traj1D(k)=0.
161     prof_mask1D_cur(k,bi,bj)=0.
162     prof_data1D(k)=0.
163     prof_weights1D(k)=0.
164     enddo
165    
166     if (vec_quantities(num_file,num_var,bi,bj).EQV..TRUE.) then
167     do K=1,ProfDepthNo(num_file,bi,bj)
168     prof_traj1D_mean(K)=0.
169     enddo
170     call active_read_profile(num_file,
171     & ProfDepthNo(num_file,bi,bj),prof_traj1D,num_var,
172     & prof_num,.false.,optimcycle,bi,bj,mythid,
173     & profiles_dummy(num_file,num_var,bi,bj))
174     call profiles_readvector(num_file,num_var,
175     & prof_ind_glob(num_file,prof_num,bi,bj),
176     & ProfDepthNo(num_file,bi,bj),prof_data1D,bi,bj,myThid)
177     call profiles_readvector(num_file,-num_var,
178     & prof_ind_glob(num_file,prof_num,bi,bj),
179     & ProfDepthNo(num_file,bi,bj),
180     & prof_weights1D,bi,bj,myThid)
181    
182     do K=1,ProfDepthNo(num_file,bi,bj)
183     prof_traj1D(K)=
184     & prof_weights1D(K)*prof_mask1D_cur(K,bi,bj)
185     & *(prof_traj1D(K)-prof_data1D(K)-prof_traj1D_mean(K))
186     & *(prof_traj1D(K)-prof_data1D(K)-prof_traj1D_mean(K))
187     enddo
188     call adactive_read_profile( num_file,profdepthno(num_file,bi,bj),
189     &num_var,prof_num, .false. ,optimcycle,bi,bj,mythid,prof_traj1d )
190     endif
191     enddo !do num_var...
192     endif !if (prof_num.LE.ProfNo(num_file,bi,bj)) then
193     enddo !do prof_num=..
194     enddo !do num_file...
195     ENDDO
196     ENDDO
197    
198     c initialise adtheta/adsalt and do adinloop
199     DO bj=jtlo,jthi
200     DO bi=itlo,ithi
201     DO k=1,Nr
202     DO j=1-OLy,sNy+OLy
203     DO i=1-OLx,sNx+OLx
204     adtheta(i,j,k,bi,bj)=0.
205     adsalt(i,j,k,bi,bj)=0.
206     ENDDO
207     ENDDO
208     ENDDO
209     ENDDO
210     ENDDO
211 gforget 1.2 _EXCH_XYZ_RL ( adtheta, myThid )
212     _EXCH_XYZ_RL ( adsalt, myThid )
213 gforget 1.1 do iloop = ntimesteps, 1, -1
214     mytime = starttime+float(iloop-1)*deltatclock
215     call adprofiles_inloop( mytime,mythid )
216     end do
217     c co adexch and exch
218     call adexch_xyz_rl( mythid,adtheta )
219     call adexch_xyz_rl( mythid,adsalt )
220     call exch_xyz_rl( adtheta,mythid )
221     call exch_xyz_rl( adsalt,mythid )
222    
223     c do smooth_diff
224     c nbt_in=wc01_nbt(smoothOpNbCur)
225     c call smooth_diff3D(adtheta,nbt_in,mythid)
226     c call smooth_diff3D(adsalt,nbt_in,mythid)
227    
228     call mdswritefield(fnamegeneric,32,.false.,'RL',
229     & nR,theta,1,1,mythid)
230     call mdswritefield(fnamegeneric,32,.false.,'RL',
231     & nR,adtheta,2,1,mythid)
232     call mdswritefield(fnamegeneric,32,.false.,'RL',
233     & nR,salt,3,1,mythid)
234     call mdswritefield(fnamegeneric,32,.false.,'RL',
235     & nR,adsalt,4,1,mythid)
236    
237     endif
238    
239     ENDDO
240    
241     stop
242    
243     #endif / * ALLOW_SMOOTH_ENSINV * /
244    
245     end

  ViewVC Help
Powered by ViewVC 1.1.22