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 10679 for NEMO/branches/2019/fix_ticket2238_solution2/src/OCE/ICB/icbutl.F90 – NEMO

Ignore:
Timestamp:
2019-02-14T14:11:43+01:00 (5 years ago)
Author:
mathiot
Message:

branch for solution 2 of ticket #2238

Location:
NEMO/branches/2019/fix_ticket2238_solution2
Files:
1 edited
1 copied

Legend:

Unmodified
Added
Removed

  • NEMO/branches/2019/fix_ticket2238_solution2/src/OCE/ICB/icbutl.F90

    r10570 r10679  
    3131   PRIVATE 
    3232 
    33    PUBLIC   icb_utl_copy          ! routine called in icbstp module 
    3433   PUBLIC   icb_utl_interp        ! routine called in icbdyn, icbthm modules 
    3534   PUBLIC   icb_utl_bilin         ! routine called in icbini, icbdyn modules 
    … …  
    5453CONTAINS 
    5554 
    56    SUBROUTINE icb_utl_copy() 
    57       !!---------------------------------------------------------------------- 
    58       !!                  ***  ROUTINE icb_utl_copy  *** 
    59       !! 
    60       !! ** Purpose :   iceberg initialization. 
    61       !! 
    62       !! ** Method  : - blah blah 
    63       !!---------------------------------------------------------------------- 
    64 #if defined key_si3 
    65       REAL(wp), DIMENSION(jpi,jpj) :: zssh_lead_m    !    ocean surface (ssh_m) if ice is not embedded 
    66       !                                              !    ocean surface in leads if ice is embedded    
    67 #endif 
    68       ! copy nemo forcing arrays into iceberg versions with extra halo 
    69       ! only necessary for variables not on T points 
    70       ! and ssh which is used to calculate gradients 
    71  
    72       uo_e(:,:) = 0._wp   ;   uo_e(1:jpi,1:jpj) = ssu_m(:,:) * umask(:,:,1) 
    73       vo_e(:,:) = 0._wp   ;   vo_e(1:jpi,1:jpj) = ssv_m(:,:) * vmask(:,:,1) 
    74       ff_e(:,:) = 0._wp   ;   ff_e(1:jpi,1:jpj) = ff_f (:,:)  
    75       tt_e(:,:) = 0._wp   ;   tt_e(1:jpi,1:jpj) = sst_m(:,:) 
    76       fr_e(:,:) = 0._wp   ;   fr_e(1:jpi,1:jpj) = fr_i (:,:) 
    77       ua_e(:,:) = 0._wp   ;   ua_e(1:jpi,1:jpj) = utau (:,:) * umask(:,:,1) ! maybe mask useless because mask applied in sbcblk 
    78       va_e(:,:) = 0._wp   ;   va_e(1:jpi,1:jpj) = vtau (:,:) * vmask(:,:,1) ! maybe mask useless because mask applied in sbcblk 
    79       ! 
    80       CALL lbc_lnk_icb( 'icbutl', uo_e, 'U', -1._wp, 1, 1 ) 
    81       CALL lbc_lnk_icb( 'icbutl', vo_e, 'V', -1._wp, 1, 1 ) 
    82       CALL lbc_lnk_icb( 'icbutl', ff_e, 'F', +1._wp, 1, 1 ) 
    83       CALL lbc_lnk_icb( 'icbutl', ua_e, 'U', -1._wp, 1, 1 ) 
    84       CALL lbc_lnk_icb( 'icbutl', va_e, 'V', -1._wp, 1, 1 ) 
    85       CALL lbc_lnk_icb( 'icbutl', fr_e, 'T', +1._wp, 1, 1 ) 
    86       CALL lbc_lnk_icb( 'icbutl', tt_e, 'T', +1._wp, 1, 1 ) 
    87 #if defined key_si3 
    88       hicth(:,:) = 0._wp ;  hicth(1:jpi,1:jpj) = hm_i (:,:)   
    89       ui_e(:,:) = 0._wp ;   ui_e(1:jpi, 1:jpj) = u_ice(:,:) 
    90       vi_e(:,:) = 0._wp ;   vi_e(1:jpi, 1:jpj) = v_ice(:,:) 
    91       !       
    92       ! compute ssh slope using ssh_lead if embedded 
    93       zssh_lead_m(:,:) = ice_var_sshdyn(ssh_m, snwice_mass, snwice_mass_b) 
    94       ssh_e(:,:) = 0._wp ;  ssh_e(1:jpi, 1:jpj) = zssh_lead_m(:,:) * tmask(:,:,1) 
    95       ! 
    96       CALL lbc_lnk_icb( 'icbutl', hicth, 'T', +1._wp, 1, 1 ) 
    97       CALL lbc_lnk_icb( 'icbutl', ui_e , 'U', -1._wp, 1, 1 ) 
    98       CALL lbc_lnk_icb( 'icbutl', vi_e , 'V', -1._wp, 1, 1 ) 
    99 #else 
    100       ssh_e(:,:) = 0._wp ;  ssh_e(1:jpi, 1:jpj) = ssh_m(:,:) * tmask(:,:,1) 
    101 #endif 
    102  
    103       !! special for ssh which is used to calculate slope 
    104       !! so fudge some numbers all the way around the boundary 
    105       ssh_e(0    ,    :) = ssh_e(1  ,  :) 
    106       ssh_e(jpi+1,    :) = ssh_e(jpi,  :) 
    107       ssh_e(:    ,    0) = ssh_e(:  ,  1) 
    108       ssh_e(:    ,jpj+1) = ssh_e(:  ,jpj) 
    109       ssh_e(0,0)         = ssh_e(1,1) 
    110       ssh_e(jpi+1,0)     = ssh_e(jpi,1) 
    111       ssh_e(0,jpj+1)     = ssh_e(1,jpj) 
    112       ssh_e(jpi+1,jpj+1) = ssh_e(jpi,jpj) 
    113       CALL lbc_lnk_icb( 'icbutl', ssh_e, 'T', +1._wp, 1, 1 ) 
    114       ! 
    115    END SUBROUTINE icb_utl_copy 
    116  
    117  
    11855   SUBROUTINE icb_utl_interp( pi, pe1, puo, pui, pua, pssh_i,   & 
    11956      &                       pj, pe2, pvo, pvi, pva, pssh_j,   & 
    … …  
    14077      ! 
    14178      REAL(wp) ::   zcd, zmod       ! local scalars 
     79      REAL(wp), DIMENSION(jpi,jpj) :: zssh 
    14280      !!---------------------------------------------------------------------- 
    14381 
    … …  
    14583      pe2 = icb_utl_bilin_e( e2t, e2u, e2v, e2f, pi, pj ) 
    14684      ! 
    147       puo  = icb_utl_bilin_h( uo_e, pi, pj, 'U' )             ! ocean velocities 
    148       pvo  = icb_utl_bilin_h( vo_e, pi, pj, 'V' ) 
    149       psst = icb_utl_bilin_h( tt_e, pi, pj, 'T' )             ! SST 
    150       pcn  = icb_utl_bilin_h( fr_e , pi, pj, 'T' )            ! ice concentration 
    151       pff  = icb_utl_bilin_h( ff_e , pi, pj, 'F' )            ! Coriolis parameter 
    152       ! 
    153       pua  = icb_utl_bilin_h( ua_e , pi, pj, 'U' )            ! 10m wind 
    154       pva  = icb_utl_bilin_h( va_e , pi, pj, 'V' )            ! here (ua,va) are stress => rough conversion from stress to speed 
     85      puo  = icb_utl_bilin( ssu_m(:,:) * umask(:,:,1), pi, pj, 'U' )             ! ocean velocities 
     86      pvo  = icb_utl_bilin( ssv_m(:,:) * vmask(:,:,1), pi, pj, 'V' ) 
     87      psst = icb_utl_bilin( sst_m(:,:), pi, pj, 'T' )             ! SST 
     88      pcn  = icb_utl_bilin( fr_i(:,:) , pi, pj, 'T' )            ! ice concentration 
     89      pff  = icb_utl_bilin( ff_f(:,:) , pi, pj, 'F' )            ! Coriolis parameter 
     90      ! 
     91      pua  = icb_utl_bilin( utau(:,:) * umask(:,:,1), pi, pj, 'U' )            ! 10m wind 
     92      pva  = icb_utl_bilin( vtau(:,:) * vmask(:,:,1), pi, pj, 'V' )            ! here (ua,va) are stress => rough conversion from stress to speed 
    15593      zcd  = 1.22_wp * 1.5e-3_wp                              ! air density * drag coefficient 
    15694      zmod = 1._wp / MAX(  1.e-20, SQRT(  zcd * SQRT( pua*pua + pva*pva)  )  ) 
    … …  
    15997 
    16098#if defined key_si3 
    161       pui = icb_utl_bilin_h( ui_e , pi, pj, 'U' )              ! sea-ice velocities 
    162       pvi = icb_utl_bilin_h( vi_e , pi, pj, 'V' ) 
    163       phi = icb_utl_bilin_h( hicth, pi, pj, 'T' )              ! ice thickness 
     99      pui = icb_utl_bilin( u_ice , pi, pj, 'U' )              ! sea-ice velocities 
     100      pvi = icb_utl_bilin( v_ice , pi, pj, 'V' ) 
     101      phi = icb_utl_bilin( hm_i  , pi, pj, 'T' )              ! ice thickness 
     102      ! Estimate SSH gradient in i- and j-direction (centred evaluation) 
     103      ! compute ssh slope using ssh_lead if embedded 
     104      zssh(:,:) = ice_var_sshdyn(ssh_m, snwice_mass, snwice_mass_b) 
    164105#else 
    165106      pui = 0._wp 
    166107      pvi = 0._wp 
    167108      phi = 0._wp 
     109      zssh(:,:) = ssh_m(:,:) 
    168110#endif 
     111      zssh(:,:) = zssh(:,:) * tmask(:,:,1) 
    169112 
    170113      ! Estimate SSH gradient in i- and j-direction (centred evaluation) 
    171       pssh_i = ( icb_utl_bilin_h( ssh_e, pi+0.1_wp, pj, 'T' ) -   & 
    172          &       icb_utl_bilin_h( ssh_e, pi-0.1_wp, pj, 'T' )  ) / ( 0.2_wp * pe1 ) 
    173       pssh_j = ( icb_utl_bilin_h( ssh_e, pi, pj+0.1_wp, 'T' ) -   & 
    174          &       icb_utl_bilin_h( ssh_e, pi, pj-0.1_wp, 'T' )  ) / ( 0.2_wp * pe2 ) 
     114      pssh_i = ( icb_utl_bilin( zssh(:,:), pi+0.1_wp, pj, 'T' ) -   & 
     115         &       icb_utl_bilin( zssh(:,:), pi-0.1_wp, pj, 'T' )  ) / ( 0.2_wp * pe1 ) 
     116      pssh_j = ( icb_utl_bilin( zssh(:,:), pi, pj+0.1_wp, 'T' ) -   & 
     117         &       icb_utl_bilin( zssh(:,:), pi, pj-0.1_wp, 'T' )  ) / ( 0.2_wp * pe2 ) 
    175118      ! 
    176119   END SUBROUTINE icb_utl_interp 
    177  
    178  
    179    REAL(wp) FUNCTION icb_utl_bilin_h( pfld, pi, pj, cd_type ) 
    180       !!---------------------------------------------------------------------- 
    181       !!                  ***  FUNCTION icb_utl_bilin  *** 
    182       !! 
    183       !! ** Purpose :   bilinear interpolation at berg location depending on the grid-point type 
    184       !!                this version deals with extra halo points 
    185       !! 
    186       !!       !!gm  CAUTION an optional argument should be added to handle 
    187       !!             the slip/no-slip conditions  ==>>> to be done later 
    188       !! 
    189       !!---------------------------------------------------------------------- 
    190       REAL(wp), DIMENSION(0:jpi+1,0:jpj+1), INTENT(in) ::   pfld      ! field to be interpolated 
    191       REAL(wp)                            , INTENT(in) ::   pi, pj    ! targeted coordinates in (i,j) referential 
    192       CHARACTER(len=1)                    , INTENT(in) ::   cd_type   ! type of pfld array grid-points: = T , U , V or F points 
    193       ! 
    194       INTEGER  ::   ii, ij   ! local integer 
    195       REAL(wp) ::   zi, zj   ! local real 
    196       !!---------------------------------------------------------------------- 
    197       ! 
    198       SELECT CASE ( cd_type ) 
    199       CASE ( 'T' ) 
    200          ! note that here there is no +0.5 added 
    201          ! since we're looking for four T points containing quadrant we're in of  
    202          ! current T cell 
    203          ii = MAX(1, INT( pi     )) 
    204          ij = MAX(1, INT( pj     ))    ! T-point 
    205          zi = pi - REAL(ii,wp) 
    206          zj = pj - REAL(ij,wp) 
    207       CASE ( 'U' ) 
    208          ii = MAX(1, INT( pi-0.5 )) 
    209          ij = MAX(1, INT( pj     ))    ! U-point 
    210          zi = pi - 0.5 - REAL(ii,wp) 
    211          zj = pj - REAL(ij,wp) 
    212       CASE ( 'V' ) 
    213          ii = MAX(1, INT( pi     )) 
    214          ij = MAX(1, INT( pj-0.5 ))    ! V-point 
    215          zi = pi - REAL(ii,wp) 
    216          zj = pj - 0.5 - REAL(ij,wp) 
    217       CASE ( 'F' ) 
    218          ii = MAX(1, INT( pi-0.5 )) 
    219          ij = MAX(1, INT( pj-0.5 ))    ! F-point 
    220          zi = pi - 0.5 - REAL(ii,wp) 
    221          zj = pj - 0.5 - REAL(ij,wp) 
    222       END SELECT 
    223       ! 
    224       ! find position in this processor. Prevent near edge problems (see #1389) 
    225       ! 
    226       IF    ( ii < mig( 1 ) ) THEN   ;   ii = 1 
    227       ELSEIF( ii > mig(jpi) ) THEN   ;   ii = jpi 
    228       ELSE                           ;   ii = mi1(ii) 
    229       ENDIF 
    230       IF    ( ij < mjg( 1 ) ) THEN   ;   ij = 1 
    231       ELSEIF( ij > mjg(jpj) ) THEN   ;   ij = jpj 
    232       ELSE                           ;   ij  = mj1(ij) 
    233       ENDIF 
    234       ! 
    235       IF( ii == jpi )   ii = ii-1       
    236       IF( ij == jpj )   ij = ij-1 
    237       ! 
    238       icb_utl_bilin_h = ( pfld(ii,ij  ) * (1.-zi) + pfld(ii+1,ij  ) * zi ) * (1.-zj)   & 
    239          &            + ( pfld(ii,ij+1) * (1.-zi) + pfld(ii+1,ij+1) * zi ) *     zj 
    240       ! 
    241    END FUNCTION icb_utl_bilin_h 
    242  
    243120 
    244121   REAL(wp) FUNCTION icb_utl_bilin( pfld, pi, pj, cd_type ) 
    … …  
    270147            zj = pj - REAL(ij,wp) 
    271148         CASE ( 'U' ) 
    272             ii = MAX(1, INT( pi-0.5 )) 
    273             ij = MAX(1, INT( pj     ))    ! U-point 
    274             zi = pi - 0.5 - REAL(ii,wp) 
     149            ii = MAX(1, INT( pi-0.5_wp )) 
     150            ij = MAX(1, INT( pj        ))    ! U-point 
     151            zi = pi - 0.5_wp - REAL(ii,wp) 
    275152            zj = pj - REAL(ij,wp) 
    276153         CASE ( 'V' ) 
    277             ii = MAX(1, INT( pi     )) 
    278             ij = MAX(1, INT( pj-0.5 ))    ! V-point 
     154            ii = MAX(1, INT( pi        )) 
     155            ij = MAX(1, INT( pj-0.5_wp ))    ! V-point 
    279156            zi = pi - REAL(ii,wp) 
    280             zj = pj - 0.5 - REAL(ij,wp) 
     157            zj = pj - 0.5_wp - REAL(ij,wp) 
    281158         CASE ( 'F' ) 
    282             ii = MAX(1, INT( pi-0.5 )) 
    283             ij = MAX(1, INT( pj-0.5 ))    ! F-point 
    284             zi = pi - 0.5 - REAL(ii,wp) 
    285             zj = pj - 0.5 - REAL(ij,wp) 
     159            ii = MAX(1, INT( pi-0.5_wp )) 
     160            ij = MAX(1, INT( pj-0.5_wp ))    ! F-point 
     161            zi = pi - 0.5_wp - REAL(ii,wp) 
     162            zj = pj - 0.5_wp - REAL(ij,wp) 
    286163      END SELECT 
    287164      ! 
    288       ! find position in this processor. Prevent near edge problems (see #1389) 
    289       IF    ( ii < mig( 1 ) ) THEN   ;   ii = 1 
    290       ELSEIF( ii > mig(jpi) ) THEN   ;   ii = jpi 
    291       ELSE                           ;   ii = mi1(ii) 
    292       ENDIF 
    293       IF    ( ij < mjg( 1 ) ) THEN   ;   ij = 1 
    294       ELSEIF( ij > mjg(jpj) ) THEN   ;   ij = jpj 
    295       ELSE                           ;   ij  = mj1(ij) 
    296       ENDIF 
    297       ! 
    298       IF( ii == jpi )   ii = ii-1       
    299       IF( ij == jpj )   ij = ij-1 
     165      ! find position in this processor. 
     166      ii = mi1(ii) ; ij  = mj1(ij) 
    300167      ! 
    301168      icb_utl_bilin = ( pfld(ii,ij  ) * (1.-zi) + pfld(ii+1,ij  ) * zi ) * (1.-zj)   & 
    … …  
    333200      zj = pj - REAL(ij,wp) 
    334201      ! 
    335       ! find position in this processor. Prevent near edge problems (see #1389) 
    336       IF    ( ii < mig( 1 ) ) THEN   ;   ii = 1 
    337       ELSEIF( ii > mig(jpi) ) THEN   ;   ii = jpi 
    338       ELSE                           ;   ii = mi1(ii) 
    339       ENDIF 
    340       IF    ( ij < mjg( 1 ) ) THEN   ;   ij = 1 
    341       ELSEIF( ij > mjg(jpj) ) THEN   ;   ij = jpj 
    342       ELSE                           ;   ij  = mj1(ij) 
    343       ENDIF 
    344       ! 
    345       IF( ii == jpi )   ii = ii-1       
    346       IF( ij == jpj )   ij = ij-1 
     202      ! find position in this processor. 
     203      ii = mi1(ii) ; ij  = mj1(ij) 
    347204      ! 
    348205      z4(1) = pfld(ii  ,ij  ) 
    … …  
    392249      zi = pi - REAL(ii,wp)          !!gm use here mig, mjg arrays 
    393250      zj = pj - REAL(ij,wp) 
    394  
    395       ! find position in this processor. Prevent near edge problems (see #1389) 
    396       IF    ( ii < mig( 1 ) ) THEN   ;   ii = 1 
    397       ELSEIF( ii > mig(jpi) ) THEN   ;   ii = jpi 
    398       ELSE                           ;   ii = mi1(ii) 
    399       ENDIF 
    400       IF    ( ij < mjg( 1 ) ) THEN   ;   ij = 1 
    401       ELSEIF( ij > mjg(jpj) ) THEN   ;   ij = jpj 
    402       ELSE                           ;   ij  = mj1(ij) 
    403       ENDIF 
    404       ! 
    405       IF( ii == jpi )   ii = ii-1       
    406       IF( ij == jpj )   ij = ij-1 
     251      ! 
     252      ! find position in this processor. 
     253      ii = mi1(ii) ; ij  = mj1(ij) 
    407254      ! 
    408255      IF(    0.0_wp <= zi .AND. zi < 0.5_wp   ) THEN 
Note: See TracChangeset for help on using the changeset viewer.