Ignore:
Timestamp:
2011-09-15T08:41:58+02:00 (9 years ago)
Author:
cbricaud
Message:

modified routine for netcdf output

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/2011/dev_r2802_MERCATOR9_floats/NEMOGCM/NEMO/OPA_SRC/FLO/flodom.F90

    r2528 r2839  
    4343      !!               the longitude (degree) and the depth (m). 
    4444      !!----------------------------------------------------------------------       
    45       LOGICAL  ::   llinmesh 
    46       INTEGER  ::   ji, jj, jk   ! DO loop index on 3 directions 
    47       INTEGER  ::   jfl, jfl1    ! number of floats    
    48       INTEGER  ::   inum         ! logical unit for file read 
     45      CHARACTER (len=21) ::  clname        ! floats initialisation filename 
     46      LOGICAL            ::   llinmesh 
     47      INTEGER            ::   ji, jj, jk   ! DO loop index on 3 directions 
     48      INTEGER            ::   jfl, jfl1    ! number of floats    
     49      INTEGER            ::   inum         ! logical unit for file read 
     50      INTEGER            ::   jtrash       ! trash var for reading   
     51      INTEGER            ::   ierr 
    4952      INTEGER, DIMENSION(jpnfl) ::   iimfl, ijmfl, ikmfl       ! index mesh of floats 
    5053      INTEGER, DIMENSION(jpnfl) ::   idomfl,  ivtest, ihtest   !   -             - 
     
    6669 
    6770         ! read of the restart file 
    68          READ(inum) ( tpifl  (jfl), jfl=1, jpnrstflo),   &  
     71         READ(inum,*)  ( tpifl  (jfl), jfl=1, jpnrstflo),   &  
    6972                        ( tpjfl  (jfl), jfl=1, jpnrstflo),   & 
    7073                        ( tpkfl  (jfl), jfl=1, jpnrstflo),   & 
     
    208211         ENDIF 
    209212      ELSE 
    210          IF(lwp) WRITE(numout,*) '                     init_float read ' 
     213 
     214         IF( ln_ariane )THEN 
     215 
     216            IF(lwp) WRITE(numout,*) '                     init_float read with ariane convention (mesh indexes)' 
     217 
     218            ! First initialisation of floats with ariane convention 
     219            !  
     220            ! The indexes are read directly from file (warning ariane 
     221            ! convention, are refered to  
     222            ! U,V,W grids - and not T-)  
     223            ! The isobar advection is managed with the sign of tpkfl ( >0 -> 3D 
     224            ! advection, <0 -> 2D)  
     225            ! Some variables are not read, as - gl         : time index; 4th 
     226            ! column         
     227            !                                 - transport  : transport ; 5th 
     228            !                                 column 
     229            ! and paste in the jtrash var 
     230            ! At the end, ones need to replace the indexes on T grid 
     231            ! RMQ : there is no float groups identification ! 
     232  
     233            clname='init_float_ariane' 
     234 
     235            nisobfl = 1 ! we assume that by default we want 3D advection 
     236             
     237            ! we check that the number of floats in the init_file are consistant 
     238            ! with the namelist 
     239            IF( lwp ) THEN  
     240               jfl1=0 
     241               OPEN( unit=inum, file=clname,status='old',access='sequential',form='formatted') 
     242               DO WHILE (ierr .GE. 0) 
     243                 jfl1=jfl1+1 
     244                 READ (inum,*, iostat=ierr) 
     245               END DO 
     246               CLOSE(inum) 
     247               IF( (jfl1-1) .NE. jpnfl )THEN 
     248                  WRITE (numout,*) ' STOP the number of floats in' ,clname,'  = ',jfl1 
     249                  WRITE (numout,*) '  is not equal to jfl= ',jpnfl  
     250                  STOP 
     251               ENDIF  
     252            ENDIF  
     253 
     254            ! we get the init values  
     255            CALL ctl_opn( inum, clname, 'OLD', 'FORMATTED', 'SEQUENTIAL',   & 
     256               &         1, numout, .TRUE., 1 ) 
     257            DO jfl = 1, jpnfl 
     258              READ(inum,*) tpifl(jfl),tpjfl(jfl),tpkfl(jfl),jtrash, jtrash 
     259              if(lwp)write(numout,*)"read : ",tpifl(jfl),tpjfl(jfl),tpkfl(jfl),jtrash, jtrash ; call flush(numout) 
     260  
     261              IF ( tpkfl(jfl) .LT. 0. ) nisobfl(jfl) = 0 !set the 2D advection according to init_float 
     262              ngrpfl(jfl)=jfl 
     263            END DO 
     264 
     265            ! conversion from ariane index to T grid index 
     266            tpkfl = abs(tpkfl)-0.5 ! reversed vertical axis 
     267            tpifl = tpifl+0.5  
     268            tpjfl = tpjfl+0.5 
     269 
     270            ! verif of non land point initialisation : no need if correct init 
     271             
     272         ELSE  
     273            IF(lwp) WRITE(numout,*) '                     init_float read ' 
    211274          
    212          ! First initialisation of floats 
    213          ! the initials positions of floats are written in a file 
    214          ! with a variable to know if it is a isobar float a number  
    215          ! to identified who want the trajectories of this float and  
    216          ! an index for the number of the float          
    217          ! open the init file  
    218          CALL ctl_opn( inum, 'init_float', 'OLD', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp ) 
    219          READ(inum) (flxx(jfl)   , jfl=1, jpnfl),   & 
    220                     (flyy(jfl)   , jfl=1, jpnfl),   & 
    221                     (flzz(jfl)   , jfl=1, jpnfl),   & 
    222                     (nisobfl(jfl), jfl=1, jpnfl),   & 
    223                     (ngrpfl(jfl) , jfl=1, jpnfl) 
    224          CLOSE(inum) 
    225              
    226          ! Test to find the grid point coordonate with the geographical position          
    227          DO jfl = 1, jpnfl 
    228             ihtest(jfl) = 0 
    229             ivtest(jfl) = 0 
    230             ikmfl(jfl) = 0 
     275            ! First initialisation of floats 
     276            ! the initials positions of floats are written in a file 
     277            ! with a variable to know if it is a isobar float a number  
     278            ! to identified who want the trajectories of this float and  
     279            ! an index for the number of the float          
     280            ! open the init file  
     281            CALL ctl_opn( inum, 'init_float', 'OLD', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp ) 
     282            READ(inum,*) (flxx(jfl)   , jfl=1, jpnfl),   & 
     283                         (flyy(jfl)   , jfl=1, jpnfl),   & 
     284                         (flzz(jfl)   , jfl=1, jpnfl),   & 
     285                         (nisobfl(jfl), jfl=1, jpnfl),   & 
     286                         (ngrpfl(jfl) , jfl=1, jpnfl) 
     287            CLOSE(inum) 
     288             
     289            ! Test to find the grid point coordonate with the geographical position          
     290            DO jfl = 1, jpnfl 
     291              ihtest(jfl) = 0 
     292              ivtest(jfl) = 0 
     293              ikmfl(jfl) = 0 
    231294# if   defined key_mpp_mpi 
    232             DO ji = MAX(nldi,2), nlei 
    233                DO jj = MAX(nldj,2), nlej   ! NO vector opt. 
    234 # else 
    235             DO ji = 2, jpi 
    236                DO jj = 2, jpj   ! NO vector opt. 
     295               DO ji = MAX(nldi,2), nlei 
     296                  DO jj = MAX(nldj,2), nlej   ! NO vector opt. 
     297# else  
     298               DO ji = 2, jpi 
     299                  DO jj = 2, jpj   ! NO vector opt. 
    237300# endif                   
    238                   ! for each float we find the indexes of the mesh  
     301                     ! for each float we find the indexes of the mesh  
    239302                   
    240                   CALL findmesh(glamf(ji-1,jj-1),gphif(ji-1,jj-1),   & 
    241                                 glamf(ji-1,jj  ),gphif(ji-1,jj  ),   & 
    242                                 glamf(ji  ,jj  ),gphif(ji  ,jj  ),   & 
    243                                 glamf(ji  ,jj-1),gphif(ji  ,jj-1),   & 
    244                                 flxx(jfl)       ,flyy(jfl)       ,   & 
    245                                 glamt(ji  ,jj  ),gphit(ji  ,jj  ), llinmesh) 
    246                   IF(llinmesh) THEN 
    247                      iimfl(jfl)  = ji 
    248                      ijmfl(jfl)  = jj 
    249                      ihtest(jfl) = ihtest(jfl)+1 
    250                      DO jk = 1, jpk-1 
    251                         IF( (fsdepw(ji,jj,jk) <= flzz(jfl)) .AND. (fsdepw(ji,jj,jk+1) >  flzz(jfl)) ) THEN 
    252                            ikmfl(jfl)  = jk 
    253                            ivtest(jfl) = ivtest(jfl) + 1 
    254                         ENDIF 
    255                      END DO 
    256                   ENDIF 
     303                     CALL findmesh(glamf(ji-1,jj-1),gphif(ji-1,jj-1),   & 
     304                                   glamf(ji-1,jj  ),gphif(ji-1,jj  ),   & 
     305                                   glamf(ji  ,jj  ),gphif(ji  ,jj  ),   & 
     306                                   glamf(ji  ,jj-1),gphif(ji  ,jj-1),   & 
     307                                   flxx(jfl)       ,flyy(jfl)       ,   & 
     308                                   glamt(ji  ,jj  ),gphit(ji  ,jj  ), llinmesh) 
     309                     IF(llinmesh) THEN 
     310                        iimfl(jfl)  = ji 
     311                        ijmfl(jfl)  = jj 
     312                        ihtest(jfl) = ihtest(jfl)+1 
     313                        DO jk = 1, jpk-1 
     314                           IF( (fsdepw(ji,jj,jk) <= flzz(jfl)) .AND. (fsdepw(ji,jj,jk+1) >  flzz(jfl)) ) THEN 
     315                              ikmfl(jfl)  = jk 
     316                              ivtest(jfl) = ivtest(jfl) + 1 
     317                           ENDIF 
     318                        END DO 
     319                     ENDIF 
     320                  END DO 
    257321               END DO 
    258             END DO 
    259              
    260             ! If the float is in a mesh computed by an other processor we put iimfl=ijmfl=-1             
    261             IF( ihtest(jfl) == 0 ) THEN 
    262                iimfl(jfl) = -1 
    263                ijmfl(jfl) = -1 
    264             ENDIF 
    265          END DO 
     322             
     323               ! If the float is in a mesh computed by an other processor we put iimfl=ijmfl=-1             
     324               IF( ihtest(jfl) == 0 ) THEN 
     325                  iimfl(jfl) = -1 
     326                  ijmfl(jfl) = -1 
     327               ENDIF 
     328            END DO 
    266329          
    267          ! A zero in the sum of the arrays "ihtest" and "ivtest"           
    268          IF( lk_mpp )   CALL mpp_sum(ihtest,jpnfl)   ! sums over the global domain 
    269          IF( lk_mpp )   CALL mpp_sum(ivtest,jpnfl) 
    270  
    271          DO jfl = 1, jpnfl 
    272             IF( (ihtest(jfl) > 1 ) .OR. ( ivtest(jfl) > 1 )) THEN 
    273                IF(lwp) WRITE(numout,*) 'THE FLOAT',jfl,' IS NOT IN ONLY ONE MESH' 
    274             ENDIF 
    275             IF( ihtest(jfl) == 0 ) THEN  
    276                IF(lwp) WRITE(numout,*)'THE FLOAT',jfl,' IS IN NO MESH' 
    277             ENDIF 
    278          END DO 
     330            ! A zero in the sum of the arrays "ihtest" and "ivtest"           
     331            IF( lk_mpp )   CALL mpp_sum(ihtest,jpnfl)   ! sums over the global domain 
     332            IF( lk_mpp )   CALL mpp_sum(ivtest,jpnfl) 
     333 
     334            DO jfl = 1, jpnfl 
     335               IF( (ihtest(jfl) > 1 ) .OR. ( ivtest(jfl) > 1 )) THEN 
     336                  IF(lwp) WRITE(numout,*) 'THE FLOAT',jfl,' IS NOT IN ONLY ONE MESH' 
     337               ENDIF 
     338               IF( ihtest(jfl) == 0 ) THEN  
     339                  IF(lwp) WRITE(numout,*)'THE FLOAT',jfl,' IS IN NO MESH' 
     340               ENDIF 
     341            END DO 
    279342         
    280          ! We compute the distance between the float and the face of  the mesh          
    281          DO jfl = 1, jpnfl 
    282             ! Made only if the float is in the domain of the processor 
    283             IF( (iimfl(jfl) >= 0 ) .AND. ( ijmfl(jfl) >= 0 ) ) THEN 
    284                 
    285                ! TEST TO KNOW IF THE FLOAT IS NOT INITIALISED IN THE COAST 
    286                 
    287                idomfl(jfl) = 0 
    288                IF( tmask(iimfl(jfl),ijmfl(jfl),ikmfl(jfl)) == 0. ) idomfl(jfl)=1 
    289                 
    290                ! Computation of the distance between the float 
    291                ! and the faces of the mesh 
    292                !            zdxab 
    293                !             . 
    294                !        B----.---------C 
    295                !        |    .         | 
    296                !        |<------>flo   | 
    297                !        |        ^     | 
    298                !        |        |.....|....zdyad 
    299                !        |        |     | 
    300                !        A--------|-----D 
    301                 
    302                zdxab = dstnce(flxx(jfl),flyy(jfl),glamf(iimfl(jfl)-1,ijmfl(jfl)-1),flyy(jfl))                 
    303                zdyad = dstnce(flxx(jfl),flyy(jfl),flxx(jfl),gphif(iimfl(jfl)-1,ijmfl(jfl)-1)) 
    304                 
    305                ! Translation of this distances (in meter) in indexes 
    306                 
    307                tpifl(jfl) = (iimfl(jfl)-0.5)+zdxab/ e1u(iimfl(jfl)-1,ijmfl(jfl))+(mig(1)-jpizoom) 
    308                tpjfl(jfl) = (ijmfl(jfl)-0.5)+zdyad/ e2v(iimfl(jfl),ijmfl(jfl)-1)+(mjg(1)-jpjzoom) 
    309                tpkfl(jfl) = (fsdepw(iimfl(jfl),ijmfl(jfl),ikmfl(jfl)+1) - flzz(jfl))*(ikmfl(jfl))                     & 
    310                           / (fsdepw(iimfl(jfl),ijmfl(jfl),ikmfl(jfl)+1) - fsdepw(iimfl(jfl),ijmfl(jfl),ikmfl(jfl)))   & 
    311                           + (flzz(jfl) - fsdepw(iimfl(jfl),ijmfl(jfl),ikmfl(jfl)))*(ikmfl(jfl)+1)                     & 
    312                           / (fsdepw(iimfl(jfl),ijmfl(jfl),ikmfl(jfl)+1) - fsdepw(iimfl(jfl),ijmfl(jfl),ikmfl(jfl))) 
    313             ELSE 
    314                tpifl (jfl) = 0.e0 
    315                tpjfl (jfl) = 0.e0 
    316                tpkfl (jfl) = 0.e0 
    317                idomfl(jfl) = 0 
    318             ENDIF 
    319          END DO 
     343            ! We compute the distance between the float and the face of  the mesh          
     344            DO jfl = 1, jpnfl 
     345               ! Made only if the float is in the domain of the processor 
     346               IF( (iimfl(jfl) >= 0 ) .AND. ( ijmfl(jfl) >= 0 ) ) THEN 
     347                
     348                  ! TEST TO KNOW IF THE FLOAT IS NOT INITIALISED IN THE COAST 
     349                
     350                  idomfl(jfl) = 0 
     351                  IF( tmask(iimfl(jfl),ijmfl(jfl),ikmfl(jfl)) == 0. ) idomfl(jfl)=1 
     352                
     353                  ! Computation of the distance between the float 
     354                  ! and the faces of the mesh 
     355                  !            zdxab 
     356                  !             . 
     357                  !        B----.---------C 
     358                  !        |    .         | 
     359                  !        |<------>flo   | 
     360                  !        |        ^     | 
     361                  !        |        |.....|....zdyad 
     362                  !        |        |     | 
     363                  !        A--------|-----D 
     364                
     365                  zdxab = dstnce(flxx(jfl),flyy(jfl),glamf(iimfl(jfl)-1,ijmfl(jfl)-1),flyy(jfl))                 
     366                  zdyad = dstnce(flxx(jfl),flyy(jfl),flxx(jfl),gphif(iimfl(jfl)-1,ijmfl(jfl)-1)) 
     367                
     368                  ! Translation of this distances (in meter) in indexes 
     369                
     370                  tpifl(jfl) = (iimfl(jfl)-0.5)+zdxab/ e1u(iimfl(jfl)-1,ijmfl(jfl))+(mig(1)-jpizoom) 
     371                  tpjfl(jfl) = (ijmfl(jfl)-0.5)+zdyad/ e2v(iimfl(jfl),ijmfl(jfl)-1)+(mjg(1)-jpjzoom) 
     372                  tpkfl(jfl) = (fsdepw(iimfl(jfl),ijmfl(jfl),ikmfl(jfl)+1) - flzz(jfl))*(ikmfl(jfl))                     & 
     373                             / (fsdepw(iimfl(jfl),ijmfl(jfl),ikmfl(jfl)+1) - fsdepw(iimfl(jfl),ijmfl(jfl),ikmfl(jfl)))   & 
     374                             + (flzz(jfl) - fsdepw(iimfl(jfl),ijmfl(jfl),ikmfl(jfl)))*(ikmfl(jfl)+1)                     & 
     375                             / (fsdepw(iimfl(jfl),ijmfl(jfl),ikmfl(jfl)+1) - fsdepw(iimfl(jfl),ijmfl(jfl),ikmfl(jfl))) 
     376               ELSE 
     377                  tpifl (jfl) = 0.e0 
     378                  tpjfl (jfl) = 0.e0 
     379                  tpkfl (jfl) = 0.e0 
     380                  idomfl(jfl) = 0 
     381               ENDIF 
     382            END DO 
    320383          
    321          ! The sum of all the arrays tpifl, tpjfl, tpkfl give 3 arrays with the positions of all the floats.  
    322          IF( lk_mpp )   CALL mpp_sum( tpifl , jpnfl )   ! sums over the global domain 
    323          IF( lk_mpp )   CALL mpp_sum( tpjfl , jpnfl ) 
    324          IF( lk_mpp )   CALL mpp_sum( tpkfl , jpnfl ) 
    325          IF( lk_mpp )   CALL mpp_sum( idomfl, jpnfl ) 
     384            ! The sum of all the arrays tpifl, tpjfl, tpkfl give 3 arrays with the positions of all the floats.  
     385            IF( lk_mpp )   CALL mpp_sum( tpifl , jpnfl )   ! sums over the global domain 
     386            IF( lk_mpp )   CALL mpp_sum( tpjfl , jpnfl ) 
     387            IF( lk_mpp )   CALL mpp_sum( tpkfl , jpnfl ) 
     388            IF( lk_mpp )   CALL mpp_sum( idomfl, jpnfl ) 
     389         ENDIF 
     390 
    326391      ENDIF 
    327392             
Note: See TracChangeset for help on using the changeset viewer.