Changeset 13357


Ignore:
Timestamp:
2020-07-30T13:20:57+02:00 (12 days ago)
Author:
mathiot
Message:

ticket #1900: add changes for Merino 2016 works. Results unchanged if flag ln_M2016 is set to .false. . Remaining step: test ln_M2016 = .true.

Location:
NEMO/branches/2020/tickets_icb_1900
Files:
10 edited

Legend:

Unmodified
Added
Removed
  • NEMO/branches/2020/tickets_icb_1900/cfgs/SHARED/namelist_ref

    r13216 r13357  
    614614   ln_use_calving          = .false. ! Use calving data even when nn_test_icebergs > 0 
    615615   rn_speed_limit          = 0.      ! CFL speed limit for a berg 
    616  
     616   ! 
     617   ln_M2016                = .false. ! use Merino et al. (2016) modification (use of 3d ocean data instead of only sea surface data) 
     618   ! 
    617619   cn_dir      = './'      !  root directory for the calving data location 
    618620   !___________!_________________________!___________________!___________!_____________!________!___________!__________________!__________!_______________! 
  • NEMO/branches/2020/tickets_icb_1900/src/OCE/ICB/icb_oce.F90

    r13265 r13357  
    5757   TYPE, PUBLIC ::   point              !: properties of an individual iceberg (position, mass, size, etc...) 
    5858      INTEGER  ::   year 
    59       REAL(wp) ::   xi , yj                                              ! iceberg coordinates in the (i,j) referential (global) 
     59      REAL(wp) ::   xi , yj , zk                                         ! iceberg coordinates in the (i,j) referential (global) and deepest level affected 
    6060      REAL(wp) ::   e1 , e2                                              ! horizontal scale factors at the iceberg position 
    6161      REAL(wp) ::   lon, lat, day                                        ! geographic position 
    6262      REAL(wp) ::   mass, thickness, width, length, uvel, vvel           ! iceberg physical properties 
    63       REAL(wp) ::   uo, vo, ui, vi, ua, va, ssh_x, ssh_y, sst, cn, hi    ! properties of iceberg environment  
     63      REAL(wp) ::   ssu, ssv, ui, vi, ua, va, ssh_x, ssh_y, sst, cn, hi  ! properties of iceberg environment  
    6464      REAL(wp) ::   mass_of_bits, heat_density 
     65      INTEGER  ::   kb                                                   ! icb bottom level 
    6566   END TYPE point 
    6667 
     
    8586   ! Extra arrays with bigger halo, needed when interpolating forcing onto iceberg position 
    8687   ! particularly for MPP when iceberg can lie inside T grid but outside U, V, or f grid 
    87    REAL(wp), PUBLIC, DIMENSION(:,:), ALLOCATABLE ::   uo_e, vo_e 
    88    REAL(wp), PUBLIC, DIMENSION(:,:), ALLOCATABLE ::   tt_e, fr_e 
     88   REAL(wp), PUBLIC, DIMENSION(:,:), ALLOCATABLE ::   ssu_e, ssv_e 
     89   REAL(wp), PUBLIC, DIMENSION(:,:), ALLOCATABLE ::   sst_e, fr_e 
    8990   REAL(wp), PUBLIC, DIMENSION(:,:), ALLOCATABLE ::   ua_e, va_e 
    9091   REAL(wp), PUBLIC, DIMENSION(:,:), ALLOCATABLE ::   ssh_e 
    9192   REAL(wp), PUBLIC, DIMENSION(:,:), ALLOCATABLE ::   tmask_e, umask_e, vmask_e 
    9293   REAl(wp), PUBLIC, DIMENSION(:,:), ALLOCATABLE ::   rlon_e, rlat_e, ff_e 
    93       ! prepration to Nacho work 
    94       ! tt3d_e 
    95       ! uo3d_e 
    96       ! vo3d_e 
    97       ! 
     94   REAl(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE ::   uoce_e, voce_e, toce_e, e3t_e 
     95   ! 
    9896#if defined key_si3 || defined key_cice 
    9997   REAL(wp), PUBLIC, DIMENSION(:,:), ALLOCATABLE ::   hi_e, ui_e, vi_e 
     
    123121   INTEGER , PUBLIC ::   nn_verbose_write                !: timesteps between verbose messages 
    124122   REAL(wp), PUBLIC ::   rn_rho_bergs                    !: Density of icebergs 
     123   REAL(wp), PUBLIC ::   rho_berg_1_oce                  !: convertion factor (thickness to draft) (rn_rho_bergs/pp_rho_seawater) 
    125124   REAL(wp), PUBLIC ::   rn_LoW_ratio                    !: Initial ratio L/W for newly calved icebergs 
    126125   REAL(wp), PUBLIC ::   rn_bits_erosion_fraction        !: Fraction of erosion melt flux to divert to bergy bits 
     
    130129   LOGICAL , PUBLIC ::   ln_time_average_weight          !: Time average the weight on the ocean    !!gm I don't understand that ! 
    131130   REAL(wp), PUBLIC ::   rn_speed_limit                  !: CFL speed limit for a berg 
     131   LOGICAL , PUBLIC ::   ln_M2016                        !: use Nacho's Merino 2016 work 
    132132   ! 
    133133   ! restart 
     
    141141   REAL(wp), DIMENSION(nclasses), PUBLIC ::   rn_initial_thickness !  Single instance of an icebergs type initialised in icebergs_init and updated in icebergs_run 
    142142   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   src_calving, src_calving_hflx    !: accumulate input ice 
     143   INTEGER , PUBLIC             , SAVE                     ::   micbkb                           !: deepest level affected by icebergs 
    143144   INTEGER , PUBLIC             , SAVE                     ::   numicb                           !: iceberg IO 
    144145   INTEGER , PUBLIC             , SAVE, DIMENSION(nkounts) ::   num_bergs                        !: iceberg counter 
     
    180181      ! 
    181182      ! expanded arrays for bilinear interpolation 
    182       ALLOCATE( uo_e(0:jpi+1,0:jpj+1) , ua_e(0:jpi+1,0:jpj+1) ,    & 
    183          &      vo_e(0:jpi+1,0:jpj+1) , va_e(0:jpi+1,0:jpj+1) ,    & 
     183      ALLOCATE( ssu_e(0:jpi+1,0:jpj+1) , ua_e(0:jpi+1,0:jpj+1) ,   & 
     184         &      ssv_e(0:jpi+1,0:jpj+1) , va_e(0:jpi+1,0:jpj+1) ,   & 
    184185#if defined key_si3 || defined key_cice 
    185186         &      ui_e(0:jpi+1,0:jpj+1) ,                            & 
     
    188189#endif 
    189190         &      fr_e(0:jpi+1,0:jpj+1) ,                            & 
    190          &      tt_e(0:jpi+1,0:jpj+1) , ssh_e(0:jpi+1,0:jpj+1) ,   & 
     191         &      sst_e(0:jpi+1,0:jpj+1) , ssh_e(0:jpi+1,0:jpj+1) ,  & 
    191192         &      first_width(nclasses) , first_length(nclasses) ,   & 
    192193         &      src_calving (jpi,jpj) ,                            & 
     
    194195      icb_alloc = icb_alloc + ill 
    195196 
    196       ! prepration to Nacho work 
    197       ! tt3d_e 
    198       ! uo3d_e 
    199       ! vo3d_e 
     197      IF ( ln_M2016 ) THEN 
     198         ALLOCATE( uoce_e(0:jpi+1,0:jpj+1,jpk), voce_e(0:jpi+1,0:jpj+1,jpk), & 
     199            &      toce_e(0:jpi+1,0:jpj+1,jpk), e3t_e(0:jpi+1,0:jpj+1,jpk) , STAT=ill ) 
     200         icb_alloc = icb_alloc + ill 
     201      END IF 
    200202      ! 
    201203      ALLOCATE( tmask_e(0:jpi+1,0:jpj+1), umask_e(0:jpi+1,0:jpj+1), vmask_e(0:jpi+1,0:jpj+1), & 
  • NEMO/branches/2020/tickets_icb_1900/src/OCE/ICB/icbclv.F90

    r13273 r13357  
    2121   USE lbclnk         ! NEMO boundary exchanges for gridded data 
    2222 
    23    USE icb_oce        ! iceberg variables 
    2423   USE icbdia         ! iceberg diagnostics 
    2524   USE icbutl         ! iceberg utility routines 
     
    142141                  newpt%mass           = rn_initial_mass     (jn) 
    143142                  newpt%thickness      = rn_initial_thickness(jn) 
     143                  newpt%kb             = 1         ! compute correctly in icbthm if needed         
    144144                  newpt%width          = first_width         (jn) 
    145145                  newpt%length         = first_length        (jn) 
  • NEMO/branches/2020/tickets_icb_1900/src/OCE/ICB/icbdyn.F90

    r13265 r13357  
    197197      ij  = mj1( ij  ) 
    198198      ! 
     199      ! assume icb is grounded if tmask(ii,ij,1) or tmask(ii,ij,ikb), depending of the option is not 0 
     200      !IF ( ln_icb_ground ) THEN 
     201      !   ! interpol needed data 
     202      !   CALL icb_utl_interp( pxi, pyj, pe3t=ze3t )   ! 3d velocities 
     203      !   
     204      !   !compute bottom level 
     205      !   CALL icb_utl_getkb( ikb, ze3t, zD ) 
     206      ! 
     207      !   IF(  tmask(ii,ij,ikb)  /=   0._wp  )   RETURN           ! berg reach a new t-cell, but an ocean one 
     208      !ELSE 
    199209      IF(  tmask(ii,ij,1)  /=   0._wp  )   RETURN           ! berg reach a new t-cell, but an ocean one 
     210      !END IF 
    200211      ! 
    201212      ! From here, berg have reach land: treat grounding/bouncing 
     
    254265      REAL(wp), PARAMETER ::   pp_Cr0       = 0.06_wp    ! 
    255266      ! 
    256       INTEGER  ::   itloop 
    257       REAL(wp) ::   zuo, zui, zua, zuwave, zssh_x, zsst, zcn, zhi 
    258       REAL(wp) ::   zvo, zvi, zva, zvwave, zssh_y 
     267      INTEGER  ::   itloop, ikb, jk 
     268      REAL(wp) ::   zuo, zssu, zui, zua, zuwave, zssh_x, zcn, zhi 
     269      REAL(wp) ::   zvo, zssv, zvi, zva, zvwave, zssh_y 
    259270      REAL(wp) ::   zff, zT, zD, zW, zL, zM, zF 
    260271      REAL(wp) ::   zdrag_ocn, zdrag_atm, zdrag_ice, zwave_rad 
    261       REAL(wp) ::   z_ocn, z_atm, z_ice 
     272      REAL(wp) ::   z_ocn, z_atm, z_ice, zdep 
    262273      REAL(wp) ::   zampl, zwmod, zCr, zLwavelength, zLcutoff, zLtop 
    263274      REAL(wp) ::   zlambda, zdetA, zA11, zA12, zaxe, zaye, zD_hi 
    264275      REAL(wp) ::   zuveln, zvveln, zus, zvs, zspeed, zloc_dx, zspeed_new 
     276      REAL(wp), DIMENSION(jpk) :: zuoce, zvoce, ze3t, zdepw 
    265277      !!---------------------------------------------------------------------- 
    266278 
    267279      ! Interpolate gridded fields to berg 
    268280      nknberg = berg%number(1) 
    269       CALL icb_utl_interp( pxi, pyj, pe1=pe1, pe2=pe2,   &   ! scale factor 
    270          &                 puo=zuo, pui=zui, pua=zua,    &   ! oce/ice/atm velocities 
    271          &                 pvo=zvo, pvi=zvi, pva=zva,    &   ! oce/ice/atm velocities 
    272          &                 pssh_i=zssh_x, pssh_j=zssh_y, &   ! ssh gradient 
    273          &                 phi=zhi, pff=zff )               ! ice thickness and coriolis 
    274  
    275       IF (lwp) THEN 
    276          WRITE(numout,*) 'icbdyn ', pxi, pyj  
    277          WRITE(numout,*) pe1, zuo, zui, zua, zssh_x 
    278          WRITE(numout,*) pe2, zvo, zvi, zva, zssh_y 
    279          WRITE(numout,*) zhi, zff  
    280          WRITE(numout,*) '' 
    281       END IF 
    282  
     281      CALL icb_utl_interp( pxi, pyj, pe1=pe1, pe2=pe2,     &   ! scale factor 
     282         &                 pssu=zssu, pui=zui, pua=zua,    &   ! oce/ice/atm velocities 
     283         &                 pssv=zssv, pvi=zvi, pva=zva,    &   ! oce/ice/atm velocities 
     284         &                 pssh_i=zssh_x, pssh_j=zssh_y,   &   ! ssh gradient 
     285         &                 phi=zhi, pff=zff)                   ! ice thickness and coriolis 
    283286 
    284287      zM = berg%current_point%mass 
    285288      zT = berg%current_point%thickness               ! total thickness 
    286       zD = ( rn_rho_bergs / pp_rho_seawater ) * zT    ! draught (keel depth) 
     289      zD = rho_berg_1_oce * zT                        ! draught (keel depth) 
    287290      zF = zT - zD                                    ! freeboard 
    288291      zW = berg%current_point%width 
     
    291294      zhi   = MIN( zhi   , zD    ) 
    292295      zD_hi = MAX( 0._wp, zD-zhi ) 
    293  
    294       ! Wave radiation 
    295       zuwave = zua - zuo   ;   zvwave = zva - zvo     ! Use wind speed rel. to ocean for wave model 
     296  
     297     ! Wave radiation 
     298      zuwave = zua - zssu   ;   zvwave = zva - zssv   ! Use wind speed rel. to ocean for wave model 
    296299      zwmod  = zuwave*zuwave + zvwave*zvwave          ! The wave amplitude and length depend on the  current; 
    297300      !                                               ! wind speed relative to the ocean. Actually wmod is wmod**2 here. 
     
    318321      IF( abs(zui) + abs(zvi) == 0._wp )   z_ice = 0._wp 
    319322 
     323      ! lateral velocities 
     324      ! default ssu and ssv 
     325      ! ln_M2016: mean velocity along the profile 
     326      IF ( ln_M2016 ) THEN 
     327         ! interpol needed data 
     328         CALL icb_utl_interp( pxi, pyj, puoce=zuoce, pvoce=zvoce, pe3t=ze3t )   ! 3d velocities 
     329         
     330         !compute bottom level 
     331         CALL icb_utl_getkb( ikb, ze3t, zD ) 
     332          
     333         ! compute mean velocity  
     334         zuo = 0.0 ; zvo = 0.0; zdep = 0.0 
     335         DO jk = 1, ikb-1 
     336            zuo = zuo + zuoce(jk) * ze3t(jk) 
     337            zvo = zvo + zvoce(jk) * ze3t(jk) 
     338            zdep = zdep + ze3t(jk) 
     339         END DO 
     340         zuo = (zuo + zuoce(ikb) * (zD - zdep)) / zD 
     341         zvo = (zvo + zvoce(ikb) * (zD - zdep)) / zD 
     342      ELSE 
     343         zuo = zssu 
     344         zvo = zssv 
     345      END IF 
     346 
    320347      zuveln = puvel   ;   zvveln = pvvel ! Copy starting uvel, vvel 
    321348      ! 
     
    330357         ! Explicit accelerations 
    331358         !zaxe= zff*pvvel -grav*zssh_x +zwave_rad*zuwave & 
    332          !    -zdrag_ocn*(puvel-zuo) -zdrag_atm*(puvel-zua) -zdrag_ice*(puvel-zui) 
     359         !    -zdrag_ocn*(puvel-zssu) -zdrag_atm*(puvel-zua) -zdrag_ice*(puvel-zui) 
    333360         !zaye=-zff*puvel -grav*zssh_y +zwave_rad*zvwave & 
    334          !    -zdrag_ocn*(pvvel-zvo) -zdrag_atm*(pvvel-zva) -zdrag_ice*(pvvel-zvi) 
     361         !    -zdrag_ocn*(pvvel-zssv) -zdrag_atm*(pvvel-zva) -zdrag_ice*(pvvel-zvi) 
    335362         zaxe = -grav * zssh_x + zwave_rad * zuwave 
    336363         zaye = -grav * zssh_y + zwave_rad * zvwave 
  • NEMO/branches/2020/tickets_icb_1900/src/OCE/ICB/icbini.F90

    r13276 r13357  
    7878      ! 
    7979      !                          ! initialised variable with extra haloes to zero 
    80       uo_e(:,:) = 0._wp   ;   vo_e(:,:) = 0._wp   ; 
    81       ua_e(:,:) = 0._wp   ;   va_e(:,:) = 0._wp   ; 
    82       ff_e(:,:) = 0._wp   ;   tt_e(:,:) = 0._wp   ; 
    83       fr_e(:,:) = 0._wp   ; 
    84       ! prepration to Nacho work 
    85       ! e3t_e 
    86       ! tt3d_e 
    87       ! uo3d_e 
    88       ! vo3d_e 
     80      ssu_e(:,:) = 0._wp   ;   ssv_e(:,:) = 0._wp   ; 
     81      ua_e(:,:)  = 0._wp   ;   va_e(:,:)  = 0._wp   ; 
     82      ff_e(:,:)  = 0._wp   ;   sst_e(:,:)  = 0._wp   ; 
     83      fr_e(:,:)  = 0._wp   ; 
     84      ! 
     85      IF ( ln_M2016 ) THEN 
     86         toce_e(:,:,:) = 0._wp 
     87         uoce_e(:,:,:) = 0._wp 
     88         voce_e(:,:,:) = 0._wp 
     89         e3t_e(:,:,:)  = 0._wp 
     90      END IF 
    8991      ! 
    9092#if defined key_si3 
     
    106108      first_width (:) = SQRT(  rn_initial_mass(:) / ( rn_LoW_ratio * rn_rho_bergs * rn_initial_thickness(:) )  ) 
    107109      first_length(:) = rn_LoW_ratio * first_width(:) 
     110      rho_berg_1_oce  = rn_rho_bergs / pp_rho_seawater  ! scale factor used for convertion thickness to draft 
     111      ! 
     112      ! deepest level affected by icebergs 
     113      ! can be tuned but the safest is this  
     114      ! (with z* and z~ the depth of each level change overtime, so the more robust micbkb is jpk) 
     115      micbkb = jpk 
    108116 
    109117      berg_grid%calving      (:,:)   = 0._wp 
     
    365373               localpt%uvel = 0._wp 
    366374               localpt%vvel = 0._wp 
     375               localpt%kb   = 1 
    367376               CALL icb_utl_incr() 
    368377               localberg%number(:) = num_bergs(:) 
     
    398407         &              rn_bits_erosion_fraction        , rn_sicn_shift       , ln_passive_mode      ,   & 
    399408         &              ln_time_average_weight          , nn_test_icebergs    , rn_test_box          ,   & 
    400          &              ln_use_calving , rn_speed_limit , cn_dir, sn_icb      ,                          & 
     409         &              ln_use_calving , rn_speed_limit , cn_dir, sn_icb      , ln_M2016             ,   & 
    401410         &              cn_icbrst_indir, cn_icbrst_in   , cn_icbrst_outdir    , cn_icbrst_out 
    402411      !!---------------------------------------------------------------------- 
  • NEMO/branches/2020/tickets_icb_1900/src/OCE/ICB/icbstp.F90

    r11536 r13357  
    5252CONTAINS 
    5353 
    54    SUBROUTINE icb_stp( kt ) 
     54   SUBROUTINE icb_stp( kt, Kmm ) 
    5555      !!---------------------------------------------------------------------- 
    5656      !!                  ***  ROUTINE icb_stp  *** 
     
    6161      !!---------------------------------------------------------------------- 
    6262      INTEGER, INTENT(in) ::   kt   ! time step index 
     63      INTEGER, INTENT(in) ::   Kmm  ! ocean time level index 
    6364      ! 
    6465      LOGICAL ::   ll_sample_traj, ll_budget, ll_verbose   ! local logical 
     
    9293      ! 
    9394      !                                   !* copy nemo forcing arrays into iceberg versions with extra halo 
    94       CALL icb_utl_copy()                 ! only necessary for variables not on T points 
     95      CALL icb_utl_copy( Kmm )                 ! only necessary for variables not on T points 
    9596      ! 
    9697      ! 
  • NEMO/branches/2020/tickets_icb_1900/src/OCE/ICB/icbthm.F90

    r13265 r13357  
    4848      INTEGER, INTENT(in) ::   kt   ! timestep number, just passed to icb_utl_print_berg 
    4949      ! 
    50       INTEGER  ::   ii, ij 
    51       REAL(wp) ::   zM, zT, zW, zL, zSST, zVol, zLn, zWn, zTn, znVol, zIC, zDn 
    52       REAL(wp) ::   zMv, zMe, zMb, zmelt, zdvo, zdva, zdM, zSs, zdMe, zdMb, zdMv 
     50      INTEGER  ::   ii, ij, jk, ikb 
     51      REAL(wp) ::   zM, zT, zW, zL, zSST, zVol, zLn, zWn, zTn, znVol, zIC, zDn, zD, zvb, zub, ztb 
     52      REAL(wp) ::   zMv, zMe, zMb, zmelt, zdvo, zdvob, zdva, zdM, zSs, zdMe, zdMb, zdMv 
    5353      REAL(wp) ::   zMnew, zMnew1, zMnew2, zheat_hcflux, zheat_latent, z1_12 
    5454      REAL(wp) ::   zMbits, znMbits, zdMbitsE, zdMbitsM, zLbits, zAbits, zMbb 
    55       REAL(wp) ::   zxi, zyj, zff, z1_rday, z1_e1e2, zdt, z1_dt, z1_dt_e1e2 
     55      REAL(wp) ::   zxi, zyj, zff, z1_rday, z1_e1e2, zdt, z1_dt, z1_dt_e1e2, zdepw 
     56      REAL(wp), DIMENSION(jpk) :: ztoce, zuoce, zvoce, ze3t, ztmp 
    5657      TYPE(iceberg), POINTER ::   this, next 
    5758      TYPE(point)  , POINTER ::   pt 
     
    8384         pt => this%current_point 
    8485         nknberg = this%number(1) 
     86 
     87         CALL icb_utl_interp( pt%xi, pt%yj,           &   ! position 
     88             &                 pssu=pt%ssu, pua=pt%ua, &   ! oce/atm velocities 
     89             &                 pssv=pt%ssv, pva=pt%va, &   ! oce/atm velocities 
     90             &                 psst=pt%sst, pcn=pt%cn ) 
     91 
    8592         IF ( nn_sample_rate > 0 .AND. MOD(kt-1,nn_sample_rate) == 0 ) THEN 
    8693            CALL icb_utl_interp( pt%xi, pt%yj, pe1=pt%e1, pe2=pt%e2,                 & 
    87                &                 puo=pt%uo, pui=pt%ui, pua=pt%ua, pssh_i=pt%ssh_x,   & 
    88                &                 pvo=pt%vo, pvi=pt%vi, pva=pt%va, pssh_j=pt%ssh_y,   & 
    89                &                 psst=pt%sst, pcn=pt%cn, phi=pt%hi,                  & 
    90                &                 plat=pt%lat, plon=pt%lon                            ) 
    91          ELSE 
    92             CALL icb_utl_interp( pt%xi, pt%yj,           &   ! position 
    93                &                 puo=pt%uo, pua=pt%ua,   &   ! oce/atm velocities 
    94                &                 pvo=pt%vo, pva=pt%va,   &   ! oce/atm velocities 
    95                &                 psst=pt%sst, pcn=pt%cn )    ! sst and ice concentration 
    96                ! preparation Nacho add t3d and uo, vo, to basal 
     94               &                 pui=pt%ui, pssh_i=pt%ssh_x, & 
     95               &                 pvi=pt%vi, pssh_j=pt%ssh_y, & 
     96               &                 phi=pt%hi,                  & 
     97               &                 plat=pt%lat, plon=pt%lon ) 
    9798         END IF 
    9899         ! 
     
    101102         zM   = pt%mass 
    102103         zT   = pt%thickness                               ! total thickness 
    103        ! D   = (rn_rho_bergs/pp_rho_seawater)*zT ! draught (keel depth) 
     104         zD   = rho_berg_1_oce * zT          ! draught (keel depth) 
    104105       ! F   = zT - D ! freeboard 
    105106         zW   = pt%width 
     
    114115 
    115116         ! Environment 
    116          zdvo = SQRT( (pt%uvel-pt%uo)**2 + (pt%vvel-pt%vo)**2 ) 
    117          zdva = SQRT( (pt%ua  -pt%uo)**2 + (pt%va  -pt%vo)**2 ) 
    118          zSs  = 1.5_wp * SQRT( zdva ) + 0.1_wp * zdva                ! Sea state      (eqn M.A9) 
    119  
     117         ! default sst, ssu and ssv 
     118         ! ln_M2016: use temp, u and v profile 
     119         IF ( ln_M2016 ) THEN 
     120 
     121            ! load t, u, v and e3 profile at icb position 
     122            CALL icb_utl_interp( pt%xi, pt%yj, ptoce=ztoce, puoce=zuoce, pvoce=zvoce, pe3t=ze3t ) 
     123             
     124            !compute bottom level 
     125            CALL icb_utl_getkb( pt%kb, ze3t, zD ) 
     126 
     127            ikb = MIN(pt%kb,mbkt(ii,ij)) 
     128            ztb = ztoce(ikb)                                                ! basal temperature 
     129            zub = zuoce(ikb) 
     130            zvb = zvoce(ikb) 
     131         ELSE 
     132            ztb = pt%sst 
     133            zub = pt%ssu 
     134            zvb = pt%ssv 
     135         END IF 
     136 
     137         zdvob = SQRT( (pt%uvel-zub)**2 + (pt%vvel-zvb)**2 )        ! relative basal velocity 
     138         zdva  = SQRT( (pt%ua  -pt%ssu)**2 + (pt%va  -pt%ssv)**2 )  ! relative wind 
     139         zSs   = 1.5_wp * SQRT( zdva ) + 0.1_wp * zdva              ! Sea state      (eqn M.A9) 
     140         ! 
    120141         ! Melt rates in m/s (i.e. division by rday) 
    121          zMv = MAX( 7.62d-3*zSST+1.29d-3*(zSST**2)                    , 0._wp ) * z1_rday   ! Buoyant convection at sides (eqn M.A10) 
    122          zMb = MAX( 0.58_wp*(zdvo**0.8_wp)*(zSST+4.0_wp)/(zL**0.2_wp) , 0._wp ) * z1_rday   ! Basal turbulent melting     (eqn M.A7 ) 
    123          zMe = MAX( z1_12*(zSST+2.)*zSs*(1._wp+COS(rpi*(zIC**3)))     , 0._wp ) * z1_rday   ! Wave erosion                (eqn M.A8 ) 
     142         IF ( ln_M2016 ) THEN 
     143            ! Buoyant convection at sides (eqn M.A10) but averaging along all the iceberg draft 
     144            ztmp(:) = MAX( 7.62d-3*ztoce(:)+1.29d-3*(ztoce(:)**2), 0._wp ) * z1_rday 
     145            zMv = 0.0; zdepw = 0.0 
     146            DO jk = 1,ikb-1 
     147               zMv = zMv + ze3t(jk)*ztmp(jk) 
     148               zdepw = zdepw + ze3t(jk) 
     149            END DO 
     150            zMv = (zMv + (zD - zdepw)*ztmp(ikb)) / zD 
     151         ELSE 
     152            zMv = MAX( 7.62d-3*zSST+1.29d-3*(zSST**2), 0._wp ) * z1_rday 
     153         END IF 
     154 
     155         ! Basal turbulent melting     (eqn M.A7 ) but using the basal temperature 
     156         zMb = MAX( 0.58_wp*(zdvob**0.8_wp)*(ztb+4.0_wp)/(zL**0.2_wp) , 0._wp ) * z1_rday 
     157 
     158         ! Wave erosion                (eqn M.A8 ) 
     159         zMe = MAX( z1_12*(zSST+2.)*zSs*(1._wp+COS(rpi*(zIC**3)))     , 0._wp ) * z1_rday 
    124160 
    125161         IF( ln_operator_splitting ) THEN      ! Operator split update of volume/mass 
     
    209245 
    210246         ! Rolling 
    211          zDn = ( rn_rho_bergs / pp_rho_seawater ) * zTn       ! draught (keel depth) 
     247         zDn = rho_berg_1_oce * zTn       ! draught (keel depth) 
    212248         IF( zDn > 0._wp .AND. MAX(zWn,zLn) < SQRT( 0.92*(zDn**2) + 58.32*zDn ) ) THEN 
    213249            zT  = zTn 
  • NEMO/branches/2020/tickets_icb_1900/src/OCE/ICB/icbtrj.F90

    r13062 r13357  
    4040   INTEGER ::   numberid, nstepid, nscaling_id 
    4141   INTEGER ::   nlonid, nlatid, nxid, nyid, nuvelid, nvvelid, nmassid 
    42    INTEGER ::   nuoid, nvoid, nuaid, nvaid, nuiid, nviid 
     42   INTEGER ::   nssuid, nssvid, nuaid, nvaid, nuiid, nviid 
    4343   INTEGER ::   nsshxid, nsshyid, nsstid, ncntid, nthkid 
    4444   INTEGER ::   nthicknessid, nwidthid, nlengthid 
     
    111111      iret = NF90_DEF_VAR( ntrajid, 'uvel'          , NF90_DOUBLE, n_dim          , nuvelid          ) 
    112112      iret = NF90_DEF_VAR( ntrajid, 'vvel'          , NF90_DOUBLE, n_dim          , nvvelid          ) 
    113       iret = NF90_DEF_VAR( ntrajid, 'uto'           , NF90_DOUBLE, n_dim          , nuoid            ) 
    114       iret = NF90_DEF_VAR( ntrajid, 'vto'           , NF90_DOUBLE, n_dim          , nvoid            ) 
     113      iret = NF90_DEF_VAR( ntrajid, 'ssu'           , NF90_DOUBLE, n_dim          , nssuid           ) 
     114      iret = NF90_DEF_VAR( ntrajid, 'ssv'           , NF90_DOUBLE, n_dim          , nssvid           ) 
    115115      iret = NF90_DEF_VAR( ntrajid, 'uta'           , NF90_DOUBLE, n_dim          , nuaid            ) 
    116116      iret = NF90_DEF_VAR( ntrajid, 'vta'           , NF90_DOUBLE, n_dim          , nvaid            ) 
     
    148148      iret = NF90_PUT_ATT( ntrajid, nvvelid         , 'long_name', 'meridional velocity' ) 
    149149      iret = NF90_PUT_ATT( ntrajid, nvvelid         , 'units'    , 'm/s' ) 
    150       iret = NF90_PUT_ATT( ntrajid, nuoid           , 'long_name', 'ocean u component' ) 
    151       iret = NF90_PUT_ATT( ntrajid, nuoid           , 'units'    , 'm/s' ) 
    152       iret = NF90_PUT_ATT( ntrajid, nvoid           , 'long_name', 'ocean v component' ) 
    153       iret = NF90_PUT_ATT( ntrajid, nvoid           , 'units'    , 'm/s' ) 
     150      iret = NF90_PUT_ATT( ntrajid, nssuid          , 'long_name', 'ocean u component' ) 
     151      iret = NF90_PUT_ATT( ntrajid, nssuid          , 'units'    , 'm/s' ) 
     152      iret = NF90_PUT_ATT( ntrajid, nssvid          , 'long_name', 'ocean v component' ) 
     153      iret = NF90_PUT_ATT( ntrajid, nssvid          , 'units'    , 'm/s' ) 
    154154      iret = NF90_PUT_ATT( ntrajid, nuaid           , 'long_name', 'atmosphere u component' ) 
    155155      iret = NF90_PUT_ATT( ntrajid, nuaid           , 'units'    , 'm/s' ) 
     
    231231         iret = NF90_PUT_VAR( ntrajid, nuvelid         , pt%uvel          , (/ jn /) ) 
    232232         iret = NF90_PUT_VAR( ntrajid, nvvelid         , pt%vvel          , (/ jn /) ) 
    233          iret = NF90_PUT_VAR( ntrajid, nuoid           , pt%uo            , (/ jn /) ) 
    234          iret = NF90_PUT_VAR( ntrajid, nvoid           , pt%vo            , (/ jn /) ) 
     233         iret = NF90_PUT_VAR( ntrajid, nssuid          , pt%ssu           , (/ jn /) ) 
     234         iret = NF90_PUT_VAR( ntrajid, nssvid          , pt%ssv           , (/ jn /) ) 
    235235         iret = NF90_PUT_VAR( ntrajid, nuaid           , pt%ua            , (/ jn /) ) 
    236236         iret = NF90_PUT_VAR( ntrajid, nvaid           , pt%va            , (/ jn /) ) 
  • NEMO/branches/2020/tickets_icb_1900/src/OCE/ICB/icbutl.F90

    r13265 r13357  
    1919   !!---------------------------------------------------------------------- 
    2020   USE par_oce                             ! ocean parameters 
     21   USE oce,    ONLY: ts, uu, vv 
    2122   USE dom_oce                             ! ocean domain 
    2223   USE in_out_manager                      ! IO parameters 
     
    3435   PRIVATE 
    3536 
     37   INTERFACE icb_utl_bilin_h 
     38      MODULE PROCEDURE icb_utl_bilin_2d_h, icb_utl_bilin_3d_h 
     39   END INTERFACE 
     40 
    3641   PUBLIC   icb_utl_copy          ! routine called in icbstp module 
     42   PUBLIC   icb_utl_getkb         ! routine called in icbdyn and icbthm modules 
    3743   PUBLIC   icb_utl_interp        ! routine called in icbdyn, icbthm modules 
    3844   PUBLIC   icb_utl_bilin_h       ! routine called in icbdyn module 
     
    5662CONTAINS 
    5763 
    58    SUBROUTINE icb_utl_copy() 
     64   SUBROUTINE icb_utl_copy( Kmm ) 
    5965      !!---------------------------------------------------------------------- 
    6066      !!                  ***  ROUTINE icb_utl_copy  *** 
     
    6470      !! ** Method  : - blah blah 
    6571      !!---------------------------------------------------------------------- 
     72      REAL(wp), DIMENSION(0:jpi+1,0:jpj+1) :: ztmp 
    6673#if defined key_si3 
    6774      REAL(wp), DIMENSION(jpi,jpj) :: zssh_lead_m    !    ocean surface (ssh_m) if ice is not embedded 
    6875      !                                              !    ocean surface in leads if ice is embedded    
    6976#endif 
     77      INTEGER :: jk   ! vertical loop index 
     78      INTEGER :: Kmm  ! ocean time levelindex 
     79      ! 
    7080      ! copy nemo forcing arrays into iceberg versions with extra halo 
    7181      ! only necessary for variables not on T points 
    7282      ! and ssh which is used to calculate gradients 
    73  
    74       uo_e(1:jpi,1:jpj) = ssu_m(:,:) * umask(:,:,1) 
    75       vo_e(1:jpi,1:jpj) = ssv_m(:,:) * vmask(:,:,1) 
    76       tt_e(1:jpi,1:jpj) = sst_m(:,:) 
    77       fr_e(1:jpi,1:jpj) = fr_i (:,:) 
    78       ua_e(1:jpi,1:jpj) = utau (:,:) * umask(:,:,1) ! maybe mask useless because mask applied in sbcblk 
    79       va_e(1:jpi,1:jpj) = vtau (:,:) * vmask(:,:,1) ! maybe mask useless because mask applied in sbcblk 
    80       ! 
    81       CALL lbc_lnk_icb( 'icbutl', uo_e, 'U', -1._wp, 1, 1 ) 
    82       CALL lbc_lnk_icb( 'icbutl', vo_e, 'V', -1._wp, 1, 1 ) 
     83      ! 
     84      ! surface forcing 
     85      ! 
     86      ssu_e(1:jpi,1:jpj) = ssu_m(:,:) * umask(:,:,1) 
     87      ssv_e(1:jpi,1:jpj) = ssv_m(:,:) * vmask(:,:,1) 
     88      sst_e(1:jpi,1:jpj) = sst_m(:,:) 
     89      fr_e (1:jpi,1:jpj) = fr_i (:,:) 
     90      ua_e (1:jpi,1:jpj) = utau (:,:) * umask(:,:,1) ! maybe mask useless because mask applied in sbcblk 
     91      va_e (1:jpi,1:jpj) = vtau (:,:) * vmask(:,:,1) ! maybe mask useless because mask applied in sbcblk 
     92      ! 
     93      CALL lbc_lnk_icb( 'icbutl', ssu_e, 'U', -1._wp, 1, 1 ) 
     94      CALL lbc_lnk_icb( 'icbutl', ssv_e, 'V', -1._wp, 1, 1 ) 
    8395      CALL lbc_lnk_icb( 'icbutl', ua_e, 'U', -1._wp, 1, 1 ) 
    8496      CALL lbc_lnk_icb( 'icbutl', va_e, 'V', -1._wp, 1, 1 ) 
     
    98110#endif 
    99111      ! 
     112      ! (PM) could be improve with a 3d lbclnk gathering both variables 
     113      ! should be done once extra haloe generalised 
     114      IF ( ln_M2016 ) THEN 
     115         DO jk = 1,jpk 
     116            ! uoce 
     117            ztmp(1:jpi,1:jpj) = uu(:,:,jk,Kmm) 
     118            CALL lbc_lnk_icb( 'icbutl', ztmp, 'U', -1._wp, 1, 1 ) 
     119            uoce_e(:,:,jk) = ztmp(:,:) 
     120            ! 
     121            ! voce 
     122            ztmp(1:jpi,1:jpj) = vv(:,:,jk,Kmm) 
     123            CALL lbc_lnk_icb( 'icbutl', ztmp, 'V', -1._wp, 1, 1 ) 
     124            voce_e(:,:,jk) = ztmp(:,:) 
     125         END DO 
     126         toce_e(1:jpi,1:jpj,1:jpk) = ts(:,:,:,1,Kmm) 
     127         e3t_e (1:jpi,1:jpj,1:jpk) = e3t(:,:,:,Kmm) 
     128      END IF 
     129      ! 
    100130   END SUBROUTINE icb_utl_copy 
    101131 
    102132 
    103    SUBROUTINE icb_utl_interp( pi, pj, pe1, puo, pui, pua, pssh_i,   & 
    104       &                               pe2, pvo, pvi, pva, pssh_j,   & 
    105       &                               psst, pcn, phi, pff, plon, plat) 
     133   SUBROUTINE icb_utl_interp( pi, pj, pe1, pssu, pui, pua, pssh_i,     & 
     134      &                               pe2, pssv, pvi, pva, pssh_j,     & 
     135      &                               psst, pcn, phi, pff, plon, plat, & 
     136      &                               ptoce, puoce, pvoce, pe3t       ) 
    106137      !!---------------------------------------------------------------------- 
    107138      !!                  ***  ROUTINE icb_utl_interp  *** 
     
    122153      REAL(wp), INTENT(in   ) ::   pi , pj                        ! position in (i,j) referential 
    123154      REAL(wp), INTENT(  out), OPTIONAL ::   pe1, pe2                       ! i- and j scale factors 
    124       REAL(wp), INTENT(  out), OPTIONAL ::   puo, pvo, pui, pvi, pua, pva   ! ocean, ice and wind speeds 
     155      REAL(wp), INTENT(  out), OPTIONAL ::   pssu, pssv, pui, pvi, pua, pva   ! ocean, ice and wind speeds 
    125156      REAL(wp), INTENT(  out), OPTIONAL ::   pssh_i, pssh_j                 ! ssh i- & j-gradients 
    126157      REAL(wp), INTENT(  out), OPTIONAL ::   psst, pcn, phi, pff            ! SST, ice concentration, ice thickness, Coriolis 
    127158      REAL(wp), INTENT(  out), OPTIONAL ::   plat, plon                     ! position 
     159      REAL(wp), DIMENSION(jpk), INTENT(  out), OPTIONAL ::   ptoce, puoce, pvoce, pe3t   ! 3D variables 
    128160      ! 
    129161      REAL(wp), DIMENSION(4) :: zwT  , zwU  , zwV  , zwF   ! interpolation weight 
     
    147179      IF ( PRESENT(plat) ) plat= icb_utl_bilin_h( rlat_e, iiT, ijT, zwT, .false. ) 
    148180      ! 
    149       IF ( PRESENT(puo ) ) puo  = icb_utl_bilin_h( uo_e, iiU, ijU, zwU        , .false. )         ! ocean velocities 
    150       IF ( PRESENT(pvo ) ) pvo  = icb_utl_bilin_h( vo_e, iiV, ijV, zwV        , .false. )         ! 
    151       IF ( PRESENT(psst) ) psst = icb_utl_bilin_h( tt_e, iiT, ijT, zwT * zmskT, .false. ) ! sst 
    152       IF ( PRESENT(pcn ) ) pcn  = icb_utl_bilin_h( fr_e, iiT, ijT, zwT * zmskT, .false. ) ! ice concentration 
    153       IF ( PRESENT(pff ) ) pff  = icb_utl_bilin_h( ff_e, iiF, ijF, zwF        , .false. )         ! Coriolis parameter 
     181      IF ( PRESENT(pssu) ) pssu = icb_utl_bilin_h( ssu_e, iiU, ijU, zwU        , .false. )         ! ocean velocities 
     182      IF ( PRESENT(pssv) ) pssv = icb_utl_bilin_h( ssv_e, iiV, ijV, zwV        , .false. )         ! 
     183      IF ( PRESENT(psst) ) psst = icb_utl_bilin_h( sst_e, iiT, ijT, zwT * zmskT, .false. ) ! sst 
     184      IF ( PRESENT(pcn ) ) pcn  = icb_utl_bilin_h( fr_e , iiT, ijT, zwT * zmskT, .false. ) ! ice concentration 
     185      IF ( PRESENT(pff ) ) pff  = icb_utl_bilin_h( ff_e , iiF, ijF, zwF        , .false. )         ! Coriolis parameter 
    154186      ! 
    155187      IF ( PRESENT(pua) .AND. PRESENT(pva) ) THEN 
     
    188220            &       icb_utl_bilin_h( ssh_e, iiTm, ijTm, zwTm*zmskTm, .false. )  ) / ( 0.2_wp * pe2 ) 
    189221      END IF 
     222      ! 
     223      ! 3d interpolation 
     224      IF ( PRESENT(puoce) .AND. PRESENT(pvoce) ) THEN 
     225         puoce(:) = icb_utl_bilin_h( uoce_e , iiU, ijU, zwU ) 
     226         pvoce(:) = icb_utl_bilin_h( voce_e , iiV, ijV, zwV ) 
     227      END IF 
     228      IF ( PRESENT(ptoce) ) ptoce(:) = icb_utl_bilin_h( toce_e , iiT, ijT, zwT * zmskT ) 
     229      IF ( PRESENT(pe3t)  ) pe3t(:)  = e3t_e(iiT,ijT,:)    ! as in Nacho tarball need to be fix once we are able to reproduce Nacho results 
    190230      ! 
    191231   END SUBROUTINE icb_utl_interp 
     
    291331   END SUBROUTINE icb_utl_pos 
    292332 
    293    REAL(wp) FUNCTION icb_utl_bilin_h( pfld, pii, pij, pw, pllon ) 
     333   REAL(wp) FUNCTION icb_utl_bilin_2d_h( pfld, pii, pij, pw, pllon ) 
    294334      !!---------------------------------------------------------------------- 
    295335      !!                  ***  FUNCTION icb_utl_bilin  *** 
     
    321361      ! 
    322362      ! compute interpolated value 
    323       icb_utl_bilin_h = ( zdat(1)*pw(1) + zdat(2)*pw(2) + zdat(3)*pw(3) + zdat(4)*pw(4) ) / MAX(1.e-20, pw(1)+pw(2)+pw(3)+pw(4))  
    324       ! 
    325       IF( pllon .AND. icb_utl_bilin_h > 180._wp ) icb_utl_bilin_h = icb_utl_bilin_h - 360._wp 
    326       ! 
    327    END FUNCTION icb_utl_bilin_h 
     363      icb_utl_bilin_2d_h = ( zdat(1)*pw(1) + zdat(2)*pw(2) + zdat(3)*pw(3) + zdat(4)*pw(4) ) / MAX(1.e-20, pw(1)+pw(2)+pw(3)+pw(4))  
     364      ! 
     365      IF( pllon .AND. icb_utl_bilin_2d_h > 180._wp ) icb_utl_bilin_2d_h = icb_utl_bilin_2d_h - 360._wp 
     366      ! 
     367   END FUNCTION icb_utl_bilin_2d_h 
     368 
     369   FUNCTION icb_utl_bilin_3d_h( pfld, pii, pij, pw ) 
     370      !!---------------------------------------------------------------------- 
     371      !!                  ***  FUNCTION icb_utl_bilin  *** 
     372      !! 
     373      !! ** Purpose :   bilinear interpolation at berg location depending on the grid-point type 
     374      !!                this version deals with extra halo points 
     375      !! 
     376      !!       !!gm  CAUTION an optional argument should be added to handle 
     377      !!             the slip/no-slip conditions  ==>>> to be done later 
     378      !! 
     379      !!---------------------------------------------------------------------- 
     380      REAL(wp), DIMENSION(0:jpi+1,0:jpj+1, jpk), INTENT(in) ::   pfld      ! field to be interpolated 
     381      REAL(wp), DIMENSION(4)                   , INTENT(in) ::   pw        ! weight 
     382      INTEGER ,                                  INTENT(in) ::   pii, pij  ! bottom left corner 
     383      REAL(wp), DIMENSION(jpk) :: icb_utl_bilin_3d_h 
     384      ! 
     385      REAL(wp), DIMENSION(4,jpk) :: zdat ! input data 
     386      !!---------------------------------------------------------------------- 
     387      ! 
     388      ! data 
     389      zdat(1,:) = pfld(pii  ,pij  ,:) 
     390      zdat(2,:) = pfld(pii+1,pij  ,:) 
     391      zdat(3,:) = pfld(pii  ,pij+1,:) 
     392      zdat(4,:) = pfld(pii+1,pij+1,:) 
     393      ! 
     394      ! compute interpolated value 
     395      icb_utl_bilin_3d_h(:) = ( zdat(1,:)*pw(1) + zdat(2,:)*pw(2) + zdat(3,:)*pw(3) + zdat(4,:)*pw(4) ) / MAX(1.e-20, pw(1)+pw(2)+pw(3)+pw(4))  
     396      ! 
     397   END FUNCTION icb_utl_bilin_3d_h 
    328398 
    329399   REAL(wp) FUNCTION icb_utl_bilin_e( pet, peu, pev, pef, pi, pj ) 
     
    405475   END FUNCTION icb_utl_bilin_e 
    406476 
     477   SUBROUTINE icb_utl_getkb( kb, pe3, pD ) 
     478      !!---------------------------------------------------------------------- 
     479      !!                ***  ROUTINE icb_utl_getkb         *** 
     480      !! 
     481      !! ** Purpose :   compute the latest level affected by icb 
     482      !! 
     483      !!---------------------------------------------------------------------- 
     484      INTEGER, INTENT(out)               :: kb 
     485      REAL(wp), DIMENSION(:), INTENT(in) :: pe3 
     486      REAL(wp),               INTENT(in) :: pD 
     487      !! 
     488      INTEGER  :: jk 
     489      REAL(wp) :: zdepw 
     490      !! 
     491      zdepw = 0.0 
     492      kb = 1 
     493      DO WHILE ( zdepw <  pD) 
     494         zdepw = zdepw + pe3(kb) 
     495         kb = kb + 1 
     496      END DO 
     497      kb = kb - 1 
     498   END SUBROUTINE 
    407499 
    408500   SUBROUTINE icb_utl_add( bergvals, ptvals ) 
     
    593685      WRITE(numicb, 9200) kt, berg%number(1), & 
    594686                   pt%xi, pt%yj, pt%lon, pt%lat, pt%uvel, pt%vvel,  & 
    595                    pt%uo, pt%vo, pt%ua, pt%va, pt%ui, pt%vi 
     687                   pt%ssu, pt%ssv, pt%ua, pt%va, pt%ui, pt%vi 
    596688      CALL flush( numicb ) 
    597689 9200 FORMAT(5x,i5,2x,i10,6(2x,2f10.4)) 
     
    619711         WRITE(numicb,'(a," pe=(",i3,")")' ) cd_label, narea 
    620712         WRITE(numicb,'(a8,4x,a6,12x,a5,15x,a7,19x,a3,17x,a5,17x,a5,17x,a5)' )   & 
    621             &         'timestep', 'number', 'xi,yj','lon,lat','u,v','uo,vo','ua,va','ui,vi' 
     713            &         'timestep', 'number', 'xi,yj','lon,lat','u,v','ssu,ssv','ua,va','ui,vi' 
    622714      ENDIF 
    623715      DO WHILE( ASSOCIATED(this) ) 
  • NEMO/branches/2020/tickets_icb_1900/src/OCE/SBC/sbcmod.F90

    r13226 r13357  
    456456 
    457457      IF( ln_icebergs    )   THEN 
    458                                      CALL icb_stp( kt )           ! compute icebergs 
     458                                     CALL icb_stp( kt, Kmm )           ! compute icebergs 
    459459         ! Icebergs do not melt over the haloes.  
    460460         ! So emp values over the haloes are no more consistent with the inner domain values.  
Note: See TracChangeset for help on using the changeset viewer.