MODULE floblk !!====================================================================== !! *** MODULE floblk *** !! Ocean floats : trajectory computation !!====================================================================== #if defined key_floats || defined key_esopa !!---------------------------------------------------------------------- !! 'key_floats' float trajectories !!---------------------------------------------------------------------- !!---------------------------------------------------------------------- !! flotblk : compute float trajectories with Blanke algorithme !!---------------------------------------------------------------------- !! * Modules used USE flo_oce ! ocean drifting floats USE oce ! ocean dynamics and tracers USE dom_oce ! ocean space and time domain USE phycst ! physical constants USE in_out_manager ! I/O manager USE lib_mpp ! distribued memory computing library IMPLICIT NONE PRIVATE !! * Accessibility PUBLIC flo_blk ! routine called by floats.F90 !! * Substitutions # include "domzgr_substitute.h90" !!---------------------------------------------------------------------- !! OPA 9.0 , LOCEAN-IPSL (2005) !! $Header$ !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt !!---------------------------------------------------------------------- CONTAINS SUBROUTINE flo_blk( kt ) !!--------------------------------------------------------------------- !! *** ROUTINE flo_blk *** !! !! ** Purpose : Compute the geographical position,latitude, longitude !! and depth of each float at each time step. !! !! ** Method : The position of a float is computed with Bruno Blanke !! algorithm. We need to know the velocity field, the old positions !! of the floats and the grid defined on the domain. !! !!---------------------------------------------------------------------- !! * arguments INTEGER, INTENT( in ) :: kt ! ocean time step !! * Local declarations INTEGER :: jfl ! dummy loop arguments INTEGER :: ind, ifin, iloop INTEGER , DIMENSION ( jpnfl ) :: & iil, ijl, ikl, & ! index of nearest mesh iiloc , ijloc, & iiinfl, ijinfl, ikinfl, & ! index of input mesh of the float. iioutfl, ijoutfl, ikoutfl ! index of output mesh of the float. REAL(wp) , DIMENSION ( jpnfl ) :: & zgifl, zgjfl, zgkfl, & ! position of floats, index on ! velocity mesh. ztxfl, ztyfl, ztzfl, & ! time for a float to quit the mesh ! across one of the face x,y and z zttfl, & ! time for a float to quit the mesh zagefl, & ! time during which, trajectorie of ! the float has been computed zagenewfl, & ! new age of float after calculation ! of new position zufl, zvfl, zwfl, & ! interpolated vel. at float position zudfl, zvdfl, zwdfl, & ! velocity diff input/output of mesh zgidfl, zgjdfl, zgkdfl ! direction index of float REAL(wp) :: & zuinfl,zvinfl,zwinfl, & ! transport across the input face zuoutfl,zvoutfl,zwoutfl, & ! transport across the ouput face zvol, & ! volume of the mesh zsurfz, & ! surface of the face of the mesh zind REAL(wp), DIMENSION ( 2 ) :: & zsurfx, zsurfy ! surface of the face of the mesh !!--------------------------------------------------------------------- IF( kt == nit000 ) THEN IF(lwp) WRITE(numout,*) IF(lwp) WRITE(numout,*) 'flo_blk : compute Blanke trajectories for floats ' IF(lwp) WRITE(numout,*) '~~~~~~~ ' ENDIF ! Initialisation of parameters DO jfl = 1, jpnfl ! ages of floats are put at zero zagefl(jfl) = 0. ! index on the velocity grid ! We considere k coordinate negative, with this transformation ! the computation in the 3 direction is the same. zgifl(jfl) = tpifl(jfl) - 0.5 zgjfl(jfl) = tpjfl(jfl) - 0.5 zgkfl(jfl) = MIN(-1.,-(tpkfl(jfl))) ! surface drift every 10 days IF( ln_argo ) THEN IF( MOD(kt,150) >= 146 .OR. MOD(kt,150) == 0 ) zgkfl(jfl) = -1. ENDIF ! index of T mesh iil(jfl) = 1 + INT(zgifl(jfl)) ijl(jfl) = 1 + INT(zgjfl(jfl)) ikl(jfl) = INT(zgkfl(jfl)) END DO iloop = 0 222 DO jfl = 1, jpnfl # if defined key_mpp_mpi || defined key_mpp_shmem IF( (iil(jfl) >= (mig(nldi)-jpizoom+1)) .AND. (iil(jfl) <= (mig(nlei)-jpizoom+1)) .AND. & (ijl(jfl) >= (mjg(nldj)-jpjzoom+1)) .AND. (ijl(jfl) <= (mjg(nlej)-jpjzoom+1)) ) THEN iiloc(jfl) = iil(jfl) - (mig(1)-jpizoom+1) + 1 ijloc(jfl) = ijl(jfl) - (mjg(1)-jpjzoom+1) + 1 # else iiloc(jfl) = iil(jfl) ijloc(jfl) = ijl(jfl) # endif ! compute the transport across the mesh where the float is. !!bug (gm) change e3t into fse3. but never checked zsurfx(1) = e2u(iiloc(jfl)-1,ijloc(jfl) ) * fse3u(iiloc(jfl)-1,ijloc(jfl) ,-ikl(jfl)) zsurfx(2) = e2u(iiloc(jfl) ,ijloc(jfl) ) * fse3u(iiloc(jfl) ,ijloc(jfl) ,-ikl(jfl)) zsurfy(1) = e1v(iiloc(jfl) ,ijloc(jfl)-1) * fse3v(iiloc(jfl) ,ijloc(jfl)-1,-ikl(jfl)) zsurfy(2) = e1v(iiloc(jfl) ,ijloc(jfl) ) * fse3v(iiloc(jfl) ,ijloc(jfl) ,-ikl(jfl)) ! for a isobar float zsurfz is put to zero. The vertical velocity will be zero too. zsurfz = e1t(iiloc(jfl),ijloc(jfl)) * e2t(iiloc(jfl),ijloc(jfl)) zvol = zsurfz * fse3t(iiloc(jfl),ijloc(jfl),-ikl(jfl)) ! zuinfl =( ub(iiloc(jfl)-1,ijloc(jfl),-ikl(jfl)) + un(iiloc(jfl)-1,ijloc(jfl),-ikl(jfl)) )/2.*zsurfx(1) zuoutfl=( ub(iiloc(jfl) ,ijloc(jfl),-ikl(jfl)) + un(iiloc(jfl) ,ijloc(jfl),-ikl(jfl)) )/2.*zsurfx(2) zvinfl =( vb(iiloc(jfl),ijloc(jfl)-1,-ikl(jfl)) + vn(iiloc(jfl),ijloc(jfl)-1,-ikl(jfl)) )/2.*zsurfy(1) zvoutfl=( vb(iiloc(jfl),ijloc(jfl) ,-ikl(jfl)) + vn(iiloc(jfl),ijloc(jfl) ,-ikl(jfl)) )/2.*zsurfy(2) zwinfl =-(wb(iiloc(jfl),ijloc(jfl),-(ikl(jfl)-1)) & & + wn(iiloc(jfl),ijloc(jfl),-(ikl(jfl)-1)) )/2. * zsurfz*nisobfl(jfl) zwoutfl=-(wb(iiloc(jfl),ijloc(jfl),- ikl(jfl) ) & & + wn(iiloc(jfl),ijloc(jfl),- ikl(jfl) ) )/2. * zsurfz*nisobfl(jfl) ! interpolation of velocity field on the float initial position zufl(jfl)= zuinfl + ( zgifl(jfl) - float(iil(jfl)-1) ) * ( zuoutfl - zuinfl) zvfl(jfl)= zvinfl + ( zgjfl(jfl) - float(ijl(jfl)-1) ) * ( zvoutfl - zvinfl) zwfl(jfl)= zwinfl + ( zgkfl(jfl) - float(ikl(jfl)-1) ) * ( zwoutfl - zwinfl) ! faces of input and output ! u-direction IF( zufl(jfl) < 0. ) THEN iioutfl(jfl) = iil(jfl) - 1. iiinfl (jfl) = iil(jfl) zind = zuinfl zuinfl = zuoutfl zuoutfl= zind ELSE iioutfl(jfl) = iil(jfl) iiinfl (jfl) = iil(jfl) - 1 ENDIF ! v-direction IF( zvfl(jfl) < 0. ) THEN ijoutfl(jfl) = ijl(jfl) - 1. ijinfl (jfl) = ijl(jfl) zind = zvinfl zvinfl = zvoutfl zvoutfl = zind ELSE ijoutfl(jfl) = ijl(jfl) ijinfl (jfl) = ijl(jfl) - 1. ENDIF ! w-direction IF( zwfl(jfl) < 0. ) THEN ikoutfl(jfl) = ikl(jfl) - 1. ikinfl (jfl) = ikl(jfl) zind = zwinfl zwinfl = zwoutfl zwoutfl = zind ELSE ikoutfl(jfl) = ikl(jfl) ikinfl (jfl) = ikl(jfl) - 1. ENDIF ! compute the time to go out the mesh across a face ! u-direction zudfl (jfl) = zuoutfl - zuinfl zgidfl(jfl) = float(iioutfl(jfl) - iiinfl(jfl)) IF( zufl(jfl)*zuoutfl <= 0. ) THEN ztxfl(jfl) = 1.E99 ELSE IF( ABS(zudfl(jfl)) >= 1.E-5 ) THEN ztxfl(jfl)= zgidfl(jfl)/zudfl(jfl) * LOG(zuoutfl/zufl (jfl)) ELSE ztxfl(jfl)=(float(iioutfl(jfl))-zgifl(jfl))/zufl(jfl) ENDIF IF( (ABS(zgifl(jfl)-float(iiinfl (jfl))) <= 1.E-7) .OR. & (ABS(zgifl(jfl)-float(iioutfl(jfl))) <= 1.E-7) ) THEN ztxfl(jfl)=(zgidfl(jfl))/zufl(jfl) ENDIF ENDIF ! v-direction zvdfl (jfl) = zvoutfl - zvinfl zgjdfl(jfl) = float(ijoutfl(jfl)-ijinfl(jfl)) IF( zvfl(jfl)*zvoutfl <= 0. ) THEN ztyfl(jfl) = 1.E99 ELSE IF( ABS(zvdfl(jfl)) >= 1.E-5 ) THEN ztyfl(jfl) = zgjdfl(jfl)/zvdfl(jfl) * LOG(zvoutfl/zvfl (jfl)) ELSE ztyfl(jfl) = (float(ijoutfl(jfl)) - zgjfl(jfl))/zvfl(jfl) ENDIF IF( (ABS(zgjfl(jfl)-float(ijinfl (jfl))) <= 1.E-7) .OR. & (ABS(zgjfl(jfl)-float(ijoutfl(jfl))) <= 1.E-7) ) THEN ztyfl(jfl) = (zgjdfl(jfl)) / zvfl(jfl) ENDIF ENDIF ! w-direction IF( nisobfl(jfl) == 1. ) THEN zwdfl (jfl) = zwoutfl - zwinfl zgkdfl(jfl) = float(ikoutfl(jfl) - ikinfl(jfl)) IF( zwfl(jfl)*zwoutfl <= 0. ) THEN ztzfl(jfl) = 1.E99 ELSE IF( ABS(zwdfl(jfl)) >= 1.E-5 ) THEN ztzfl(jfl) = zgkdfl(jfl)/zwdfl(jfl) * LOG(zwoutfl/zwfl (jfl)) ELSE ztzfl(jfl) = (float(ikoutfl(jfl)) - zgkfl(jfl))/zwfl(jfl) ENDIF IF( (ABS(zgkfl(jfl)-float(ikinfl (jfl))) <= 1.E-7) .OR. & (ABS(zgkfl(jfl)-float(ikoutfl(jfl))) <= 1.E-7) ) THEN ztzfl(jfl) = (zgkdfl(jfl)) / zwfl(jfl) ENDIF ENDIF ENDIF ! the time to go leave the mesh is the smallest time IF( nisobfl(jfl) == 1. ) THEN zttfl(jfl) = MIN(ztxfl(jfl),ztyfl(jfl),ztzfl(jfl)) ELSE zttfl(jfl) = MIN(ztxfl(jfl),ztyfl(jfl)) ENDIF ! new age of the FLOAT zagenewfl(jfl) = zagefl(jfl) + zttfl(jfl)*zvol ! test to know if the "age" of the float is not bigger than the ! time step IF( zagenewfl(jfl) > rdt ) THEN zttfl(jfl) = (rdt-zagefl(jfl)) / zvol zagenewfl(jfl) = rdt ENDIF ! In the "minimal" direction we compute the index of new mesh ! on i-direction IF( ztxfl(jfl) <= zttfl(jfl) ) THEN zgifl(jfl) = float(iioutfl(jfl)) ind = iioutfl(jfl) IF( iioutfl(jfl) >= iiinfl(jfl) ) THEN iioutfl(jfl) = iioutfl(jfl) + 1 ELSE iioutfl(jfl) = iioutfl(jfl) - 1 ENDIF iiinfl(jfl) = ind ELSE IF( ABS(zudfl(jfl)) >= 1.E-5 ) THEN zgifl(jfl) = zgifl(jfl) + zgidfl(jfl)*zufl(jfl) & & * ( EXP( zudfl(jfl)/zgidfl(jfl)*zttfl(jfl) ) - 1. ) / zudfl(jfl) ELSE zgifl(jfl) = zgifl(jfl) + zufl(jfl) * zttfl(jfl) ENDIF ENDIF ! on j-direction IF( ztyfl(jfl) <= zttfl(jfl) ) THEN zgjfl(jfl) = float(ijoutfl(jfl)) ind = ijoutfl(jfl) IF( ijoutfl(jfl) >= ijinfl(jfl) ) THEN ijoutfl(jfl) = ijoutfl(jfl) + 1 ELSE ijoutfl(jfl) = ijoutfl(jfl) - 1 ENDIF ijinfl(jfl) = ind ELSE IF( ABS(zvdfl(jfl)) >= 1.E-5 ) THEN zgjfl(jfl) = zgjfl(jfl)+zgjdfl(jfl)*zvfl(jfl) & & * ( EXP(zvdfl(jfl)/zgjdfl(jfl)*zttfl(jfl)) - 1. ) / zvdfl(jfl) ELSE zgjfl(jfl) = zgjfl(jfl)+zvfl(jfl)*zttfl(jfl) ENDIF ENDIF ! on k-direction IF( nisobfl(jfl) == 1. ) THEN IF( ztzfl(jfl) <= zttfl(jfl) ) THEN zgkfl(jfl) = float(ikoutfl(jfl)) ind = ikoutfl(jfl) IF( ikoutfl(jfl) >= ikinfl(jfl) ) THEN ikoutfl(jfl) = ikoutfl(jfl)+1 ELSE ikoutfl(jfl) = ikoutfl(jfl)-1 ENDIF ikinfl(jfl) = ind ELSE IF( ABS(zwdfl(jfl)) >= 1.E-5 ) THEN zgkfl(jfl) = zgkfl(jfl)+zgkdfl(jfl)*zwfl(jfl) & & * ( EXP(zwdfl(jfl)/zgkdfl(jfl)*zttfl(jfl)) - 1. ) / zwdfl(jfl) ELSE zgkfl(jfl) = zgkfl(jfl)+zwfl(jfl)*zttfl(jfl) ENDIF ENDIF ENDIF ! coordinate of the new point on the temperature grid iil(jfl) = MAX(iiinfl(jfl),iioutfl(jfl)) ijl(jfl) = MAX(ijinfl(jfl),ijoutfl(jfl)) IF( nisobfl(jfl) == 1 ) ikl(jfl) = MAX(ikinfl(jfl),ikoutfl(jfl)) !!Alexcadm write(*,*)'PE ',narea, !!Alexcadm . iiinfl(jfl),iioutfl(jfl),ijinfl(jfl) !!Alexcadm . ,ijoutfl(jfl),ikinfl(jfl), !!Alexcadm . ikoutfl(jfl),ztxfl(jfl),ztyfl(jfl) !!Alexcadm . ,ztzfl(jfl),zgifl(jfl), !!Alexcadm . zgjfl(jfl) !!Alexcadm IF (jfl == 910) write(*,*)'Flotteur 910', !!Alexcadm . iiinfl(jfl),iioutfl(jfl),ijinfl(jfl) !!Alexcadm . ,ijoutfl(jfl),ikinfl(jfl), !!Alexcadm . ikoutfl(jfl),ztxfl(jfl),ztyfl(jfl) !!Alexcadm . ,ztzfl(jfl),zgifl(jfl), !!Alexcadm . zgjfl(jfl) ! reinitialisation of the age of FLOAT zagefl(jfl) = zagenewfl(jfl) # if defined key_mpp_mpi || defined key_mpp_shmem ELSE ! we put zgifl, zgjfl, zgkfl, zagefl zgifl (jfl) = 0. zgjfl (jfl) = 0. zgkfl (jfl) = 0. zagefl(jfl) = 0. iil(jfl) = 0 ijl(jfl) = 0 ENDIF # endif END DO ! synchronisation IF( lk_mpp ) CALL mpp_sum( zgifl , jpnfl ) ! sums over the global domain IF( lk_mpp ) CALL mpp_sum( zgjfl , jpnfl ) IF( lk_mpp ) CALL mpp_sum( zgkfl , jpnfl ) IF( lk_mpp ) CALL mpp_sum( zagefl, jpnfl ) IF( lk_mpp ) CALL mpp_sum( iil , jpnfl ) IF( lk_mpp ) CALL mpp_sum( ijl , jpnfl ) ! in the case of open boundaries we need to test if the floats don't ! go out of the domain. If it goes out, the float is put at the ! middle of the mesh in the domain but the trajectory isn't compute ! more time. # if defined key_obc DO jfl = 1, jpnfl IF( lp_obc_east ) THEN IF( jped <= zgjfl(jfl) .AND. zgjfl(jfl) <= jpef .AND. nieob-1 <= zgifl(jfl) ) THEN zgifl (jfl) = INT(zgifl(jfl)) + 0.5 zgjfl (jfl) = INT(zgjfl(jfl)) + 0.5 zagefl(jfl) = rdt END IF END IF IF( lp_obc_west ) THEN IF( jpwd <= zgjfl(jfl) .AND. zgjfl(jfl) <= jpwf .AND. niwob >= zgifl(jfl) ) THEN zgifl (jfl) = INT(zgifl(jfl)) + 0.5 zgjfl (jfl) = INT(zgjfl(jfl)) + 0.5 zagefl(jfl) = rdt END IF END IF IF( lp_obc_north ) THEN IF( jpnd <= zgifl(jfl) .AND. zgifl(jfl) <= jpnf .AND. njnob-1 >= zgjfl(jfl) ) THEN zgifl (jfl) = INT(zgifl(jfl)) + 0.5 zgjfl (jfl) = INT(zgjfl(jfl)) + 0.5 zagefl(jfl) = rdt END IF END IF IF( lp_obc_south ) THEN IF( jpsd <= zgifl(jfl) .AND. zgifl(jfl) <= jpsf .AND. njsob >= zgjfl(jfl) ) THEN zgifl (jfl) = INT(zgifl(jfl)) + 0.5 zgjfl (jfl) = INT(zgjfl(jfl)) + 0.5 zagefl(jfl) = rdt END IF END IF END DO #endif ! Test to know if a float hasn't integrated enought time IF( ln_argo ) THEN ifin = 1 DO jfl = 1, jpnfl IF( zagefl(jfl) < rdt ) ifin = 0 tpifl(jfl) = zgifl(jfl) + 0.5 tpjfl(jfl) = zgjfl(jfl) + 0.5 END DO ELSE ifin = 1 DO jfl = 1, jpnfl IF( zagefl(jfl) < rdt ) ifin = 0 tpifl(jfl) = zgifl(jfl) + 0.5 tpjfl(jfl) = zgjfl(jfl) + 0.5 IF( nisobfl(jfl) == 1 ) tpkfl(jfl) = -(zgkfl(jfl)) END DO ENDIF !!Alexcadm IF (lwp) write(numout,*) '---------' !!Alexcadm IF (lwp) write(numout,*) 'before Erika:',tpifl(880),tpjfl(880), !!Alexcadm . tpkfl(880),zufl(880),zvfl(880),zwfl(880) !!Alexcadm IF (lwp) write(numout,*) 'first Erika:',tpifl(900),tpjfl(900), !!Alexcadm . tpkfl(900),zufl(900),zvfl(900),zwfl(900) !!Alexcadm IF (lwp) write(numout,*) 'last Erika:',tpifl(jpnfl),tpjfl(jpnfl), !!Alexcadm . tpkfl(jpnfl),zufl(jpnfl),zvfl(jpnfl),zwfl(jpnfl) IF( ifin == 0 ) THEN iloop = iloop + 1 GO TO 222 ENDIF END SUBROUTINE flo_blk # else !!---------------------------------------------------------------------- !! Default option Empty module !!---------------------------------------------------------------------- CONTAINS SUBROUTINE flo_blk ! Empty routine END SUBROUTINE flo_blk #endif !!====================================================================== END MODULE floblk