New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
Changeset 3051 for branches/2011/dev_MERCATOR_2011_MERGE/NEMOGCM/NEMO/OPA_SRC/FLO/flodom.F90 – NEMO

Ignore:
Timestamp:
2011-11-07T16:03:29+01:00 (12 years ago)
Author:
cbricaud
Message:

add changes from dev_r2802_MERCATOR9_floats

File:
1 edited

Legend:

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

    r2528 r3051  
    44   !! Ocean floats :   domain 
    55   !!====================================================================== 
    6    !! History :  OPA  ! 1998-07 (Y.Drillet, CLIPPER)  Original code 
     6   !! History :  OPA          ! 1998-07 (Y.Drillet, CLIPPER)  Original code 
     7   !!            NEMO_3.3.1   ! 2011-09 (C.Bricaud,S.Law-Chune Mercator-Ocean):  
     8                              ! add Ariane convention, Comsecitc changes 
    79   !!---------------------------------------------------------------------- 
    810#if   defined key_floats   ||   defined key_esopa 
     
    1012   !!   'key_floats'                                     float trajectories 
    1113   !!---------------------------------------------------------------------- 
    12    !!   flo_dom        : initialization of floats 
    13    !!   findmesh       : compute index of position  
    14    !!   dstnce         : compute distance between face mesh and floats  
     14   !!   flo_dom               : initialization of floats 
     15   !!   add_new_floats        : add new floats (long/lat/depth) 
     16   !!   add_new_ariane_floats : add new floats with araine convention (i/j/k) 
     17   !!   findmesh              : compute index of position  
     18   !!   dstnce                : compute distance between face mesh and floats  
    1519   !!---------------------------------------------------------------------- 
    1620   USE oce             ! ocean dynamics and tracers 
     
    2327   PRIVATE 
    2428 
    25    PUBLIC   flo_dom    ! routine called by floats.F90 
     29   PUBLIC   flo_dom         ! routine called by floats.F90 
     30   PUBLIC   flo_dom_alloc   ! Routine called in floats.F90 
     31 
     32   CHARACTER (len=21) ::  clname1 = 'init_float'              ! floats initialisation filename 
     33   CHARACTER (len=21) ::  clname2 = 'init_float_ariane'       ! ariane floats initialisation filename 
     34 
     35 
     36   INTEGER , ALLOCATABLE, DIMENSION(:) ::   iimfl, ijmfl, ikmfl       ! index mesh of floats 
     37   INTEGER , ALLOCATABLE, DIMENSION(:) ::   idomfl, ivtest, ihtest    !   -      
     38   REAL(wp), ALLOCATABLE, DIMENSION(:) ::   zgifl, zgjfl,  zgkfl      ! distances in indexes 
    2639 
    2740   !! * Substitutions 
     
    4356      !!               the longitude (degree) and the depth (m). 
    4457      !!----------------------------------------------------------------------       
    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 
    49       INTEGER, DIMENSION(jpnfl) ::   iimfl, ijmfl, ikmfl       ! index mesh of floats 
    50       INTEGER, DIMENSION(jpnfl) ::   idomfl,  ivtest, ihtest   !   -             - 
    51       REAL(wp) ::   zdxab, zdyad 
    52       REAL(wp), DIMENSION(jpnnewflo+1)  :: zgifl, zgjfl,  zgkfl 
     58      INTEGER            ::   jfl    ! dummy loop   
     59      INTEGER            ::   inum   ! logical unit for file read 
    5360      !!--------------------------------------------------------------------- 
    5461       
     
    5966      IF(lwp) WRITE(numout,*) '           jpnfl = ',jpnfl 
    6067       
    61       IF(ln_rstflo) THEN 
     68      !-------------------------! 
     69      ! FLOAT RESTART FILE READ ! 
     70      !-------------------------! 
     71      IF( ln_rstflo )THEN 
     72 
    6273         IF(lwp) WRITE(numout,*) '        float restart file read' 
    6374          
    6475         ! open the restart file  
     76         !---------------------- 
    6577         CALL ctl_opn( inum, 'restart_float', 'OLD', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp ) 
    6678 
    6779         ! read of the restart file 
    68          READ(inum) ( tpifl  (jfl), jfl=1, jpnrstflo),   &  
     80         READ(inum,*)  ( tpifl  (jfl), jfl=1, jpnrstflo),   &  
    6981                        ( tpjfl  (jfl), jfl=1, jpnrstflo),   & 
    7082                        ( tpkfl  (jfl), jfl=1, jpnrstflo),   & 
     
    7486 
    7587         ! if we want a  surface drift  ( like PROVOR floats ) 
    76          IF( ln_argo ) THEN 
    77             DO jfl = 1, jpnrstflo 
    78                nisobfl(jfl) = 0 
    79             END DO 
    80          ENDIF 
    81  
    82          IF(lwp) WRITE(numout,*)' flo_dom: END of florstlec' 
     88         IF( ln_argo ) nisobfl(1:jpnrstflo) = 0 
    8389          
    8490         ! It is possible to add new floats.           
    85          IF(lwp) WRITE(numout,*)' flo_dom:jpnfl jpnrstflo ',jpnfl,jpnrstflo 
    86          IF( jpnfl > jpnrstflo ) THEN 
    87             ! open the init file  
    88             CALL ctl_opn( inum, 'init_float', 'OLD', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp ) 
    89             DO jfl = jpnrstflo+1, jpnfl 
    90                READ(inum,*) flxx(jfl),flyy(jfl),flzz(jfl), nisobfl(jfl),ngrpfl(jfl),jfl1 
    91             END DO 
    92             CLOSE(inum) 
    93             IF(lwp) WRITE(numout,*)' flodom: END reading init_float file' 
     91         !--------------------------------- 
     92         IF( jpnfl > jpnrstflo )THEN 
     93 
     94            IF(lwp) WRITE(numout,*) '        add new floats' 
     95 
     96            IF( ln_ariane )THEN  !Add new floats with ariane convention 
     97                CALL flo_add_new_ariane_floats(jpnrstflo+1,jpnfl)  
     98            ELSE                 !Add new floats with long/lat convention 
     99                CALL flo_add_new_floats(jpnrstflo+1,jpnfl) 
     100            ENDIF 
     101         ENDIF 
     102 
     103      !--------------------------------------! 
     104      ! FLOAT INITILISATION: NO RESTART FILE ! 
     105      !--------------------------------------! 
     106      ELSE    !ln_rstflo 
     107 
     108         IF( ln_ariane )THEN       !Add new floats with ariane convention 
     109            CALL flo_add_new_ariane_floats(1,jpnfl) 
     110         ELSE                      !Add new floats with long/lat convention 
     111            CALL flo_add_new_floats(1,jpnfl) 
     112         ENDIF 
     113 
     114      ENDIF 
    94115             
    95             ! Test to find the grid point coordonate with the geographical position             
    96             DO jfl = jpnrstflo+1, jpnfl 
    97                ihtest(jfl) = 0 
    98                ivtest(jfl) = 0 
    99                ikmfl(jfl) = 0 
     116   END SUBROUTINE flo_dom 
     117 
     118   SUBROUTINE flo_add_new_floats(kfl_start, kfl_end) 
     119      !! ------------------------------------------------------------- 
     120      !!                 ***  SUBROUTINE add_new_arianefloats  *** 
     121      !!           
     122      !! ** Purpose :    
     123      !! 
     124      !!       First initialisation of floats 
     125      !!       the initials positions of floats are written in a file 
     126      !!       with a variable to know if it is a isobar float a number  
     127      !!       to identified who want the trajectories of this float and  
     128      !!       an index for the number of the float          
     129      !!       open the init file  
     130      !!                
     131      !! ** Method  :  
     132      !!---------------------------------------------------------------------- 
     133      INTEGER, INTENT(in) :: kfl_start, kfl_end 
     134      !! 
     135      INTEGER           :: inum ! file unit 
     136      INTEGER           :: jfl,ji, jj, jk ! dummy loop indices 
     137      INTEGER           :: itrash         ! trash var for reading 
     138      INTEGER           :: ifl            ! number of floats to read 
     139      REAL(wp)          :: zdxab, zdyad 
     140      LOGICAL           :: llinmesh 
     141      CHARACTER(len=80) :: cltmp 
     142      !!--------------------------------------------------------------------- 
     143      ifl = kfl_end-kfl_start+1 
     144 
     145      ! we get the init values  
     146      !----------------------- 
     147      CALL ctl_opn( inum , clname1, 'OLD', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp ) 
     148      DO jfl = kfl_start,kfl_end 
     149         READ(inum,*) flxx(jfl),flyy(jfl),flzz(jfl), nisobfl(jfl),ngrpfl(jfl),itrash 
     150         if(lwp)write(numout,*)'read:',jfl,flxx(jfl),flyy(jfl),flzz(jfl), nisobfl(jfl),ngrpfl(jfl),itrash ; call flush(numout) 
     151      END DO 
     152      CLOSE(inum) 
     153             
     154      ! Test to find the grid point coordonate with the geographical position             
     155      !---------------------------------------------------------------------- 
     156      DO jfl = kfl_start,kfl_end 
     157         ihtest(jfl) = 0 
     158         ivtest(jfl) = 0 
     159         ikmfl(jfl) = 0 
    100160# if   defined key_mpp_mpi 
    101                DO ji = MAX(nldi,2), nlei 
    102                   DO jj = MAX(nldj,2), nlej   ! NO vector opt. 
    103 # else 
    104                DO ji = 2, jpi 
    105                   DO jj = 2, jpj   ! NO vector opt. 
     161         DO ji = MAX(nldi,2), nlei 
     162            DO jj = MAX(nldj,2), nlej   ! NO vector opt. 
     163# else          
     164         DO ji = 2, jpi 
     165            DO jj = 2, jpj   ! NO vector opt. 
    106166# endif                      
    107                      ! For each float we find the indexes of the mesh                       
    108                      CALL findmesh(glamf(ji-1,jj-1),gphif(ji-1,jj-1),   & 
    109                                    glamf(ji-1,jj  ),gphif(ji-1,jj  ),   & 
    110                                    glamf(ji  ,jj  ),gphif(ji  ,jj  ),   & 
    111                                    glamf(ji  ,jj-1),gphif(ji  ,jj-1),   & 
    112                                    flxx(jfl)       ,flyy(jfl)       ,   & 
    113                                    glamt(ji  ,jj  ),gphit(ji  ,jj  ), llinmesh) 
    114                      IF(llinmesh) THEN 
    115                         iimfl(jfl) = ji 
    116                         ijmfl(jfl) = jj 
    117                         ihtest(jfl) = ihtest(jfl)+1 
    118                         DO jk = 1, jpk-1 
    119                            IF( (fsdepw(ji,jj,jk) <= flzz(jfl)) .AND. (fsdepw(ji,jj,jk+1) > flzz(jfl)) ) THEN 
    120                               ikmfl(jfl) = jk 
    121                               ivtest(jfl) = ivtest(jfl) + 1 
    122                            ENDIF 
    123                         END DO 
     167               ! For each float we find the indexes of the mesh                       
     168               CALL flo_findmesh(glamf(ji-1,jj-1),gphif(ji-1,jj-1),   & 
     169                                 glamf(ji-1,jj  ),gphif(ji-1,jj  ),   & 
     170                                 glamf(ji  ,jj  ),gphif(ji  ,jj  ),   & 
     171                                 glamf(ji  ,jj-1),gphif(ji  ,jj-1),   & 
     172                                 flxx(jfl)       ,flyy(jfl)       ,   & 
     173                                 glamt(ji  ,jj  ),gphit(ji  ,jj  ), llinmesh) 
     174               IF( llinmesh )THEN 
     175                  iimfl(jfl) = ji 
     176                  ijmfl(jfl) = jj 
     177                  ihtest(jfl) = ihtest(jfl)+1 
     178                  DO jk = 1, jpk-1 
     179                     IF( (fsdepw(ji,jj,jk) <= flzz(jfl)) .AND. (fsdepw(ji,jj,jk+1) > flzz(jfl)) ) THEN 
     180                        ikmfl(jfl) = jk 
     181                        ivtest(jfl) = ivtest(jfl) + 1 
    124182                     ENDIF 
    125183                  END DO 
    126                END DO 
    127                IF(lwp) WRITE(numout,*)'   flo_dom: END findmesh' 
    128                 
    129                ! If the float is in a mesh computed by an other processor we put iimfl=ijmfl=-1                
    130                IF( ihtest(jfl) ==  0 ) THEN 
    131                   iimfl(jfl) = -1 
    132                   ijmfl(jfl) = -1 
    133184               ENDIF 
    134185            END DO 
     186         END DO 
     187 
     188         ! If the float is in a mesh computed by an other processor we put iimfl=ijmfl=-1                
     189         IF( ihtest(jfl) ==  0 ) THEN 
     190            iimfl(jfl) = -1 
     191            ijmfl(jfl) = -1 
     192         ENDIF 
     193      END DO 
     194 
     195      !Test if each float is in one and only one proc 
     196      !---------------------------------------------- 
     197      IF( lk_mpp )   THEN  
     198         CALL mpp_sum(ihtest,jpnfl) 
     199         CALL mpp_sum(ivtest,jpnfl) 
     200      ENDIF 
     201      DO jfl = kfl_start,kfl_end 
     202 
     203         IF( (ihtest(jfl) > 1 ) .OR. ( ivtest(jfl) > 1) ) THEN 
     204             WRITE(cltmp,'(A10,i4.4,A20)' )'THE FLOAT',jfl,' IS NOT IN ONLY ONE MESH' 
     205             CALL ctl_stop('STOP',TRIM(cltmp) ) 
     206         ENDIF 
     207         IF( (ihtest(jfl) == 0) ) THEN 
     208             WRITE(cltmp,'(A10,i4.4,A20)' )'THE FLOAT',jfl,' IS IN NO MESH' 
     209             CALL ctl_stop('STOP',TRIM(cltmp) ) 
     210         ENDIF 
     211      END DO 
     212 
     213      ! We compute the distance between the float and the face of the mesh             
     214      !------------------------------------------------------------------- 
     215      DO jfl = kfl_start,kfl_end 
     216 
     217         ! Made only if the float is in the domain of the processor               
     218         IF( (iimfl(jfl) >= 0) .AND. (ijmfl(jfl) >= 0) ) THEN 
     219 
     220            ! TEST TO KNOW IF THE FLOAT IS NOT INITIALISED IN THE COAST 
     221            idomfl(jfl) = 0 
     222            IF( tmask(iimfl(jfl),ijmfl(jfl),ikmfl(jfl)) == 0. ) idomfl(jfl) = 1 
     223 
     224            ! Computation of the distance between the float and the faces of the mesh 
     225            !            zdxab 
     226            !             . 
     227            !        B----.---------C 
     228            !        |    .         | 
     229            !        |<------>flo   | 
     230            !        |        ^     | 
     231            !        |        |.....|....zdyad 
     232            !        |        |     | 
     233            !        A--------|-----D 
     234            ! 
     235            zdxab = flo_dstnce( flxx(jfl), flyy(jfl), glamf(iimfl(jfl)-1,ijmfl(jfl)-1), flyy(jfl) ) 
     236            zdyad = flo_dstnce( flxx(jfl), flyy(jfl), flxx(jfl), gphif(iimfl(jfl)-1,ijmfl(jfl)-1) ) 
     237 
     238            ! Translation of this distances (in meter) in indexes 
     239            zgifl(jfl)= (iimfl(jfl)-0.5) + zdxab/e1u(iimfl(jfl)-1,ijmfl(jfl)) + (mig(1)-jpizoom) 
     240            zgjfl(jfl)= (ijmfl(jfl)-0.5) + zdyad/e2v(iimfl(jfl),ijmfl(jfl)-1) + (mjg(1)-jpjzoom) 
     241            zgkfl(jfl) = (( fsdepw(iimfl(jfl),ijmfl(jfl),ikmfl(jfl)+1) - flzz(jfl) )* ikmfl(jfl))   & 
     242               &                 / (  fsdepw(iimfl(jfl),ijmfl(jfl),ikmfl(jfl)+1)                              & 
     243               &                    - fsdepw(iimfl(jfl),ijmfl(jfl),ikmfl(jfl) ) )                             & 
     244               &                 + (( flzz(jfl)-fsdepw(iimfl(jfl),ijmfl(jfl),ikmfl(jfl)) ) *(ikmfl(jfl)+1))   & 
     245               &                 / (  fsdepw(iimfl(jfl),ijmfl(jfl),ikmfl(jfl)+1)                              & 
     246               &                    - fsdepw(iimfl(jfl),ijmfl(jfl),ikmfl(jfl)) ) 
     247         ELSE 
     248            zgifl(jfl) = 0.e0 
     249            zgjfl(jfl) = 0.e0 
     250            zgkfl(jfl) = 0.e0 
     251         ENDIF 
     252 
     253      END DO 
     254                   
     255      ! The sum of all the arrays zgifl, zgjfl, zgkfl give 3 arrays with the positions of all the floats. 
     256      IF( lk_mpp )   THEN  
     257         CALL mpp_sum( zgjfl, ifl )   ! sums over the global domain 
     258         CALL mpp_sum( zgkfl, ifl ) 
     259      ENDIF 
    135260             
    136             ! A zero in the sum of the arrays "ihtest" and "ivtest"              
    137 # if   defined key_mpp_mpi 
    138             CALL mpp_sum(ihtest,jpnfl) 
    139             CALL mpp_sum(ivtest,jpnfl) 
    140 # endif  
    141             DO jfl = jpnrstflo+1, jpnfl 
    142                IF( (ihtest(jfl) > 1 ) .OR. ( ivtest(jfl) > 1) ) THEN 
    143                   IF(lwp) WRITE(numout,*) 'THE FLOAT',jfl,' IS NOT IN ONLY ONE MESH' 
    144                   STOP 
    145                ENDIF 
    146                IF( (ihtest(jfl) == 0) ) THEN 
    147                   IF(lwp) WRITE(numout,*)'THE FLOAT',jfl,' IS IN NO MESH' 
    148                   STOP 
    149                ENDIF 
    150             END DO 
    151              
    152             ! We compute the distance between the float and the face of the mesh             
    153             DO jfl = jpnrstflo+1, jpnfl                
    154                ! Made only if the float is in the domain of the processor               
    155                IF( (iimfl(jfl) >= 0) .AND. (ijmfl(jfl) >= 0) ) THEN 
    156                    
    157                   ! TEST TO KNOW IF THE FLOAT IS NOT INITIALISED IN THE COAST 
    158                    
    159                   idomfl(jfl) = 0 
    160                   IF( tmask(iimfl(jfl),ijmfl(jfl),ikmfl(jfl)) == 0. ) idomfl(jfl) = 1 
    161                                             
    162                   ! Computation of the distance between the float and the faces of the mesh 
    163                   !            zdxab 
    164                   !             . 
    165                   !        B----.---------C 
    166                   !        |    .         | 
    167                   !        |<------>flo   | 
    168                   !        |        ^     | 
    169                   !        |        |.....|....zdyad 
    170                   !        |        |     | 
    171                   !        A--------|-----D 
    172                   ! 
    173               
    174                   zdxab = dstnce( flxx(jfl), flyy(jfl), glamf(iimfl(jfl)-1,ijmfl(jfl)-1), flyy(jfl) ) 
    175                   zdyad = dstnce( flxx(jfl), flyy(jfl), flxx(jfl), gphif(iimfl(jfl)-1,ijmfl(jfl)-1) ) 
    176                    
    177                   ! Translation of this distances (in meter) in indexes 
    178                    
    179                   zgifl(jfl-jpnrstflo)= (iimfl(jfl)-0.5) + zdxab/e1u(iimfl(jfl)-1,ijmfl(jfl)) + (mig(1)-jpizoom) 
    180                   zgjfl(jfl-jpnrstflo)= (ijmfl(jfl)-0.5) + zdyad/e2v(iimfl(jfl),ijmfl(jfl)-1) + (mjg(1)-jpjzoom) 
    181                   zgkfl(jfl-jpnrstflo) = (( fsdepw(iimfl(jfl),ijmfl(jfl),ikmfl(jfl)+1) - flzz(jfl) )* ikmfl(jfl))   & 
    182                      &                 / (  fsdepw(iimfl(jfl),ijmfl(jfl),ikmfl(jfl)+1)                              & 
    183                      &                    - fsdepw(iimfl(jfl),ijmfl(jfl),ikmfl(jfl) ) )                             & 
    184                      &                 + (( flzz(jfl)-fsdepw(iimfl(jfl),ijmfl(jfl),ikmfl(jfl)) ) *(ikmfl(jfl)+1))   & 
    185                      &                 / (  fsdepw(iimfl(jfl),ijmfl(jfl),ikmfl(jfl)+1)                              & 
    186                      &                    - fsdepw(iimfl(jfl),ijmfl(jfl),ikmfl(jfl)) ) 
    187                ELSE 
    188                   zgifl(jfl-jpnrstflo) = 0.e0 
    189                   zgjfl(jfl-jpnrstflo) = 0.e0 
    190                   zgkfl(jfl-jpnrstflo) = 0.e0 
    191                ENDIF 
    192             END DO 
    193              
    194             ! The sum of all the arrays zgifl, zgjfl, zgkfl give 3 arrays with the positions of all the floats. 
    195             IF( lk_mpp )   THEN 
    196                CALL mpp_sum( zgjfl, jpnnewflo )   ! sums over the global domain 
    197                CALL mpp_sum( zgkfl, jpnnewflo ) 
    198                IF(lwp) WRITE(numout,*) (zgifl(jfl),jfl=1,jpnnewflo) 
    199                IF(lwp) WRITE(numout,*) (zgjfl(jfl),jfl=1,jpnnewflo) 
    200                IF(lwp) WRITE(numout,*) (zgkfl(jfl),jfl=1,jpnnewflo)  
    201             ENDIF 
    202             
    203             DO jfl = jpnrstflo+1, jpnfl 
    204                tpifl(jfl) = zgifl(jfl-jpnrstflo) 
    205                tpjfl(jfl) = zgjfl(jfl-jpnrstflo) 
    206                tpkfl(jfl) = zgkfl(jfl-jpnrstflo) 
    207             END DO 
    208          ENDIF 
    209       ELSE 
    210          IF(lwp) WRITE(numout,*) '                     init_float read ' 
    211           
    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 
    231 # 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. 
    237 # endif                   
    238                   ! for each float we find the indexes of the mesh  
    239                    
    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 
    257                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 
    266           
    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 
    279          
    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 
    320           
    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 ) 
    326       ENDIF 
    327              
    328       ! Print the initial positions of the floats 
     261      DO jfl = kfl_start,kfl_end 
     262         tpifl(jfl) = zgifl(jfl) 
     263         tpjfl(jfl) = zgjfl(jfl) 
     264         tpkfl(jfl) = zgkfl(jfl) 
     265      END DO 
     266 
     267      ! WARNING : initial position not in the sea          
    329268      IF( .NOT. ln_rstflo ) THEN  
    330          ! WARNING : initial position not in the sea          
    331          DO jfl = 1, jpnfl 
     269         DO jfl =  kfl_start,kfl_end 
    332270            IF( idomfl(jfl) == 1 ) THEN 
    333271               IF(lwp) WRITE(numout,*)'*****************************' 
     
    341279      ENDIF 
    342280 
    343    END SUBROUTINE flo_dom 
    344  
    345  
    346    SUBROUTINE findmesh( pax, pay, pbx, pby,   & 
    347                         pcx, pcy, pdx, pdy,   & 
    348                         px  ,py  ,ptx, pty, ldinmesh ) 
     281   END SUBROUTINE flo_add_new_floats 
     282 
     283   SUBROUTINE flo_add_new_ariane_floats(kfl_start, kfl_end) 
     284      !! ------------------------------------------------------------- 
     285      !!                 ***  SUBROUTINE add_new_arianefloats  *** 
     286      !!           
     287      !! ** Purpose :    
     288      !!       First initialisation of floats with ariane convention 
     289      !!        
     290      !!       The indexes are read directly from file (warning ariane 
     291      !!       convention, are refered to  
     292      !!       U,V,W grids - and not T-)  
     293      !!       The isobar advection is managed with the sign of tpkfl ( >0 -> 3D 
     294      !!       advection, <0 -> 2D)  
     295      !!       Some variables are not read, as - gl         : time index; 4th 
     296      !!       column         
     297      !!                                       - transport  : transport ; 5th 
     298      !!                                       column 
     299      !!       and paste in the jtrash var 
     300      !!       At the end, ones need to replace the indexes on T grid 
     301      !!       RMQ : there is no float groups identification ! 
     302      !! 
     303      !!                
     304      !! ** Method  :  
     305      !!---------------------------------------------------------------------- 
     306      INTEGER, INTENT(in) :: kfl_start, kfl_end 
     307      !! 
     308      INTEGER  :: inum         ! file unit 
     309      INTEGER  :: ierr, ifl 
     310      INTEGER  :: jfl, jfl1    ! dummy loop indices 
     311      INTEGER  :: itrash       ! trash var for reading   
     312      CHARACTER(len=80) :: cltmp 
     313 
     314      !!---------------------------------------------------------------------- 
     315      nisobfl(kfl_start:kfl_end) = 1 ! we assume that by default we want 3D advection 
     316 
     317      ifl = kfl_end - kfl_start + 1  ! number of floats to read   
     318 
     319      ! we check that the number of floats in the init_file are consistant with the namelist 
     320      IF( lwp ) THEN 
     321 
     322         jfl1=0 
     323         ierr=0 
     324         CALL ctl_opn( inum, clname2, 'OLD', 'FORMATTED', 'SEQUENTIAL',  1, numout, .TRUE., 1 ) 
     325         DO WHILE (ierr .EQ. 0) 
     326            jfl1=jfl1+1 
     327            READ(inum,*, iostat=ierr) 
     328         END DO 
     329         CLOSE(inum) 
     330         IF( (jfl1-1) .NE. ifl )THEN  
     331            WRITE(cltmp,'(A25,A20,A3,i4.4,A10,i4.4)')"the number of floats in ",TRIM(clname2), & 
     332                                                     " = ",jfl1," is not equal to jfl= ",ifl 
     333            CALL ctl_stop('STOP',TRIM(cltmp) ) 
     334         ENDIF 
     335 
     336      ENDIF 
     337             
     338      ! we get the init values  
     339      CALL ctl_opn( inum, clname2, 'OLD', 'FORMATTED', 'SEQUENTIAL', 1, numout, .TRUE., 1 ) 
     340      DO jfl = kfl_start, kfl_end 
     341          READ(inum,*) tpifl(jfl),tpjfl(jfl),tpkfl(jfl),itrash, itrash 
     342               
     343          IF ( tpkfl(jfl) .LT. 0. ) nisobfl(jfl) = 0 !set the 2D advection according to init_float 
     344          ngrpfl(jfl)=jfl 
     345      END DO 
     346 
     347      ! conversion from ariane index to T grid index 
     348      tpkfl(kfl_start:kfl_end) = abs(tpkfl)-0.5 ! reversed vertical axis 
     349      tpifl(kfl_start:kfl_end) = tpifl+0.5 
     350      tpjfl(kfl_start:kfl_end) = tpjfl+0.5 
     351 
     352 
     353   END SUBROUTINE flo_add_new_ariane_floats 
     354 
     355 
     356   SUBROUTINE flo_findmesh( pax, pay, pbx, pby,   & 
     357                            pcx, pcy, pdx, pdy,   & 
     358                            px  ,py  ,ptx, pty, ldinmesh ) 
    349359      !! ------------------------------------------------------------- 
    350360      !!                ***  ROUTINE findmesh  *** 
     
    402412      ENDIF 
    403413      ! 
    404    END SUBROUTINE findmesh 
    405  
    406  
    407    FUNCTION dstnce( pla1, phi1, pla2, phi2 ) 
     414   END SUBROUTINE flo_findmesh 
     415 
     416 
     417   FUNCTION flo_dstnce( pla1, phi1, pla2, phi2 ) 
    408418      !! ------------------------------------------------------------- 
    409419      !!                 ***  Function dstnce  *** 
     
    415425      REAL(wp), INTENT(in) ::   pla1, phi1, pla2, phi2   ! ??? 
    416426      !! 
    417       REAL(wp) ::   dly1, dly2, dlx1, dlx2, dlx, dls, dld, dpi 
    418       REAL(wp) ::   dstnce 
     427      REAL(wp) :: dly1, dly2, dlx1, dlx2, dlx, dls, dld, dpi 
     428      REAL(wp) :: flo_dstnce 
    419429      !!--------------------------------------------------------------------- 
    420430      ! 
    421       dpi  = 2.* ASIN(1.) 
    422       dls  = dpi / 180. 
     431      dpi  = 2._wp * ASIN(1._wp) 
     432      dls  = dpi / 180._wp 
    423433      dly1 = phi1 * dls 
    424434      dly2 = phi2 * dls 
     
    428438      dlx = SIN(dly1) * SIN(dly2) + COS(dly1) * COS(dly2) * COS(dlx2-dlx1) 
    429439      ! 
    430       IF( ABS(dlx) > 1.0 ) dlx = 1.0 
    431       ! 
    432       dld = ATAN(DSQRT( ( 1-dlx )/( 1+dlx ) )) * 222.24 / dls 
    433       dstnce = dld * 1000. 
    434       ! 
    435    END FUNCTION dstnce 
    436  
    437 #  else 
     440      IF( ABS(dlx) > 1.0_wp ) dlx = 1.0_wp 
     441      ! 
     442      dld = ATAN(DSQRT( 1._wp * ( 1._wp-dlx )/( 1._wp+dlx ) )) * 222.24_wp / dls 
     443      flo_dstnce = dld * 1000._wp 
     444      ! 
     445   END FUNCTION flo_dstnce 
     446 
     447   INTEGER FUNCTION flo_dom_alloc() 
     448      !!---------------------------------------------------------------------- 
     449      !!                 ***  FUNCTION flo_dom_alloc  *** 
     450      !!---------------------------------------------------------------------- 
     451 
     452      ALLOCATE( iimfl(jpnfl) , ijmfl(jpnfl) , ikmfl(jpnfl) ,                          &   
     453                idomfl(jpnfl), ivtest(jpnfl), ihtest(jpnfl),                 & 
     454                zgifl(jpnfl) , zgjfl(jpnfl) , zgkfl(jpnfl)   , STAT=flo_dom_alloc ) 
     455      ! 
     456      IF( lk_mpp             )   CALL mpp_sum ( flo_dom_alloc ) 
     457      IF( flo_dom_alloc /= 0 )   CALL ctl_warn('flo_dom_alloc: failed to allocate arrays') 
     458   END FUNCTION flo_dom_alloc 
     459 
     460 
     461#else 
    438462   !!---------------------------------------------------------------------- 
    439463   !!   Default option                                         Empty module 
     
    441465CONTAINS 
    442466   SUBROUTINE flo_dom                 ! Empty routine 
     467         WRITE(*,*) 'flo_dom: : You should not have seen this print! error?' 
    443468   END SUBROUTINE flo_dom 
    444469#endif 
Note: See TracChangeset for help on using the changeset viewer.