nc_rar.F 14.3 KB
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 386 387 388 389 390 391 392 393 394 395 396 397 398 399 400 401 402 403 404 405 406 407 408 409 410 411 412 413 414 415 416 417 418 419 420 421 422 423 424 425 426 427 428 429 430 431 432 433 434 435 436 437 438 439 440 441 442 443 444 445 446 447 448 449 450 451
*
* nc_rar.F
*
* James Caveen - UQAR
* Janvier 2007
*
* nc_rar.F: Reconstruire les axes reduits
**          Decompresse un champ de points mouilles et
*           reconstruit le champ 2D ou 3D correspondant
*
*           Reduced Axes Reconstruction
*           Rebuild a 2D or 3D field from data that has been compressed
*           using the Compression by gathering algorithm
*
*           See (Voir): http://www.cgd.ucar.edu/cms/eaton/cf-metadata/CF-current.html#gath
*                       for reference
*


* $Id: nc_rar.F,v 1.5 2009-06-09 17:25:15 caveenj Exp $
*
* $Log: not supported by cvs2svn $
* Revision 1.4.2.4  2009/06/09 17:20:11  caveenj
* Remplacement de -111 par ef_unspecified_int4 pour verifier si un axe est defini ou non
*
* Revision 1.4.2.3  2007/12/05 15:14:27  caveenj
* Correction de coquilles et ajout de documentation en anglais
*
* Revision 1.4.2.2  2007/11/28 18:50:47  caveenj
* Correction d'un bogue dans nc_rar pour le calcul des dimx et dimy
* On utilise maintenant la fonction ef_get_arg_extremes pour obtenir
* l'etendue complete du domaine.
* Ajout de code permettamt de traiter des variables 2D avec un indexwet2d
* et des variables 3D avec un indexwet3d.
*
* Ajout de la fonction rom_flip qui retourne un champ selon
* l'axe des Y (IROM) afin de faciliter l'affichage.
*
* Revision 1.4.2.1  2007/08/31 19:43:26  caveenj
* Nouvelle version de nc_rar qui utiise une variable au lieu d'un axe
* pour conserver les coordonnees composites des points de grille
* Cette nouvelle version elimine le besion de definir un axe de type
* spacing uneven pour les points mouilles
*
* Revision 1.3  2007/02/14 14:43:34  caveenj
* Correction d'un bogue de calcul d'indices par l'ajout d'un test permettant
* de s'assurer que les limites min-max des axes du tableau de sortie soient
* respectees avant d'y ajouter une valeur.
*
*
*
*
* In this subroutine we provide information about
* the function.  The user configurable information 
* consists of the following:
*
* descr              Text description of the function
*
* num_args           Required number of arguments
*
* axis_inheritance   Type of axis for the result
*                       ( CUSTOM, IMPLIED_BY_ARGS, NORMAL, ABSTRACT )
*                       CUSTOM          - user defined axis
*                       IMPLIED_BY_ARGS - same axis as the incoming argument
*                       NORMAL          - the result is normal to this axis
*                       ABSTRACT        - an axis which only has index values
*
* piecemeal_ok       For memory optimization:
*                       axes where calculation may be performed piecemeal
*                       ( YES, NO )
* 
*
* For each argument we provide the following information:
*
* name               Text name for an argument
*
* unit               Text units for an argument
*
* desc               Text description of an argument
*
* axis_influence     Are this argument's axes the same as the result grid?
*                       ( YES, NO )
*
* axis_extend       How much does Ferret need to extend arg limits relative to result 
*


      SUBROUTINE nc_rar_init(id)

      INCLUDE 'ferret_cmn/EF_Util.cmn'

      INTEGER id, arg

************************************************************************
*                                            USER CONFIGURABLE PORTION |
*                                                                      |
*                                                                      V

      CALL ef_set_desc(id,'Reconstruire des Axes Reduits'//
     $     ' - Reduced Axes Reconstruction' )

      CALL ef_set_num_args(id, 3)
      CALL ef_set_axis_inheritance(id, IMPLIED_BY_ARGS, 
     .     IMPLIED_BY_ARGS,IMPLIED_BY_ARGS, IMPLIED_BY_ARGS)
      CALL ef_set_piecemeal_ok(id, NO, NO, NO, NO)

      arg = 1
      CALL ef_set_arg_name(id, arg, 'A')
      CALL ef_set_arg_desc(id, arg,'Champ decrivant la grille de '//
     $      'sortie - Field describing output grid')
      CALL ef_set_axis_influence(id, arg, YES, YES, YES, YES)

      arg = 2
      CALL ef_set_arg_name(id, arg, 'B')
      call ef_set_arg_desc(id,arg,'Champ a reconstruire - ' //
     $     'Field to rebuild')
      CALL ef_set_axis_influence(id, arg, NO, NO, NO, YES)


      arg = 3
      CALL ef_set_arg_name(id, arg, 'C')
      call ef_set_arg_desc(id,arg,'Champ des indices pts mouilles - ' //
     $     'Wet point indices field')
      CALL ef_set_axis_influence(id, arg, NO, NO, NO, NO)


*                                                                      ^
*                                                                      |
*                                            USER CONFIGURABLE PORTION |
************************************************************************

      RETURN 
      END




*
* Fonction calculant le resultat
*
* In this subroutine we compute the result
*
*
*Parameters:
*           arg_1: dummy variable used to describe the target (2)3D grig
*           arg_2: 2d or 3d field to uncompress
*           arg_3: field containing the composite indices of the 2D 3D field
*           result: the resulting (2)3D field


      SUBROUTINE nc_rar_compute(id, arg_1, arg_2, arg_3, result)

      implicit none

      INCLUDE 'ferret_cmn/EF_Util.cmn'
      INCLUDE 'ferret_cmn/EF_mem_subsc.cmn'

      INTEGER id

      REAL bad_flag(EF_MAX_ARGS), bad_flag_result

      REAL arg_1(mem1lox:mem1hix, mem1loy:mem1hiy, 
     .     mem1loz:mem1hiz, mem1lot:mem1hit)

      REAL arg_2(mem2lox:mem2hix, mem2loy:mem2hiy, 
     .     mem2loz:mem2hiz, mem2lot:mem2hit)


      REAL arg_3(mem3lox:mem3hix, mem3loy:mem3hiy, 
     .     mem3loz:mem3hiz, mem3lot:mem3hit)

      REAL result(memreslox:memreshix, memresloy:memreshiy, 
     .     memresloz:memreshiz, memreslot:memreshit)

* After initialization, the 'res_' arrays contain indexing information 
* for the result axes.  The 'arg_' arrays will contain the indexing 
* information for each variable's axes. 

      INTEGER res_lo_ss(4), res_hi_ss(4), res_incr(4)
      INTEGER arg_lo_ss(4,EF_MAX_ARGS), arg_hi_ss(4,EF_MAX_ARGS),
     $          arg_incr(4,EF_MAX_ARGS)

      INTEGER ss_min(4,EF_MAX_ARGS), ss_max(4,EF_MAX_ARGS)



************************************************************************
*                                            USER CONFIGURABLE PORTION |
*                                                                      |
*                                                                      V

      INTEGER i,j,k,l
      INTEGER i2, j2, k2, l2
      INTEGER i3, j3, k3, l3

      integer arg
      character*100 errtxt

      integer dimxy,dimx,indx,indy,indz,indice,dimy

     
      CALL ef_get_res_subscripts(id, res_lo_ss, res_hi_ss, res_incr)
      CALL ef_get_arg_subscripts(id, arg_lo_ss, arg_hi_ss, arg_incr)
      CALL ef_get_bad_flags(id, bad_flag, bad_flag_result)

C
C     Ici on va chercher les extremes de TOUT le domaine
C     afin de pouvoir calculer dimx et dimy qui servent a recalculer
C     les indices i,j,k a partir de l'indice composite de indexwet
C
C     Here we fetch the output domain extremes in order to
C     calculate dimx and dimy wich are used in the calculation
C     of the i,j,k indices from the "composite" wet index
      call ef_get_arg_ss_extremes(id,1,ss_min,ss_max)




#if defined DEBUG
      write(6,*)'extremes'

      write(6,*)'ss_min,ss_max X:',ss_min(X_AXIS,1),ss_max(X_AXIS,1)
      write(6,*)'ss_min,ss_max Y:',ss_min(Y_AXIS,1),ss_max(Y_AXIS,1)
      write(6,*)'ss_min ss_max Z',ss_min(Z_AXIS,1),ss_max(Z_AXIS,1)


      write(6,*)'res_lo_ss(X_AXIS), res_hi_ss(X_AXIS)',
     &     res_lo_ss(X_AXIS), res_hi_ss(X_AXIS)

      write(6,*)'res_lo_ss(Y_AXIS), res_hi_ss(Y_AXIS)',
     &     res_lo_ss(Y_AXIS), res_hi_ss(Y_AXIS)

      write(6,*)'res_lo_ss(Z_AXIS), res_hi_ss(Z_AXIS)',
     &     res_lo_ss(Z_AXIS), res_hi_ss(Z_AXIS)

      write(6,*)'res_lo_ss(T_AXIS), res_hi_ss(T_AXIS)',
     &     res_lo_ss(T_AXIS), res_hi_ss(T_AXIS)


 
*
*Calcul de l'indice de depart j et k du champ compresse
*cet indice demeure constant dans le calcul
* et est probablement = ef_unspecified_int4 ( i.e., non defini)
        j2 = arg_lo_ss(Y_AXIS,ARG2)
        k2 = arg_lo_ss(Z_AXIS,ARG2)
        l2 = arg_lo_ss(T_AXIS,ARG2)

        j3 = arg_lo_ss(Y_AXIS,ARG3)
        k3 = arg_lo_ss(Z_AXIS,ARG3)
        l3 = arg_lo_ss(T_AXIS,ARG3)



       write(6,*)'mem1lox:mem1hix, mem1loy:mem1hiy,'// 
     &     'mem1loz:mem1hiz, mem1lot:mem1hit',
     &     mem1lox,mem1hix, mem1loy,mem1hiy, 
     &     mem1loz,mem1hiz, mem1lot,mem1hit


       i2=arg_lo_ss(X_AXIS,ARG2)
        do i = 1,12
           write(6,*)'arg_2:',arg_2(i2,j2,k2,l2)
            i2 = i2 + arg_incr(X_AXIS,ARG2)
         enddo

       i3=arg_lo_ss(X_AXIS,ARG3)
         do i = 1,12
            write(6,*)'arg_3:',arg_3(i3,j3,k3,l3)
            i3 = i3 + arg_incr(X_AXIS,ARG3)
         enddo

        write(6,*)'arg_lo_ss(X_AXIS,ARG1), arg_hi_ss(X_AXIS,ARG1),
     &             arg_incr(X_AXIS,ARG1)',arg_lo_ss(X_AXIS,ARG1), 
     &             arg_hi_ss(X_AXIS,ARG1), arg_incr(X_AXIS,ARG1)

        write(6,*)'arg_lo_ss(Y_AXIS,ARG1), arg_hi_ss(Y_AXIS,ARG1),
     &             arg_incr(Y_AXIS,ARG1)',arg_lo_ss(Y_AXIS,ARG1), 
     &             arg_hi_ss(Y_AXIS,ARG1), arg_incr(Y_AXIS,ARG1)

        write(6,*)'arg_lo_ss(Z_AXIS,ARG1), arg_hi_ss(Z_AXIS,ARG1),
     &             arg_incr(Z_AXIS,ARG1)',arg_lo_ss(Z_AXIS,ARG1), 
     &             arg_hi_ss(Z_AXIS,ARG1), arg_incr(Z_AXIS,ARG1)

        write(6,*)'arg_lo_ss(T_AXIS,ARG1), arg_hi_ss(T_AXIS,ARG1),
     &             arg_incr(T_AXIS,ARG1)',arg_lo_ss(T_AXIS,ARG1), 
     &             arg_hi_ss(T_AXIS,ARG1), arg_incr(T_AXIS,ARG1)

#endif
C
C     Initialiser le champ de sortie a des manquants
C     Initialise output field to missing
C

      do  i=res_lo_ss(X_AXIS), res_hi_ss(X_AXIS)
         do  j=res_lo_ss(Y_AXIS), res_hi_ss(Y_AXIS)
            do  k=res_lo_ss(Z_AXIS), res_hi_ss(Z_AXIS)
               do  l=res_lo_ss(T_AXIS), res_hi_ss(T_AXIS)
                  result(i,j,k,l) = bad_flag_result
               enddo
            enddo
         enddo
      enddo




C
C    Reconstruction des axes de la grille 3D
C    et des valeurs du champ sur cette grille
C 
C    Rebuild axes for (2)3D grid and scatter values
C    of the field unto the grid   
C



*
* Calcul des dimensions necessaires au calcul des coordonnees
* Compute  constants required to decompress
*

        dimx = ss_max(X_AXIS,1) - ss_min(X_AXIS,1) + 1
        dimy = ss_max(Y_AXIS,1) - ss_min(Y_AXIS,1) + 1
        dimxy = dimx * dimy



*
*Calcul de l'indice de depart j et k du champ compresse
*cet indice demeure constant dans le calcul
* et est probablement = ef_unspecified_int4 ( i.e., non defini)

*
* Compute the starting index for j and k (all data is stored
* along x and t in compressed format)
* Theses indices remain constant throughout calculation and are
* most likely == ef_unspecified_int4 (i.e., undefined)

        j2 = arg_lo_ss(Y_AXIS,ARG2)
        k2 = arg_lo_ss(Z_AXIS,ARG2)
        l2 = arg_lo_ss(T_AXIS,ARG2)

        j3 = arg_lo_ss(Y_AXIS,ARG3)
        k3 = arg_lo_ss(Z_AXIS,ARG3)
        l3 = arg_lo_ss(T_AXIS,ARG3)


* Ici, on boucle sur les coordonnes des points mouilles
* et on recalcule les coordonnes i,j,k du champ 3D
* a partir de la coordonnee 1D 
* puis on assigne a result(i,j,k) la valeur
* correspondante dans le champ qui avait ete cree lors  
* de la compression des axes (arg_2)
*
* Here, we loop on the compressed field coordinates in order to 
* compute the i,j,k coordinates of the 3D field
* from the 1D coordinates
* We then assing to result(i,j,k) the corresponding value in the
* compressed field
*
*       3D 
        if(ss_min(Z_AXIS,1) .ne. ef_unspecified_int4 ) then
           i2 = arg_lo_ss(X_AXIS,ARG2)
           do i3 = arg_lo_ss(X_AXIS,ARG2),arg_hi_ss(X_AXIS,ARG2)
              indice = int( arg_3( i3, j3, k3, l3) )
              indz = (indice - 1)/dimxy
              indice = indice - indz*dimxy
              indy = (indice - 1)/dimx
              indx = indice - indy*dimx
              indz = indz + 1 
              indy = indy + 1
*     
*     On s'assure que les indices calcules sont a
*     l'interieur du domaine demande 
*     
*     
*     Making sure that the computed indices fall within
*     the requested region
*     
              if( (indx .ge. res_lo_ss(X_AXIS) .and. 
     &             indx .le. res_hi_ss(X_AXIS)) .and.
     &             (indy .ge. res_lo_ss(Y_AXIS) .and. 
     &             indy .le. res_hi_ss(Y_AXIS)) .and.
     &             (indz .ge. res_lo_ss(Z_AXIS) .and. 
     &             indz .le. res_hi_ss(Z_AXIS))) then

                 
*     
*     On repete l'operation pour chaque pas de temps demande
*     Repeat for all requested time step
*     
                 l2 = arg_lo_ss(T_AXIS,ARG2) 

                 do l = res_lo_ss(T_AXIS),res_hi_ss(T_AXIS)
                    result(indx,indy,indz,l) = arg_2(i2,j2,k2,l2)
                    l2 = l2 + arg_incr(T_AXIS,ARG2)
                 enddo
              endif
              i2 = i2 + arg_incr(X_AXIS,ARG2)
           enddo

C     2D
        else

           indz = ef_unspecified_int4
           i2 = arg_lo_ss(X_AXIS,ARG2)
           do i3 = arg_lo_ss(X_AXIS,ARG2),arg_hi_ss(X_AXIS,ARG2)
              indice = int( arg_3( i3, j3, k3, l3) )
              indy = (indice - 1)/dimx
              indx = indice - indy*dimx
              indy = indy + 1
*     
*     On s'assure que les indices calcules sont a
*     l'interieur du domaine demande 
*     
*     
*     Making sure that the computed indices fall within
*     the requested region
*     
              if( (indx .ge. res_lo_ss(X_AXIS) .and. 
     &             indx .le. res_hi_ss(X_AXIS)) .and.
     &             (indy .ge. res_lo_ss(Y_AXIS) .and. 
     &             indy .le. res_hi_ss(Y_AXIS)) ) then

                 
*     
*     On repete l'operation pour chaque pas de temps demande
*     Repeat for all requested time step
*     
                 l2 = arg_lo_ss(T_AXIS,ARG2) 

                 do l = res_lo_ss(T_AXIS),res_hi_ss(T_AXIS)
                    result(indx,indy,indz,l) = arg_2(i2,j2,k2,l2)
                    l2 = l2 + arg_incr(T_AXIS,ARG2)
                 enddo
              endif
              i2 = i2 + arg_incr(X_AXIS,ARG2)
           enddo

        endif


      
*                                                                      ^
*                                                                      |
*                                            USER CONFIGURABLE PORTION |
************************************************************************

      RETURN 
      END