/[MITgcm]/MITgcm/pkg/ecco/cost_theta.F
ViewVC logotype

Diff of /MITgcm/pkg/ecco/cost_theta.F

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

revision 1.2 by heimbach, Sat Dec 20 20:11:14 2003 UTC revision 1.8 by heimbach, Sat Oct 15 18:27:32 2005 UTC
# Line 3  C $Header$ Line 3  C $Header$
3  #include "COST_CPPOPTIONS.h"  #include "COST_CPPOPTIONS.h"
4    
5    
6        subroutine cost_Theta(        subroutine cost_theta( myiter, mytime, mythid )
      I                       myiter,  
      I                       mytime,  
      I                       mythid  
      &                     )  
7    
8  c     ==================================================================  c     ==================================================================
9  c     SUBROUTINE cost_Theta  c     SUBROUTINE cost_theta
10  c     ==================================================================  c     ==================================================================
11  c  c
12  c     o Evaluate cost function contribution of temperature.  c     o Evaluate cost function contribution of temperature.
# Line 27  c Line 23  c
23  c              - set ladinit to .true. to initialise adtbar file  c              - set ladinit to .true. to initialise adtbar file
24  c  c
25  c     ==================================================================  c     ==================================================================
26  c     SUBROUTINE cost_Theta  c     SUBROUTINE cost_theta
27  c     ==================================================================  c     ==================================================================
28    
29        implicit none        implicit none
# Line 36  c     == global variables == Line 32  c     == global variables ==
32    
33  #include "EEPARAMS.h"  #include "EEPARAMS.h"
34  #include "SIZE.h"  #include "SIZE.h"
35    #include "PARAMS.h"
36  #include "GRID.h"  #include "GRID.h"
37  #include "DYNVARS.h"  #include "DYNVARS.h"
38    
# Line 53  c     == routine arguments == Line 50  c     == routine arguments ==
50    
51  c     == local variables ==  c     == local variables ==
52    
       _RS        one_rs  
       parameter( one_rs = 1. )  
   
53        integer bi,bj        integer bi,bj
54        integer i,j,k        integer i,j,k
55        integer itlo,ithi        integer itlo,ithi
56        integer jtlo,jthi        integer jtlo,jthi
57        integer jmin,jmax        integer jmin,jmax
58        integer imin,imax        integer imin,imax
59        integer irec        integer irec, irectmp
60        integer levmon        integer levmon
61        integer levoff        integer levoff
62        integer iltheta        integer iltheta
# Line 72  c     == local variables == Line 66  c     == local variables ==
66    
67        _RL cmask (1-olx:snx+olx,1-oly:sny+oly)          _RL cmask (1-olx:snx+olx,1-oly:sny+oly)  
68        _RL spval        _RL spval
69          _RL spmax
70    
71        character*(80) fnametheta        character*(80) fnametheta
72    
# Line 79  c     == local variables == Line 74  c     == local variables ==
74        logical ladinit        logical ladinit
75    
76        character*(MAX_LEN_MBUF) msgbuf        character*(MAX_LEN_MBUF) msgbuf
   
77  #ifdef GENERIC_BAR_MONTH  #ifdef GENERIC_BAR_MONTH
78        integer mrec, nyears, iyear, tmpmon        integer mrec, nyears, iyear
79  #endif  #endif      
80    
81          _RL diagnosfld3d(1-olx:snx+olx,1-oly:sny+oly,nr,nsx,nsy)
82          
83  c     == external functions ==  c     == external functions ==
84    
85        integer  ilnblnk        integer  ilnblnk
# Line 100  c     == end of interface == Line 96  c     == end of interface ==
96        imin = 1        imin = 1
97        imax = snx        imax = snx
98    
99        spval = -9990.        spval = -1.8
100          spmax = 40.
101    
102  c--   Read tiled data.  c--   Read tiled data.
103        doglobalread = .false.        doglobalread = .false.
# Line 108  c--   Read tiled data. Line 105  c--   Read tiled data.
105    
106  #ifdef ALLOW_THETA_COST_CONTRIBUTION  #ifdef ALLOW_THETA_COST_CONTRIBUTION
107    
 #ifdef ECCO_VERBOSE  
       _BEGIN_MASTER( mythid )  
         write(msgbuf,'(a)') ' '  
         call print_message( msgbuf, standardmessageunit,  
      &                      SQUEEZE_RIGHT , mythid)  
         write(msgbuf,'(a,i8.8)')  
      &    ' cost_Theta: number of records to process = ',nmonsrec  
         call print_message( msgbuf, standardmessageunit,  
      &                      SQUEEZE_RIGHT , mythid)  
         write(msgbuf,'(a)') ' '  
         call print_message( msgbuf, standardmessageunit,  
      &                      SQUEEZE_RIGHT , mythid)  
       _END_MASTER( mythid )  
 #endif  
   
108        if (optimcycle .ge. 0) then        if (optimcycle .ge. 0) then
109          iltheta = ilnblnk( tbarfile )          iltheta = ilnblnk( tbarfile )
110          write(fnametheta(1:80),'(2a,i10.10)')          write(fnametheta(1:80),'(2a,i10.10)')
# Line 132  c--   Read tiled data. Line 114  c--   Read tiled data.
114        fcthread = 0. _d 0        fcthread = 0. _d 0
115    
116  #ifdef GENERIC_BAR_MONTH  #ifdef GENERIC_BAR_MONTH
   
117  c--   Loop over month  c--   Loop over month
118        if (nmonsrec.GT.12) then        do irec = 1,min(nmonsrec,12)
          tmpmon = 12  
       else  
          tmpmon = nmonsrec  
       endif  
       do irec = 1,tmpmon  
119           nyears=int((nmonsrec-irec)/12)+1           nyears=int((nmonsrec-irec)/12)+1
          if(nyears.gt.0) then  
120              do iyear=1,nyears              do iyear=1,nyears
121              mrec=irec+(iyear-1)*12              mrec=irec+(iyear-1)*12
122            print *, 'ph-rec irec, iyear, mrec ',              irectmp=mrec
      &               irec, iyear, mrec  
123  c--     Read time averages and the monthly mean data.  c--     Read time averages and the monthly mean data.
124              call active_read_xyz( fnametheta, tbar, mrec,              call active_read_xyz( fnametheta, tbar, mrec,
125       &           doglobalread, ladinit,       &           doglobalread, ladinit,
# Line 162  c--     Read time averages and the month Line 136  c--     Read time averages and the month
136                                tbar(i,j,k,bi,bj) =(tbar_gen(i,j,k,bi,bj)                                tbar(i,j,k,bi,bj) =(tbar_gen(i,j,k,bi,bj)
137       $                             +tbar(i,j,k,bi,bj))/float(nyears)       $                             +tbar(i,j,k,bi,bj))/float(nyears)
138                             else                             else
139                            tbar_gen(i,j,k,bi,bj) =tbar_gen(i,j,k,bi,bj)                             tbar_gen(i,j,k,bi,bj) =tbar_gen(i,j,k,bi,bj)
140       $                             +tbar(i,j,k,bi,bj)       $                             +tbar(i,j,k,bi,bj)
141                             endif                             endif
142                          enddo                          enddo
# Line 170  c--     Read time averages and the month Line 144  c--     Read time averages and the month
144                    enddo                    enddo
145                 enddo                 enddo
146              enddo              enddo
147           enddo              enddo
   
148  #else  #else
   
149  c--   Loop over records.  c--   Loop over records.
150        do irec = 1,nmonsrec        do irec = 1,nmonsrec
151    
152            irectmp = irec
153  c--     Read time averages and the monthly mean data.  c--     Read time averages and the monthly mean data.
154          call active_read_xyz( fnametheta, tbar, irec,          call active_read_xyz( fnametheta, tbar, irec,
155       &                        doglobalread, ladinit,       &                        doglobalread, ladinit,
156       &                        optimcycle, mythid,       &                        optimcycle, mythid,
157       &                        xx_tbar_mean_dummy )       &                        xx_tbar_mean_dummy )
   
158  #endif  #endif
   
159  c--     Determine the month to be read.  c--     Determine the month to be read.
160          levoff = mod(modelstartdate(1)/100,100)          levoff = mod(modelstartdate(1)/100,100)
161          levmon = (irec-1) + levoff          levmon = (irectmp-1) + levoff
162          levmon = mod(levmon-1,12)+1          levmon = mod(levmon-1,12)+1
163    
         print *, 'ph-rec irec, levmon ', irec, levmon  
164          call mdsreadfield( tdatfile, cost_iprec, cost_yftype,          call mdsreadfield( tdatfile, cost_iprec, cost_yftype,
165       &                     nr, tdat, levmon, mythid)       &                     nr, tdat, levmon, mythid)
166    
# Line 204  c--         Loop over the model layers Line 174  c--         Loop over the model layers
174  c--           Determine the mask on weights  c--           Determine the mask on weights
175                do j = jmin,jmax                do j = jmin,jmax
176                  do i = imin,imax                  do i = imin,imax
177                    cmask(i,j) = 1. _d 0                    cmask(i,j) = cosphi(i,j,bi,bj)
178                    if (tdat(i,j,k,bi,bj) .eq. 0.) then                    if (tdat(i,j,k,bi,bj) .eq. 0.) then
179                      cmask(i,j) = 0. _d 0                      cmask(i,j) = 0. _d 0
180                    endif                    else if (tdat(i,j,k,bi,bj) .lt. spval) then
                   
                   if (tdat(i,j,k,bi,bj) .lt. spval) then  
181                      cmask(i,j) = 0. _d 0                      cmask(i,j) = 0. _d 0
182                    endif                    else if (tdat(i,j,k,bi,bj) .gt. spmax) then
                     
 cph(  
 cph               print *, 'WARNING: SPECIFIC SETUP FOR ECCO'  
 cph               below statement could be replaced by following  
 cph               to make it independnet of Nr:  
 cph  
 cph               if ( rC(K) .GT. -1000. ) then  
 cph)  
 c                 set cmask=0 in areas shallower than 1000m  
                   if (_hFacC(i,j,13,bi,bj) .eq. 0.) then  
183                      cmask(i,j) = 0. _d 0                      cmask(i,j) = 0. _d 0
184                    endif                    endif
185                  enddo                  enddo
# Line 231  c--           Compute model data misfit Line 189  c--           Compute model data misfit
189  c             the temperature field.  c             the temperature field.
190                do j = jmin,jmax                do j = jmin,jmax
191                  do i = imin,imax                  do i = imin,imax
192                    if (_hFacC(i,j,k,bi,bj) .ne. 0.) then                    if ( _hFacC(i,j,k,bi,bj) .ne. 0. ) then
193                       fctile = fctile +                       fctile = fctile +
194       &                    (wtheta(k,bi,bj)*cosphi(i,j,bi,bj)*cmask(i,j)*       &                    (wthetaLev(i,j,k,bi,bj)*cmask(i,j)*
195         &                    (tbar(i,j,k,bi,bj) - tdat(i,j,k,bi,bj))*
196         &                    (tbar(i,j,k,bi,bj) - tdat(i,j,k,bi,bj)) )
197                         if ( wthetaLev(i,j,k,bi,bj)*cmask(i,j) .ne. 0. )
198         &                    num_temp(bi,bj) = num_temp(bi,bj) + 1. _d 0
199                         diagnosfld3d(i,j,k,bi,bj) =
200         &                    (wthetaLev(i,j,k,bi,bj)*cmask(i,j)*
201       &                    (tbar(i,j,k,bi,bj) - tdat(i,j,k,bi,bj))*       &                    (tbar(i,j,k,bi,bj) - tdat(i,j,k,bi,bj))*
202       &                    (tbar(i,j,k,bi,bj) - tdat(i,j,k,bi,bj))      )       &                    (tbar(i,j,k,bi,bj) - tdat(i,j,k,bi,bj)) )
203                      else
204                         diagnosfld3d(i,j,k,bi,bj) = 0.
205                    endif                    endif
206                  enddo                  enddo
207                enddo                enddo
# Line 243  c             the temperature field. Line 209  c             the temperature field.
209              enddo              enddo
210  c--         End of loop over layers.  c--         End of loop over layers.
211    
212                call mdswritefield( 'DiagnosCost_ClimTheta',
213         &           writeBinaryPrec, globalfiles, 'RL', Nr,
214         &           diagnosfld3d, irec, optimcycle, mythid )
215    
216              fcthread         = fcthread         + fctile              fcthread         = fcthread         + fctile
217              objf_temp(bi,bj) = objf_temp(bi,bj) + fctile              objf_temp(bi,bj) = objf_temp(bi,bj) + fctile
218    
 #ifdef ECCO_VERBOSE  
 c--         Print cost function for each tile in each thread.  
             write(msgbuf,'(a)') ' '  
             call print_message( msgbuf, standardmessageunit,  
      &                          SQUEEZE_RIGHT , mythid)  
             write(msgbuf,'(a,i8.8,1x,i3.3,1x,i3.3)')  
      &        ' cost_Theta: irec,bi,bj          =  ',irec,bi,bj  
             call print_message( msgbuf, standardmessageunit,  
      &                          SQUEEZE_RIGHT , mythid)  
             write(msgbuf,'(a,d22.15)')  
      &        '     cost function (temperature) = ',  
      &        fctile  
             call print_message( msgbuf, standardmessageunit,  
      &                      SQUEEZE_RIGHT , mythid)  
             write(msgbuf,'(a)') ' '  
             call print_message( msgbuf, standardmessageunit,  
      &                          SQUEEZE_RIGHT , mythid)  
 #endif  
   
219            enddo            enddo
220          enddo          enddo
221    
 #ifdef ECCO_VERBOSE  
 c--     Print cost function for all tiles.  
         _GLOBAL_SUM_R8( fcthread , myThid )  
         write(msgbuf,'(a)') ' '  
         call print_message( msgbuf, standardmessageunit,  
      &                      SQUEEZE_RIGHT , mythid)  
         write(msgbuf,'(a,i8.8)')  
      &    ' cost_Theta: irec = ',irec  
         call print_message( msgbuf, standardmessageunit,  
      &                      SQUEEZE_RIGHT , mythid)  
         write(msgbuf,'(a,a,d22.15)')  
      &    ' global cost function value',  
      &    ' (temperature) = ',fcthread  
         call print_message( msgbuf, standardmessageunit,  
      &                      SQUEEZE_RIGHT , mythid)  
         write(msgbuf,'(a)') ' '  
         call print_message( msgbuf, standardmessageunit,  
      &                      SQUEEZE_RIGHT , mythid)  
 #endif  
   
 #ifdef GENERIC_BAR_MONTH  
         endif  
 #endif  
   
222        enddo        enddo
223  c--   End of loop over records.  c--   End of loop over records.
224    
Legend:
Removed from v.1.2  
changed lines
  Added in v.1.8
  ViewVC Help
Powered by ViewVC 1.1.22