Changeset 3359
- Timestamp:
- 2012-04-18T12:42:56+02:00 (12 years ago)
- Location:
- branches/2012/dev_r3337_NOCS10_ICB/NEMOGCM/NEMO/OPA_SRC
- Files:
-
- 12 edited
-
DIA/diawri.F90 (modified) (2 diffs)
-
ICB/icb_oce.F90 (modified) (3 diffs)
-
ICB/icbclv.F90 (modified) (9 diffs)
-
ICB/icbdia.F90 (modified) (18 diffs)
-
ICB/icbdyn.F90 (modified) (13 diffs)
-
ICB/icbini.F90 (modified) (20 diffs)
-
ICB/icblbc.F90 (modified) (31 diffs)
-
ICB/icbrst.F90 (modified) (7 diffs)
-
ICB/icbrun.F90 (modified) (9 diffs)
-
ICB/icbthm.F90 (modified) (6 diffs)
-
ICB/icbtrj.F90 (modified) (6 diffs)
-
ICB/icbutl.F90 (modified) (12 diffs)
Legend:
- Unmodified
- Added
- Removed
-
branches/2012/dev_r3337_NOCS10_ICB/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90
r3340 r3359 442 442 CALL histdef( nid_T, "berg_melt" , "Melt rate of icebergs" , "kg/m2/s", & 443 443 & jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout ) 444 CALL histdef( nid_T, "berg_ melt_buoy" , "Buoyancy component of iceberg melt rate" , "kg/m2/s", &444 CALL histdef( nid_T, "berg_buoy_melt" , "Buoyancy component of iceberg melt rate" , "kg/m2/s", & 445 445 & jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout ) 446 CALL histdef( nid_T, "berg_ melt_eros" , "Erosion component of iceberg melt rate" , "kg/m2/s", &446 CALL histdef( nid_T, "berg_eros_melt" , "Erosion component of iceberg melt rate" , "kg/m2/s", & 447 447 & jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout ) 448 CALL histdef( nid_T, "berg_ melt_conv" , "Convective component of iceberg melt rate", "kg/m2/s", &448 CALL histdef( nid_T, "berg_conv_melt" , "Convective component of iceberg melt rate", "kg/m2/s", & 449 449 & jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout ) 450 450 CALL histdef( nid_T, "berg_virtual_area" , "Virtual coverage by icebergs" , "m2" , & … … 628 628 IF( ln_bergdia ) THEN 629 629 CALL histwrite( nid_T, "berg_melt" , it, berg_melt , ndim_hT, ndex_hT ) 630 CALL histwrite( nid_T, "berg_ melt_buoy" , it, melt_buoy, ndim_hT, ndex_hT )631 CALL histwrite( nid_T, "berg_ melt_eros" , it, melt_eros, ndim_hT, ndex_hT )632 CALL histwrite( nid_T, "berg_ melt_conv" , it, melt_conv, ndim_hT, ndex_hT )630 CALL histwrite( nid_T, "berg_buoy_melt" , it, buoy_melt , ndim_hT, ndex_hT ) 631 CALL histwrite( nid_T, "berg_eros_melt" , it, eros_melt , ndim_hT, ndex_hT ) 632 CALL histwrite( nid_T, "berg_conv_melt" , it, conv_melt , ndim_hT, ndex_hT ) 633 633 CALL histwrite( nid_T, "berg_virtual_area" , it, virtual_area , ndim_hT, ndex_hT ) 634 CALL histwrite( nid_T, "bits_src" , it, bits_src, ndim_hT, ndex_hT )635 CALL histwrite( nid_T, "bits_melt" , it, bits_melt, ndim_hT, ndex_hT )636 CALL histwrite( nid_T, "bits_mass" , it, bits_mass, ndim_hT, ndex_hT )634 CALL histwrite( nid_T, "bits_src" , it, bits_src , ndim_hT, ndex_hT ) 635 CALL histwrite( nid_T, "bits_melt" , it, bits_melt , ndim_hT, ndex_hT ) 636 CALL histwrite( nid_T, "bits_mass" , it, bits_mass , ndim_hT, ndex_hT ) 637 637 CALL histwrite( nid_T, "berg_mass" , it, berg_mass , ndim_hT, ndex_hT ) 638 638 ! -
branches/2012/dev_r3337_NOCS10_ICB/NEMOGCM/NEMO/OPA_SRC/ICB/icb_oce.F90
r3339 r3359 88 88 INTEGER :: current_year 89 89 REAL(wp) :: current_yearday ! 1.00-365.99 90 REAL(wp), DIMENSION(:), POINTER :: initial_width, initial_length91 LOGICAL :: restarted_bergs=.false. ! Indicate whether we read state from a restart or not90 REAL(wp), DIMENSION(:), POINTER :: first_width, first_length 91 LOGICAL :: l_restarted_bergs=.FALSE. ! Indicate whether we read state from a restart or not 92 92 ! ! arbitrary numbers for diawri entry 93 93 REAL(wp), DIMENSION(nclasses), PUBLIC :: class_num=(/ 1, 2, 3, 4, 5, 6, 7, 8, 9, 10 /) … … 107 107 108 108 !!gm almost all those PARAM ARE defined in NEMO 109 REAL(wp), PUBLIC, PARAMETER :: rho_ice = 916.7_wp ! Density of fresh ice @ 0oC (kg/m^3)110 REAL(wp), PUBLIC, PARAMETER :: rho_water = 999.8_wp ! Density of fresh water @ 0oC (kg/m^3)111 REAL(wp), PUBLIC, PARAMETER :: rho_air = 1.1_wp ! Density of air @ 0oC (kg/m^3) ???112 REAL(wp), PUBLIC, PARAMETER :: rho_seawater = 1025._wp ! Approx. density of surface sea water @ 0oC (kg/m^3)109 REAL(wp), PUBLIC, PARAMETER :: pp_rho_ice = 916.7_wp ! Density of fresh ice @ 0oC (kg/m^3) 110 REAL(wp), PUBLIC, PARAMETER :: pp_rho_water = 999.8_wp ! Density of fresh water @ 0oC (kg/m^3) 111 REAL(wp), PUBLIC, PARAMETER :: pp_rho_air = 1.1_wp ! Density of air @ 0oC (kg/m^3) ??? 112 REAL(wp), PUBLIC, PARAMETER :: pp_rho_seawater = 1025._wp ! Approx. density of surface sea water @ 0oC (kg/m^3) 113 113 !!gm end 114 REAL(wp), PUBLIC, PARAMETER :: Cd_av = 1.3_wp ! (Vertical) Drag coefficient between bergs and atmos (?)115 REAL(wp), PUBLIC, PARAMETER :: Cd_ah = 0.0055_wp ! (Horizontal) Drag coefficient between bergs and atmos (?)116 REAL(wp), PUBLIC, PARAMETER :: Cd_wv = 0.9_wp ! (Vertical) Drag coefficient between bergs and ocean (?)117 REAL(wp), PUBLIC, PARAMETER :: Cd_wh = 0.0012_wp ! (Horizontal) Drag coefficient between bergs and ocean (?)118 REAL(wp), PUBLIC, PARAMETER :: Cd_iv = 0.9_wp ! (Vertical) Drag coefficient between bergs and sea-ice (?)119 !TOM> no horizontal drag for sea ice! real, PARAMETER :: Cd_ih=0.0012 ! (Horizontal) Drag coefficient between bergs and sea-ice (?)114 REAL(wp), PUBLIC, PARAMETER :: pp_Cd_av = 1.3_wp ! (Vertical) Drag coefficient between bergs and atmos (?) 115 REAL(wp), PUBLIC, PARAMETER :: pp_Cd_ah = 0.0055_wp ! (Horizontal) Drag coefficient between bergs and atmos (?) 116 REAL(wp), PUBLIC, PARAMETER :: pp_Cd_wv = 0.9_wp ! (Vertical) Drag coefficient between bergs and ocean (?) 117 REAL(wp), PUBLIC, PARAMETER :: pp_Cd_wh = 0.0012_wp ! (Horizontal) Drag coefficient between bergs and ocean (?) 118 REAL(wp), PUBLIC, PARAMETER :: pp_Cd_iv = 0.9_wp ! (Vertical) Drag coefficient between bergs and sea-ice (?) 119 !TOM> no horizontal drag for sea ice! real, PARAMETER :: pp_Cd_ih=0.0012 ! (Horizontal) Drag coefficient between bergs and sea-ice (?) 120 120 121 121 ! !** namberg namelist parameters (and defaults) ** … … 147 147 148 148 ! Single instance of an icebergs type initialised in icebergs_init and updated in icebergs_run 149 REAL(wp), ALLOCATABLE, PUBLIC, SAVE, DIMENSION(:,:) :: p_calving, p_calving_hflx149 REAL(wp), ALLOCATABLE, PUBLIC, SAVE, DIMENSION(:,:) :: src_calving, src_calving_hflx 150 150 INTEGER , PUBLIC, SAVE :: numicb 151 INTEGER , PUBLIC, SAVE, DIMENSION(nkounts) :: kount_bergs152 INTEGER , PUBLIC, SAVE :: icbdi, icbei, icbdj,icbej153 REAL(wp) , PUBLIC, SAVE :: icb_left,icb_right154 INTEGER , PUBLIC, SAVE :: icbpack155 INTEGER , PUBLIC, SAVE :: ktberg,knberg156 INTEGER, ALLOCATABLE , PUBLIC, SAVE, DIMENSION(:) :: icbfldpts157 INTEGER, ALLOCATABLE , PUBLIC, SAVE, DIMENSION(:) :: icbflddest158 INTEGER, ALLOCATABLE , PUBLIC, SAVE, DIMENSION(:) :: icbfldproc151 INTEGER , PUBLIC, SAVE, DIMENSION(nkounts) :: num_bergs 152 INTEGER , PUBLIC, SAVE :: nicbdi, nicbei, nicbdj, nicbej 153 REAL(wp) , PUBLIC, SAVE :: ricb_left, ricb_right 154 INTEGER , PUBLIC, SAVE :: nicbpack 155 INTEGER , PUBLIC, SAVE :: nktberg, nknberg 156 INTEGER, ALLOCATABLE , PUBLIC, SAVE, DIMENSION(:) :: nicbfldpts 157 INTEGER, ALLOCATABLE , PUBLIC, SAVE, DIMENSION(:) :: nicbflddest 158 INTEGER, ALLOCATABLE , PUBLIC, SAVE, DIMENSION(:) :: nicbfldproc 159 159 160 160 TYPE(FLD), ALLOCATABLE, PUBLIC , DIMENSION(:) :: sf_icb ! structure: file information, fields read -
branches/2012/dev_r3337_NOCS10_ICB/NEMOGCM/NEMO/OPA_SRC/ICB/icbclv.F90
r3339 r3359 14 14 !! accumulate_calving : 15 15 !! icb_gen : generate test icebergs 16 !! icb_nam : read iceberg namelist17 16 !!---------------------------------------------------------------------- 18 17 USE par_oce ! NEMO parameters … … 29 28 PRIVATE 30 29 31 PUBLIC accumulate_calving ! routine called in xxx.F90 module32 PUBLIC calve_icebergs ! routine called in xxx.F90 module30 PUBLIC accumulate_calving ! routine called in icbrun.F90 module 31 PUBLIC calve_icebergs ! routine called in icbrun.F90 module 33 32 34 33 CONTAINS … … 43 42 !!---------------------------------------------------------------------- 44 43 ! 45 INTEGER :: kt 46 REAL(wp) :: calving_used, rdist, ufact 47 INTEGER :: jn, ji, jj, nmax 48 LOGICAL, SAVE :: first_call = .TRUE. 44 INTEGER, INTENT(in) :: kt 45 ! 46 REAL(wp) :: zcalving_used, zdist, zfact 47 INTEGER :: jn, ji, jj ! loop counters 48 INTEGER :: imx ! temporary integer for max berg class 49 LOGICAL, SAVE :: ll_first_call = .TRUE. 49 50 !!---------------------------------------------------------------------- 50 51 ! … … 52 53 ! Use interior mask: so no bergs in overlap areas and convert from km^3/year to kg/s 53 54 ! this assumes that input is given as equivalent water flux so that pure water density is appropriate 54 ufact = ( (1000._wp)**3 / ( NINT(rday) * nyear_len(1) ) ) * 1000._wp55 berg_grid%calving(:,:) = p_calving(:,:) * tmask_i(:,:) * ufact55 zfact = ( (1000._wp)**3 / ( NINT(rday) * nyear_len(1) ) ) * 1000._wp 56 berg_grid%calving(:,:) = src_calving(:,:) * tmask_i(:,:) * zfact 56 57 57 58 ! Heat in units of W/m2, and mask (just in case) 58 berg_grid%calving_hflx(:,:) = p_calving_hflx(:,:) * tmask_i(:,:)59 berg_grid%calving_hflx(:,:) = src_calving_hflx(:,:) * tmask_i(:,:) 59 60 60 IF( first_call .AND. .NOT.restarted_bergs) THEN ! This is a hack to simplify initialization61 first_call = .FALSE.61 IF( ll_first_call .AND. .NOT. l_restarted_bergs) THEN ! This is a hack to simplify initialization 62 ll_first_call = .FALSE. 62 63 !do jn=1, nclasses 63 64 ! where (berg_grid%calving==0.) berg_grid%stored_ice(:,:,jn)=0. … … 77 78 DO jj = 1,jpj 78 79 DO ji = 1,jpi 79 nmax = berg_grid%maxclass(ji,jj)80 rdist = SUM( rn_distribution(1:nclasses) ) / SUM( rn_distribution(1:nmax) )81 DO jn = 1, nmax80 imx = berg_grid%maxclass(ji,jj) 81 zdist = SUM( rn_distribution(1:nclasses) ) / SUM( rn_distribution(1:imx) ) 82 DO jn = 1, imx 82 83 berg_grid%stored_ice(ji,jj,jn) = berg_grid%stored_ice(ji,jj,jn) + & 83 berg_dt * berg_grid%calving(ji,jj) * rn_distribution(jn) * rdist84 berg_dt * berg_grid%calving(ji,jj) * rn_distribution(jn) * zdist 84 85 END DO 85 86 END DO … … 87 88 88 89 ! before changing the calving, save the amount we're about to use and do budget 89 calving_used = SUM( berg_grid%calving(:,:) )90 zcalving_used = SUM( berg_grid%calving(:,:) ) 90 91 berg_grid%tmp(:,:) = berg_dt * berg_grid%calving_hflx(:,:) * e1e2t(:,:) * tmask_i(:,:) 91 92 berg_grid%stored_heat (:,:) = berg_grid%stored_heat (:,:) + berg_grid%tmp(:,:) 92 CALL incoming_budget( kt, calving_used, berg_grid%tmp )93 CALL incoming_budget( kt, zcalving_used, berg_grid%tmp ) 93 94 ! 94 95 END SUBROUTINE accumulate_calving 95 96 ! ##############################################################################97 96 98 97 SUBROUTINE calve_icebergs() … … 114 113 TYPE(iceberg) :: newberg 115 114 TYPE(point) :: newpt 116 LOGICAL :: lret 117 INTEGER :: ji, jj, jn ! dummy loop indices 118 REAL(wp) :: xi, yj, ddt, calved_to_berg, heat_to_berg 115 INTEGER :: ji, jj, jn ! dummy loop indices 116 REAL(wp) :: zddt, zcalved_to_berg, zheat_to_berg 119 117 !!---------------------------------------------------------------------- 120 118 … … 122 120 123 121 DO jn = 1, nclasses 124 DO jj = icbdj,icbej125 DO ji = icbdi,icbei122 DO jj = nicbdj, nicbej 123 DO ji = nicbdi, nicbei 126 124 ! 127 ddt = 0._wp125 zddt = 0._wp 128 126 icnt = 0 129 127 ! … … 140 138 newpt%mass = rn_initial_mass (jn) 141 139 newpt%thickness = rn_initial_thickness(jn) 142 newpt%width = initial_width (jn)143 newpt%length = initial_length (jn)140 newpt%width = first_width (jn) 141 newpt%length = first_length (jn) 144 142 newberg%mass_scaling = rn_mass_scaling (jn) 145 143 newpt%mass_of_bits = 0._wp ! no bergy 146 144 ! 147 145 newpt%year = current_year 148 newpt%day = current_yearday + ddt / rday146 newpt%day = current_yearday + zddt / rday 149 147 newpt%heat_density = berg_grid%stored_heat(ji,jj) / berg_grid%stored_ice(ji,jj,jn) ! This is in J/kg 150 148 ! 151 149 CALL increment_kounter() 152 newberg%number(:) = kount_bergs(:)150 newberg%number(:) = num_bergs(:) 153 151 ! 154 152 CALL add_new_berg_to_list( newberg, newpt ) 155 153 ! 156 calved_to_berg = rn_initial_mass(jn) * rn_mass_scaling(jn) ! Units of kg154 zcalved_to_berg = rn_initial_mass(jn) * rn_mass_scaling(jn) ! Units of kg 157 155 ! ! Heat content 158 heat_to_berg =calved_to_berg * newpt%heat_density ! Units of J159 berg_grid%stored_heat(ji,jj) = berg_grid%stored_heat(ji,jj) - heat_to_berg156 zheat_to_berg = zcalved_to_berg * newpt%heat_density ! Units of J 157 berg_grid%stored_heat(ji,jj) = berg_grid%stored_heat(ji,jj) - zheat_to_berg 160 158 ! ! Stored mass 161 berg_grid%stored_ice(ji,jj,jn) = berg_grid%stored_ice(ji,jj,jn) - calved_to_berg159 berg_grid%stored_ice(ji,jj,jn) = berg_grid%stored_ice(ji,jj,jn) - zcalved_to_berg 162 160 ! 163 ddt =ddt + berg_dt * 2._wp / 17._wp ! Minor offset to start day !!gm why???161 zddt = zddt + berg_dt * 2._wp / 17._wp ! Minor offset to start day !!gm why??? 164 162 icnt = icnt + 1 165 163 ! 166 CALL calving_budget(ji, jj, jn, calved_to_berg,heat_to_berg )164 CALL calving_budget(ji, jj, jn, zcalved_to_berg, zheat_to_berg ) 167 165 END DO 168 166 icntmax = MAX( icntmax, icnt ) -
branches/2012/dev_r3337_NOCS10_ICB/NEMOGCM/NEMO/OPA_SRC/ICB/icbdia.F90
r3339 r3359 30 30 31 31 REAL(wp), DIMENSION(:,:) , POINTER, PUBLIC :: berg_melt=>NULL() ! Melting+erosion rate of icebergs (kg/s/m^2) 32 REAL(wp), DIMENSION(:,:) , POINTER, PUBLIC :: melt_buoy=>NULL() ! Buoyancy component of melting rate (kg/s/m^2)33 REAL(wp), DIMENSION(:,:) , POINTER, PUBLIC :: melt_eros=>NULL() ! Erosion component of melting rate (kg/s/m^2)34 REAL(wp), DIMENSION(:,:) , POINTER, PUBLIC :: melt_conv=>NULL() ! Convective component of melting rate (kg/s/m^2)35 REAL(wp), DIMENSION(:,:) , POINTER, PUBLIC :: bits_src=>NULL() ! Mass flux from berg erosion into bergy bits (kg/s/m^2)36 REAL(wp), DIMENSION(:,:) , POINTER, PUBLIC :: bits_melt=>NULL() ! Melting rate of bergy bits (kg/s/m^2)37 REAL(wp), DIMENSION(:,:) , POINTER, PUBLIC :: bits_mass=>NULL() ! Mass distribution of bergy bits (kg/s/m^2)32 REAL(wp), DIMENSION(:,:) , POINTER, PUBLIC :: buoy_melt=>NULL() ! Buoyancy component of melting rate (kg/s/m^2) 33 REAL(wp), DIMENSION(:,:) , POINTER, PUBLIC :: eros_melt=>NULL() ! Erosion component of melting rate (kg/s/m^2) 34 REAL(wp), DIMENSION(:,:) , POINTER, PUBLIC :: conv_melt=>NULL() ! Convective component of melting rate (kg/s/m^2) 35 REAL(wp), DIMENSION(:,:) , POINTER, PUBLIC :: bits_src=>NULL() ! Mass flux from berg erosion into bergy bits (kg/s/m^2) 36 REAL(wp), DIMENSION(:,:) , POINTER, PUBLIC :: bits_melt=>NULL() ! Melting rate of bergy bits (kg/s/m^2) 37 REAL(wp), DIMENSION(:,:) , POINTER, PUBLIC :: bits_mass=>NULL() ! Mass distribution of bergy bits (kg/s/m^2) 38 38 REAL(wp), DIMENSION(:,:) , POINTER, PUBLIC :: virtual_area=>NULL() ! Virtual surface coverage by icebergs (m^2) 39 REAL(wp), DIMENSION(:,:) , POINTER, PUBLIC :: berg_mass=>NULL() 39 REAL(wp), DIMENSION(:,:) , POINTER, PUBLIC :: berg_mass=>NULL() ! Mass distribution (kg/m^2) 40 40 REAL(wp), DIMENSION(:,:,:), POINTER, PUBLIC :: real_calving=>NULL() ! Calving rate into iceberg class at calving locations (kg/s) 41 41 REAL(wp), DIMENSION(:,:) , POINTER, PRIVATE :: tmpc=>NULL() ! Temporary work space 42 REAL(wp), DIMENSION(:) , POINTER, PRIVATE :: zsumbuf=>NULL() ! Temporary work space43 INTEGER , DIMENSION(:) , POINTER, PRIVATE :: isumbuf=>NULL() ! Temporary work space44 45 REAL(wp) , PRIVATE :: net_berg_melt46 REAL(wp) , PRIVATE :: net_bits_src47 REAL(wp) , PRIVATE :: net_bits_melt42 REAL(wp), DIMENSION(:) , POINTER, PRIVATE :: rsumbuf=>NULL() ! Temporary work space to reduce mpp exchanges 43 INTEGER , DIMENSION(:) , POINTER, PRIVATE :: nsumbuf=>NULL() ! Temporary work space to reduce mpp exchanges 44 45 REAL(wp) , PRIVATE :: berg_melt_net 46 REAL(wp) , PRIVATE :: bits_src_net 47 REAL(wp) , PRIVATE :: bits_melt_net 48 48 REAL(wp) , PRIVATE :: bits_mass_start, bits_mass_end 49 49 REAL(wp) , PRIVATE :: floating_heat_start, floating_heat_end 50 50 REAL(wp) , PRIVATE :: floating_mass_start, floating_mass_end 51 REAL(wp) , PRIVATE :: icebergs_mass_start, icebergs_mass_end51 REAL(wp) , PRIVATE :: bergs_mass_start, bergs_mass_end 52 52 REAL(wp) , PRIVATE :: stored_start, stored_heat_start 53 53 REAL(wp) , PRIVATE :: stored_end , stored_heat_end 54 REAL(wp) , PRIVATE :: net_incoming_calving, net_outgoing_calving55 REAL(wp) , PRIVATE :: net_incoming_calving_heat, net_outgoing_calving_heat56 REAL(wp) , PRIVATE :: net_incoming_calving_heat_used57 REAL(wp) , PRIVATE :: net_calving_received, net_calving_returned, net_calving_used58 REAL(wp) , PRIVATE :: net_heat_to_bergs, net_heat_to_ocean, net_melt59 REAL(wp) , PRIVATE :: net_calving_to_bergs54 REAL(wp) , PRIVATE :: calving_src_net, calving_out_net 55 REAL(wp) , PRIVATE :: calving_src_heat_net, calving_out_heat_net 56 REAL(wp) , PRIVATE :: calving_src_heat_used_net 57 REAL(wp) , PRIVATE :: calving_rcv_net, calving_ret_net, calving_used_net 58 REAL(wp) , PRIVATE :: heat_to_bergs_net, heat_to_ocean_net, melt_net 59 REAL(wp) , PRIVATE :: calving_to_bergs_net 60 60 61 61 INTEGER , PRIVATE :: nbergs_start, nbergs_end, nbergs_calved … … 64 64 INTEGER , DIMENSION(nclasses) , PRIVATE :: nbergs_calved_by_class 65 65 66 PUBLIC icb_budget_end ! routine called in xxx.F90 module67 PUBLIC icb_budget_init ! routine called in xxx.F90 module68 PUBLIC icb_budget ! routine called in xxx.F90 module69 PUBLIC icb_budget_step ! routine called in xxx.F90 module70 PUBLIC icb_budget_put ! routine called in xxx.F90 module71 PUBLIC melt_budget ! routine called in xxx.F90 module72 PUBLIC size_budget ! routine called in xxx.F90 module73 PUBLIC speed_budget ! routine called in xxx.F90 module74 PUBLIC calving_budget ! routine called in xxx.F90 module75 PUBLIC incoming_budget ! routine called in xxx.F90 module66 PUBLIC icb_budget_end ! routine called in icbrun.F90 module 67 PUBLIC icb_budget_init ! routine called in icbini.F90 module 68 PUBLIC icb_budget ! routine called in icbrun.F90 module 69 PUBLIC icb_budget_step ! routine called in icbrun.F90 module 70 PUBLIC icb_budget_put ! routine called in icbrun.F90 module 71 PUBLIC melt_budget ! routine called in icbthm.F90 module 72 PUBLIC size_budget ! routine called in icbthm.F90 module 73 PUBLIC speed_budget ! routine called in icbdyn.F90 module 74 PUBLIC calving_budget ! routine called in icbclv.F90 module 75 PUBLIC incoming_budget ! routine called in icbclv.F90 module 76 76 77 77 CONTAINS … … 83 83 IF( .NOT. ln_bergdia ) RETURN 84 84 DEALLOCATE( berg_melt ) 85 DEALLOCATE( melt_buoy)86 DEALLOCATE( melt_eros)87 DEALLOCATE( melt_conv)85 DEALLOCATE( buoy_melt ) 86 DEALLOCATE( eros_melt ) 87 DEALLOCATE( conv_melt ) 88 88 DEALLOCATE( bits_src ) 89 89 DEALLOCATE( bits_melt ) … … 94 94 DEALLOCATE( tmpc ) 95 95 IF( lk_mpp ) THEN 96 DEALLOCATE( zsumbuf )97 DEALLOCATE( isumbuf )96 DEALLOCATE( rsumbuf ) 97 DEALLOCATE( nsumbuf ) 98 98 ENDIF 99 99 … … 106 106 IF( .NOT. ln_bergdia ) RETURN 107 107 ALLOCATE( berg_melt (jpi,jpj) ) ; berg_melt (:,:) = 0._wp 108 ALLOCATE( melt_buoy (jpi,jpj) ) ; melt_buoy(:,:) = 0._wp109 ALLOCATE( melt_eros (jpi,jpj) ) ; melt_eros(:,:) = 0._wp110 ALLOCATE( melt_conv (jpi,jpj) ) ; melt_conv(:,:) = 0._wp111 ALLOCATE( bits_src (jpi,jpj) ) ; bits_src (:,:) = 0._wp112 ALLOCATE( bits_melt (jpi,jpj) ) ; bits_melt (:,:) = 0._wp113 ALLOCATE( bits_mass (jpi,jpj) ) ; bits_mass (:,:) = 0._wp108 ALLOCATE( buoy_melt (jpi,jpj) ) ; buoy_melt (:,:) = 0._wp 109 ALLOCATE( eros_melt (jpi,jpj) ) ; eros_melt (:,:) = 0._wp 110 ALLOCATE( conv_melt (jpi,jpj) ) ; conv_melt (:,:) = 0._wp 111 ALLOCATE( bits_src (jpi,jpj) ) ; bits_src (:,:) = 0._wp 112 ALLOCATE( bits_melt (jpi,jpj) ) ; bits_melt (:,:) = 0._wp 113 ALLOCATE( bits_mass (jpi,jpj) ) ; bits_mass (:,:) = 0._wp 114 114 ALLOCATE( virtual_area (jpi,jpj) ) ; virtual_area (:,:) = 0._wp 115 115 ALLOCATE( berg_mass (jpi,jpj) ) ; berg_mass (:,:) = 0._wp … … 117 117 ALLOCATE( tmpc(jpi,jpj) ) ; tmpc(:,:)=0. 118 118 119 nbergs_start 120 nbergs_end 121 stored_end 122 nbergs_start 123 stored_start 124 nbergs_melted 125 nbergs_calved 126 nbergs_calved_by_class(:) 127 nspeeding_tickets 128 stored_heat_end 129 floating_heat_end 130 floating_mass_end 131 icebergs_mass_end= 0._wp132 bits_mass_end 133 stored_heat_start 134 floating_heat_start 135 floating_mass_start 136 icebergs_mass_start= 0._wp137 bits_mass_start 138 bits_mass_end 139 net_calving_used= 0._wp140 net_calving_to_bergs= 0._wp141 net_heat_to_bergs= 0._wp142 net_heat_to_ocean= 0._wp143 net_calving_received= 0._wp144 net_calving_returned= 0._wp145 net_incoming_calving= 0._wp146 net_outgoing_calving= 0._wp147 net_incoming_calving_heat = 0._wp148 net_incoming_calving_heat_used= 0._wp149 net_outgoing_calving_heat = 0._wp150 net_melt= 0._wp151 net_berg_melt= 0._wp152 net_bits_melt= 0._wp153 net_bits_src= 0._wp154 155 floating_mass_start 156 icebergs_mass_start= sum_mass( first_berg, justbergs=.true. )157 bits_mass_start 119 nbergs_start = 0 120 nbergs_end = 0 121 stored_end = 0._wp 122 nbergs_start = 0._wp 123 stored_start = 0._wp 124 nbergs_melted = 0 125 nbergs_calved = 0 126 nbergs_calved_by_class(:) = 0 127 nspeeding_tickets = 0 128 stored_heat_end = 0._wp 129 floating_heat_end = 0._wp 130 floating_mass_end = 0._wp 131 bergs_mass_end = 0._wp 132 bits_mass_end = 0._wp 133 stored_heat_start = 0._wp 134 floating_heat_start = 0._wp 135 floating_mass_start = 0._wp 136 bergs_mass_start = 0._wp 137 bits_mass_start = 0._wp 138 bits_mass_end = 0._wp 139 calving_used_net = 0._wp 140 calving_to_bergs_net = 0._wp 141 heat_to_bergs_net = 0._wp 142 heat_to_ocean_net = 0._wp 143 calving_rcv_net = 0._wp 144 calving_ret_net = 0._wp 145 calving_src_net = 0._wp 146 calving_out_net = 0._wp 147 calving_src_heat_net = 0._wp 148 calving_src_heat_used_net = 0._wp 149 calving_out_heat_net = 0._wp 150 melt_net = 0._wp 151 berg_melt_net = 0._wp 152 bits_melt_net = 0._wp 153 bits_src_net = 0._wp 154 155 floating_mass_start = sum_mass( first_berg ) 156 bergs_mass_start = sum_mass( first_berg, justbergs=.true. ) 157 bits_mass_start = sum_mass( first_berg, justbits=.true. ) 158 158 IF( lk_mpp ) THEN 159 ALLOCATE( zsumbuf(23) ) ; zsumbuf(:) = 0._wp160 ALLOCATE( isumbuf(4+nclasses) ) ; isumbuf(:) = 0161 zsumbuf(1) = floating_mass_start162 zsumbuf(2) = icebergs_mass_start163 zsumbuf(3) = bits_mass_start164 CALL mpp_sum( zsumbuf(1:3), 3 )165 floating_mass_start = zsumbuf(1)166 icebergs_mass_start = zsumbuf(2)167 bits_mass_start = zsumbuf(3)159 ALLOCATE( rsumbuf(23) ) ; rsumbuf(:) = 0._wp 160 ALLOCATE( nsumbuf(4+nclasses) ) ; nsumbuf(:) = 0 161 rsumbuf(1) = floating_mass_start 162 rsumbuf(2) = bergs_mass_start 163 rsumbuf(3) = bits_mass_start 164 CALL mpp_sum( rsumbuf(1:3), 3 ) 165 floating_mass_start = rsumbuf(1) 166 bergs_mass_start = rsumbuf(2) 167 bits_mass_start = rsumbuf(3) 168 168 ENDIF 169 169 … … 172 172 !!------------------------------------------------------------------------- 173 173 174 SUBROUTINE icb_budget( l budge )174 SUBROUTINE icb_budget( ld_budge ) 175 175 ! Arguments 176 LOGICAL :: lbudge176 LOGICAL, INTENT(in) :: ld_budge 177 177 ! Local variables 178 INTEGER :: k179 REAL(wp) :: unused_calving, tmpsum, grdd_berg_mass,grdd_bits_mass180 181 IF( .NOT. ln_bergdia ) RETURN 182 183 unused_calving= SUM( berg_grid%calving(:,:) )184 tmpsum= SUM( berg_grid%floating_melt(:,:) * e1e2t(:,:) * tmask_i(:,:) )185 net_melt = net_melt +tmpsum * berg_dt186 net_outgoing_calving = net_outgoing_calving + ( unused_calving +tmpsum ) * berg_dt187 tmpsum= SUM( berg_melt(:,:) * e1e2t(:,:) * tmask_i(:,:) )188 net_berg_melt = net_berg_melt +tmpsum * berg_dt189 tmpsum= SUM( bits_src(:,:) * e1e2t(:,:) * tmask_i(:,:) )190 net_bits_src = net_bits_src +tmpsum * berg_dt191 tmpsum= SUM( bits_melt(:,:) * e1e2t(:,:) * tmask_i(:,:) )192 net_bits_melt = net_bits_melt +tmpsum * berg_dt193 tmpsum = SUM( p_calving(:,:) * tmask_i(:,:) )194 net_calving_returned = net_calving_returned+tmpsum * berg_dt195 tmpsum= SUM( berg_grid%calving_hflx(:,:) * e1e2t(:,:) * tmask_i(:,:) )196 net_outgoing_calving_heat = net_outgoing_calving_heat +tmpsum * berg_dt ! Units of J197 198 IF( l budge ) THEN178 INTEGER :: ik 179 REAL(wp) :: zunused_calving, ztmpsum, zgrdd_berg_mass, zgrdd_bits_mass 180 181 IF( .NOT. ln_bergdia ) RETURN 182 183 zunused_calving = SUM( berg_grid%calving(:,:) ) 184 ztmpsum = SUM( berg_grid%floating_melt(:,:) * e1e2t(:,:) * tmask_i(:,:) ) 185 melt_net = melt_net + ztmpsum * berg_dt 186 calving_out_net = calving_out_net + ( zunused_calving + ztmpsum ) * berg_dt 187 ztmpsum = SUM( berg_melt(:,:) * e1e2t(:,:) * tmask_i(:,:) ) 188 berg_melt_net = berg_melt_net + ztmpsum * berg_dt 189 ztmpsum = SUM( bits_src(:,:) * e1e2t(:,:) * tmask_i(:,:) ) 190 bits_src_net = bits_src_net + ztmpsum * berg_dt 191 ztmpsum = SUM( bits_melt(:,:) * e1e2t(:,:) * tmask_i(:,:) ) 192 bits_melt_net = bits_melt_net + ztmpsum * berg_dt 193 ztmpsum = SUM( src_calving(:,:) * tmask_i(:,:) ) 194 calving_ret_net = calving_ret_net + ztmpsum * berg_dt 195 ztmpsum = SUM( berg_grid%calving_hflx(:,:) * e1e2t(:,:) * tmask_i(:,:) ) 196 calving_out_heat_net = calving_out_heat_net + ztmpsum * berg_dt ! Units of J 197 198 IF( ld_budge ) THEN 199 199 stored_end = SUM( berg_grid%stored_ice(:,:,:) ) 200 200 stored_heat_end = SUM( berg_grid%stored_heat(:,:) ) 201 201 floating_mass_end = sum_mass( first_berg ) 202 icebergs_mass_end= sum_mass( first_berg,justbergs=.true. )203 bits_mass_end = sum_mass( first_berg,justbits=.true. )202 bergs_mass_end = sum_mass( first_berg,justbergs=.true. ) 203 bits_mass_end = sum_mass( first_berg,justbits=.true. ) 204 204 floating_heat_end = sum_heat( first_berg ) 205 205 206 206 nbergs_end = count_bergs() 207 grdd_berg_mass= SUM( berg_mass (:,:)*e1e2t(:,:)*tmask_i(:,:) )208 grdd_bits_mass = SUM( bits_mass(:,:)*e1e2t(:,:)*tmask_i(:,:) )207 zgrdd_berg_mass = SUM( berg_mass (:,:)*e1e2t(:,:)*tmask_i(:,:) ) 208 zgrdd_bits_mass = SUM( bits_mass(:,:)*e1e2t(:,:)*tmask_i(:,:) ) 209 209 210 210 IF( lk_mpp ) THEN 211 zsumbuf( 1) = stored_end212 zsumbuf( 2) = stored_heat_end213 zsumbuf( 3) = floating_mass_end214 zsumbuf( 4) = icebergs_mass_end215 zsumbuf( 5) = bits_mass_end216 zsumbuf( 6) = floating_heat_end217 zsumbuf( 7) = net_calving_returned218 zsumbuf( 8) = net_outgoing_calving219 zsumbuf( 9) = net_calving_received220 zsumbuf(10) = net_incoming_calving221 zsumbuf(11) = net_incoming_calving_heat222 zsumbuf(12) = net_incoming_calving_heat_used223 zsumbuf(13) = net_outgoing_calving_heat224 zsumbuf(14) = net_calving_used225 zsumbuf(15) = net_calving_to_bergs226 zsumbuf(16) = net_heat_to_bergs227 zsumbuf(17) = net_heat_to_ocean228 zsumbuf(18) = net_melt229 zsumbuf(19) = net_berg_melt230 zsumbuf(20) = net_bits_src231 zsumbuf(21) = net_bits_melt232 zsumbuf(22) =grdd_berg_mass233 zsumbuf(23) =grdd_bits_mass234 235 CALL mpp_sum( zsumbuf(1:23), 23)236 237 stored_end = zsumbuf( 1)238 stored_heat_end = zsumbuf( 2)239 floating_mass_end = zsumbuf( 3)240 icebergs_mass_end = zsumbuf( 4)241 bits_mass_end = zsumbuf( 5)242 floating_heat_end = zsumbuf( 6)243 net_calving_returned = zsumbuf( 7)244 net_outgoing_calving = zsumbuf( 8)245 net_calving_received = zsumbuf( 9)246 net_incoming_calving = zsumbuf(10)247 net_incoming_calving_heat = zsumbuf(11)248 net_incoming_calving_heat_used = zsumbuf(12)249 net_outgoing_calving_heat = zsumbuf(13)250 net_calving_used = zsumbuf(14)251 net_calving_to_bergs = zsumbuf(15)252 net_heat_to_bergs = zsumbuf(16)253 net_heat_to_ocean = zsumbuf(17)254 net_melt = zsumbuf(18)255 net_berg_melt = zsumbuf(19)256 net_bits_src = zsumbuf(20)257 net_bits_melt = zsumbuf(21)258 grdd_berg_mass = zsumbuf(22)259 grdd_bits_mass = zsumbuf(23)260 261 isumbuf(1) = nbergs_end262 isumbuf(2) = nbergs_calved263 isumbuf(3) = nbergs_melted264 isumbuf(4) = nspeeding_tickets265 DO k = 1,nclasses266 isumbuf(4+k) = nbergs_calved_by_class(k)211 rsumbuf( 1) = stored_end 212 rsumbuf( 2) = stored_heat_end 213 rsumbuf( 3) = floating_mass_end 214 rsumbuf( 4) = bergs_mass_end 215 rsumbuf( 5) = bits_mass_end 216 rsumbuf( 6) = floating_heat_end 217 rsumbuf( 7) = calving_ret_net 218 rsumbuf( 8) = calving_out_net 219 rsumbuf( 9) = calving_rcv_net 220 rsumbuf(10) = calving_src_net 221 rsumbuf(11) = calving_src_heat_net 222 rsumbuf(12) = calving_src_heat_used_net 223 rsumbuf(13) = calving_out_heat_net 224 rsumbuf(14) = calving_used_net 225 rsumbuf(15) = calving_to_bergs_net 226 rsumbuf(16) = heat_to_bergs_net 227 rsumbuf(17) = heat_to_ocean_net 228 rsumbuf(18) = melt_net 229 rsumbuf(19) = berg_melt_net 230 rsumbuf(20) = bits_src_net 231 rsumbuf(21) = bits_melt_net 232 rsumbuf(22) = zgrdd_berg_mass 233 rsumbuf(23) = zgrdd_bits_mass 234 235 CALL mpp_sum( rsumbuf(1:23), 23) 236 237 stored_end = rsumbuf( 1) 238 stored_heat_end = rsumbuf( 2) 239 floating_mass_end = rsumbuf( 3) 240 bergs_mass_end = rsumbuf( 4) 241 bits_mass_end = rsumbuf( 5) 242 floating_heat_end = rsumbuf( 6) 243 calving_ret_net = rsumbuf( 7) 244 calving_out_net = rsumbuf( 8) 245 calving_rcv_net = rsumbuf( 9) 246 calving_src_net = rsumbuf(10) 247 calving_src_heat_net = rsumbuf(11) 248 calving_src_heat_used_net = rsumbuf(12) 249 calving_out_heat_net = rsumbuf(13) 250 calving_used_net = rsumbuf(14) 251 calving_to_bergs_net = rsumbuf(15) 252 heat_to_bergs_net = rsumbuf(16) 253 heat_to_ocean_net = rsumbuf(17) 254 melt_net = rsumbuf(18) 255 berg_melt_net = rsumbuf(19) 256 bits_src_net = rsumbuf(20) 257 bits_melt_net = rsumbuf(21) 258 zgrdd_berg_mass = rsumbuf(22) 259 zgrdd_bits_mass = rsumbuf(23) 260 261 nsumbuf(1) = nbergs_end 262 nsumbuf(2) = nbergs_calved 263 nsumbuf(3) = nbergs_melted 264 nsumbuf(4) = nspeeding_tickets 265 DO ik = 1,nclasses 266 nsumbuf(4+ik) = nbergs_calved_by_class(ik) 267 267 ENDDO 268 268 269 CALL mpp_sum( isumbuf(1:nclasses+4), nclasses+4 )270 271 nbergs_end = isumbuf(1)272 nbergs_calved = isumbuf(2)273 nbergs_melted = isumbuf(3)274 nspeeding_tickets = isumbuf(4)275 DO k = 1,nclasses276 nbergs_calved_by_class( k)= isumbuf(4+k)269 CALL mpp_sum( nsumbuf(1:nclasses+4), nclasses+4 ) 270 271 nbergs_end = nsumbuf(1) 272 nbergs_calved = nsumbuf(2) 273 nbergs_melted = nsumbuf(3) 274 nspeeding_tickets = nsumbuf(4) 275 DO ik = 1,nclasses 276 nbergs_calved_by_class(ik)= nsumbuf(4+ik) 277 277 ENDDO 278 278 279 279 ENDIF 280 280 281 CALL report_state( 'stored ice','kg','',stored_start,'',stored_end,'')282 CALL report_state( 'floating','kg','',floating_mass_start,'',floating_mass_end,'',nbergs_end)283 CALL report_state( 'icebergs','kg','',icebergs_mass_start,'',icebergs_mass_end,'')284 CALL report_state( 'bits','kg','',bits_mass_start,'',bits_mass_end,'')285 CALL report_istate( 'berg #','',nbergs_start,'',nbergs_end,'')286 CALL report_ibudget( 'berg #','calved',nbergs_calved, &287 'melted',nbergs_melted, &288 '#',nbergs_start,nbergs_end)289 CALL report_budget( 'stored mass','kg','calving used',net_calving_used, &290 'bergs',net_calving_to_bergs, &291 292 CALL report_budget( 'floating mass','kg','calving used',net_calving_to_bergs, &293 'bergs',net_melt, &294 295 CALL report_budget( 'berg mass','kg','calving',net_calving_to_bergs, &296 'melt+eros',net_berg_melt, &297 'berg mass',icebergs_mass_start,icebergs_mass_end)298 CALL report_budget( 'bits mass','kg','eros used',net_bits_src, &299 'bergs',net_bits_melt, &300 'stored mass',bits_mass_start,bits_mass_end)301 CALL report_budget( 'net mass','kg','recvd',net_calving_received, &302 'rtrnd',net_calving_returned, &303 'net mass',stored_start+floating_mass_start, &304 stored_end+floating_mass_end)305 CALL report_consistant( 'iceberg mass','kg','gridded',grdd_berg_mass,'bergs',icebergs_mass_end)306 CALL report_consistant( 'bits mass','kg','gridded',grdd_bits_mass,'bits',bits_mass_end)307 CALL report_state( 'net heat','J','',stored_heat_start+floating_heat_start,'', &308 stored_heat_end+floating_heat_end,'')309 CALL report_state( 'stored heat','J','',stored_heat_start,'',stored_heat_end,'')310 CALL report_state( 'floating heat','J','',floating_heat_start,'',floating_heat_end,'')311 CALL report_budget( 'net heat','J','net heat',net_incoming_calving_heat, &312 'net heat',net_outgoing_calving_heat, &313 'net heat',stored_heat_start+floating_heat_start, &314 stored_heat_end+floating_heat_end)315 CALL report_budget( 'stored heat','J','calving used',net_incoming_calving_heat_used, &316 'bergs',net_heat_to_bergs, &317 'net heat',stored_heat_start,stored_heat_end)318 CALL report_budget( 'flting heat','J','calved',net_heat_to_bergs, &319 'melt',net_heat_to_ocean, &320 'net heat',floating_heat_start,floating_heat_end)281 CALL report_state( 'stored ice','kg','',stored_start,'',stored_end,'') 282 CALL report_state( 'floating','kg','',floating_mass_start,'',floating_mass_end,'',nbergs_end) 283 CALL report_state( 'icebergs','kg','',bergs_mass_start,'',bergs_mass_end,'') 284 CALL report_state( 'bits','kg','',bits_mass_start,'',bits_mass_end,'') 285 CALL report_istate( 'berg #','',nbergs_start,'',nbergs_end,'') 286 CALL report_ibudget( 'berg #','calved',nbergs_calved, & 287 'melted',nbergs_melted, & 288 '#',nbergs_start,nbergs_end) 289 CALL report_budget( 'stored mass','kg','calving used',calving_used_net, & 290 'bergs',calving_to_bergs_net, & 291 'stored mass',stored_start,stored_end) 292 CALL report_budget( 'floating mass','kg','calving used',calving_to_bergs_net, & 293 'bergs',melt_net, & 294 'stored mass',floating_mass_start,floating_mass_end) 295 CALL report_budget( 'berg mass','kg','calving',calving_to_bergs_net, & 296 'melt+eros',berg_melt_net, & 297 'berg mass',bergs_mass_start,bergs_mass_end) 298 CALL report_budget( 'bits mass','kg','eros used',bits_src_net, & 299 'bergs',bits_melt_net, & 300 'stored mass',bits_mass_start,bits_mass_end) 301 CALL report_budget( 'net mass','kg','recvd',calving_rcv_net, & 302 'rtrnd',calving_ret_net, & 303 'net mass',stored_start+floating_mass_start, & 304 stored_end+floating_mass_end) 305 CALL report_consistant( 'iceberg mass','kg','gridded',zgrdd_berg_mass,'bergs',bergs_mass_end) 306 CALL report_consistant( 'bits mass','kg','gridded',zgrdd_bits_mass,'bits',bits_mass_end) 307 CALL report_state( 'net heat','J','',stored_heat_start+floating_heat_start,'', & 308 stored_heat_end+floating_heat_end,'') 309 CALL report_state( 'stored heat','J','',stored_heat_start,'',stored_heat_end,'') 310 CALL report_state( 'floating heat','J','',floating_heat_start,'',floating_heat_end,'') 311 CALL report_budget( 'net heat','J','net heat',calving_src_heat_net, & 312 'net heat',calving_out_heat_net, & 313 'net heat',stored_heat_start+floating_heat_start, & 314 stored_heat_end+floating_heat_end) 315 CALL report_budget( 'stored heat','J','calving used',calving_src_heat_used_net, & 316 'bergs',heat_to_bergs_net, & 317 'net heat',stored_heat_start,stored_heat_end) 318 CALL report_budget( 'flting heat','J','calved',heat_to_bergs_net, & 319 'melt',heat_to_ocean_net, & 320 'net heat',floating_heat_start,floating_heat_end) 321 321 IF (nn_verbose_level >= 1) THEN 322 CALL report_consistant( 'top interface','kg','from SIS',net_incoming_calving, &323 'received',net_calving_received)324 CALL report_consistant( 'bot interface','kg','sent',net_outgoing_calving, &325 'returned',net_calving_returned)322 CALL report_consistant( 'top interface','kg','from SIS',calving_src_net, & 323 'received',calving_rcv_net) 324 CALL report_consistant( 'bot interface','kg','sent',calving_out_net, & 325 'returned',calving_ret_net) 326 326 ENDIF 327 WRITE( numicb,'("calved by class = ",i6,20(",",i6))') (nbergs_calved_by_class(k),k=1,nclasses)328 IF ( nspeeding_tickets .GT. 0 ) WRITE(numicb,'("speeding tickets issued = ",i6)') nspeeding_tickets329 330 nbergs_start 331 stored_start 332 nbergs_melted 333 nbergs_calved 334 nbergs_calved_by_class(:) 335 nspeeding_tickets 336 stored_heat_start 337 floating_heat_start 338 floating_mass_start 339 icebergs_mass_start = icebergs_mass_end340 bits_mass_start 341 net_calving_used= 0._wp342 net_calving_to_bergs= 0._wp343 net_heat_to_bergs= 0._wp344 net_heat_to_ocean= 0._wp345 net_calving_received= 0._wp346 net_calving_returned= 0._wp347 net_incoming_calving= 0._wp348 net_outgoing_calving= 0._wp349 net_incoming_calving_heat = 0._wp350 net_incoming_calving_heat_used= 0._wp351 net_outgoing_calving_heat = 0._wp352 net_melt= 0._wp353 net_berg_melt= 0._wp354 net_bits_melt= 0._wp355 net_bits_src= 0._wp327 WRITE( numicb, '("calved by class = ",i6,20(",",i6))') (nbergs_calved_by_class(ik),ik=1,nclasses) 328 IF ( nspeeding_tickets > 0 ) WRITE( numicb, '("speeding tickets issued = ",i6)') nspeeding_tickets 329 330 nbergs_start = nbergs_end 331 stored_start = stored_end 332 nbergs_melted = 0 333 nbergs_calved = 0 334 nbergs_calved_by_class(:) = 0 335 nspeeding_tickets = 0 336 stored_heat_start = stored_heat_end 337 floating_heat_start = floating_heat_end 338 floating_mass_start = floating_mass_end 339 bergs_mass_start = bergs_mass_end 340 bits_mass_start = bits_mass_end 341 calving_used_net = 0._wp 342 calving_to_bergs_net = 0._wp 343 heat_to_bergs_net = 0._wp 344 heat_to_ocean_net = 0._wp 345 calving_rcv_net = 0._wp 346 calving_ret_net = 0._wp 347 calving_src_net = 0._wp 348 calving_out_net = 0._wp 349 calving_src_heat_net = 0._wp 350 calving_src_heat_used_net = 0._wp 351 calving_out_heat_net = 0._wp 352 melt_net = 0._wp 353 berg_melt_net = 0._wp 354 bits_melt_net = 0._wp 355 bits_src_net = 0._wp 356 356 ENDIF 357 357 … … 366 366 IF( .NOT. ln_bergdia ) RETURN 367 367 berg_melt (:,:) = 0._wp 368 melt_buoy(:,:) = 0._wp369 melt_eros(:,:) = 0._wp370 melt_conv(:,:) = 0._wp371 bits_src (:,:) = 0._wp372 bits_melt (:,:) = 0._wp373 bits_mass (:,:) = 0._wp368 buoy_melt (:,:) = 0._wp 369 eros_melt (:,:) = 0._wp 370 conv_melt (:,:) = 0._wp 371 bits_src (:,:) = 0._wp 372 bits_melt (:,:) = 0._wp 373 bits_mass (:,:) = 0._wp 374 374 berg_mass (:,:) = 0._wp 375 375 virtual_area (:,:) = 0._wp … … 384 384 IF( .NOT. ln_bergdia ) RETURN 385 385 CALL iom_put( "berg_melt" , berg_melt (:,:) ) ! 'Melt rate of icebergs' , 'kg/m2/s' 386 CALL iom_put( "berg_ melt_buoy" , melt_buoy(:,:) ) ! 'Buoyancy component of iceberg melt rate' , 'kg/m2/s'387 CALL iom_put( "berg_ melt_eros" , melt_eros(:,:) ) ! 'Erosion component of iceberg melt rate' , 'kg/m2/s'388 CALL iom_put( "berg_ melt_conv" , melt_conv(:,:) ) ! 'Convective component of iceberg melt rate', 'kg/m2/s'386 CALL iom_put( "berg_buoy_melt" , buoy_melt (:,:) ) ! 'Buoyancy component of iceberg melt rate' , 'kg/m2/s' 387 CALL iom_put( "berg_eros_melt" , eros_melt (:,:) ) ! 'Erosion component of iceberg melt rate' , 'kg/m2/s' 388 CALL iom_put( "berg_conv_melt" , conv_melt (:,:) ) ! 'Convective component of iceberg melt rate', 'kg/m2/s' 389 389 CALL iom_put( "berg_virtual_area" , virtual_area (:,:) ) ! 'Virtual coverage by icebergs' , 'm2' 390 CALL iom_put( "bits_src" , bits_src (:,:) ) ! 'Mass source of bergy bits' , 'kg/m2/s'391 CALL iom_put( "bits_melt" , bits_melt (:,:) ) ! 'Melt rate of bergy bits' , 'kg/m2/s'392 CALL iom_put( "bits_mass" , bits_mass (:,:) ) ! 'Bergy bit density field' , 'kg/m2'390 CALL iom_put( "bits_src" , bits_src (:,:) ) ! 'Mass source of bergy bits' , 'kg/m2/s' 391 CALL iom_put( "bits_melt" , bits_melt (:,:) ) ! 'Melt rate of bergy bits' , 'kg/m2/s' 392 CALL iom_put( "bits_mass" , bits_mass (:,:) ) ! 'Bergy bit density field' , 'kg/m2' 393 393 CALL iom_put( "berg_mass" , berg_mass (:,:) ) ! 'Iceberg density field' , 'kg/m2' 394 394 CALL iom_put( "berg_real_calving" , real_calving (:,:,:) ) ! 'Calving into iceberg class' , 'kg/s' … … 398 398 !!------------------------------------------------------------------------- 399 399 400 SUBROUTINE calving_budget( ji, jj, jn, calved,heated )401 INTEGER :: ji, jj, jn402 REAL(wp) ::calved403 REAL(wp) ::heated404 405 IF( .NOT. ln_bergdia ) RETURN 406 real_calving( ji,jj,jn) = real_calving(ji,jj,jn) +calved / berg_dt400 SUBROUTINE calving_budget( ki, kj, kn, pcalved, pheated ) 401 INTEGER, INTENT(in) :: ki, kj, kn 402 REAL(wp), INTENT(in) :: pcalved 403 REAL(wp), INTENT(in) :: pheated 404 405 IF( .NOT. ln_bergdia ) RETURN 406 real_calving(ki,kj,kn) = real_calving(ki,kj,kn) + pcalved / berg_dt 407 407 nbergs_calved = nbergs_calved + 1 408 nbergs_calved_by_class( jn) = nbergs_calved_by_class(jn) + 1409 net_calving_to_bergs = net_calving_to_bergs +calved410 net_heat_to_bergs = net_heat_to_bergs +heated408 nbergs_calved_by_class(kn) = nbergs_calved_by_class(kn) + 1 409 calving_to_bergs_net = calving_to_bergs_net + pcalved 410 heat_to_bergs_net = heat_to_bergs_net + pheated 411 411 412 412 END SUBROUTINE calving_budget … … 414 414 !!------------------------------------------------------------------------- 415 415 416 SUBROUTINE incoming_budget( kt, calving_used,heat_used )416 SUBROUTINE incoming_budget( kt, pcalving_used, pheat_used ) 417 417 INTEGER , INTENT(in) :: kt 418 REAL(wp), INTENT(in) :: calving_used419 REAL(wp), DIMENSION(:,:), INTENT(in) :: heat_used418 REAL(wp), INTENT(in) :: pcalving_used 419 REAL(wp), DIMENSION(:,:), INTENT(in) :: pheat_used 420 420 421 421 IF( .NOT. ln_bergdia ) RETURN … … 431 431 ENDIF 432 432 433 net_calving_received = net_calving_received+ SUM( berg_grid%calving(:,:) ) * berg_dt434 net_incoming_calving = net_calving_received435 net_incoming_calving_heat = net_incoming_calving_heat + &433 calving_rcv_net = calving_rcv_net + SUM( berg_grid%calving(:,:) ) * berg_dt 434 calving_src_net = calving_rcv_net 435 calving_src_heat_net = calving_src_heat_net + & 436 436 SUM( berg_grid%calving_hflx(:,:) * e1e2t(:,:) ) * berg_dt ! Units of J 437 net_calving_used = net_calving_used +calving_used * berg_dt438 net_incoming_calving_heat_used = net_incoming_calving_heat_used + SUM(heat_used(:,:) )437 calving_used_net = calving_used_net + pcalving_used * berg_dt 438 calving_src_heat_used_net = calving_src_heat_used_net + SUM( pheat_used(:,:) ) 439 439 440 440 END SUBROUTINE incoming_budget … … 442 442 !!------------------------------------------------------------------------- 443 443 444 SUBROUTINE size_budget(ji, jj, Wn, Ln, Abits, mass_scale, Mnew, nMbits, z1_e1e2) 445 INTEGER :: ji, jj 446 REAL(wp) :: Wn, Ln, Abits, mass_scale, Mnew, nMbits, z1_e1e2 447 448 IF( .NOT. ln_bergdia ) RETURN 449 virtual_area(ji,jj) = virtual_area(ji,jj)+(Wn*Ln+Abits)*mass_scale ! m^2 450 berg_mass(ji,jj) = berg_mass(ji,jj) + Mnew * z1_e1e2 ! kg/m2 451 bits_mass(ji,jj) = bits_mass(ji,jj) + nMbits * z1_e1e2 ! kg/m2 444 SUBROUTINE size_budget(ki, kj, pWn, pLn, pAbits, & 445 & pmass_scale, pMnew, pnMbits, pz1_e1e2) 446 INTEGER, INTENT(in) :: ki, kj 447 REAL(wp), INTENT(in) :: pWn, pLn, pAbits, pmass_scale, pMnew, pnMbits, pz1_e1e2 448 449 IF( .NOT. ln_bergdia ) RETURN 450 virtual_area(ki,kj) = virtual_area(ki,kj) + ( pWn * pLn + pAbits ) * pmass_scale ! m^2 451 berg_mass(ki,kj) = berg_mass(ki,kj) + pMnew * pz1_e1e2 ! kg/m2 452 bits_mass(ki,kj) = bits_mass(ki,kj) + pnMbits * pz1_e1e2 ! kg/m2 452 453 453 454 END SUBROUTINE size_budget … … 464 465 !!------------------------------------------------------------------------- 465 466 466 SUBROUTINE melt_budget(ji, jj, mnew, heat, mass_scale, dM, dMbitsE, dMbitsM, dMb, dMe, dMv, z1_dt_e1e2) 467 468 INTEGER :: ji, jj 469 REAL(wp), INTENT(in) :: mnew, heat, mass_scale 470 REAL(wp), INTENT(in) :: dM, dMbitsE, dMbitsM, dMb, dMe, dMv, z1_dt_e1e2 471 472 IF( .NOT. ln_bergdia ) RETURN 473 474 berg_melt (ji,jj) = berg_melt (ji,jj) + dM * z1_dt_e1e2 ! kg/m2/s 475 bits_src (ji,jj) = bits_src (ji,jj) + dMbitsE * z1_dt_e1e2 ! mass flux into bergy bitskg/m2/s 476 bits_melt (ji,jj) = bits_melt (ji,jj) + dMbitsM * z1_dt_e1e2 ! melt rate of bergy bits kg/m2/s 477 melt_buoy (ji,jj) = melt_buoy (ji,jj) + dMb * z1_dt_e1e2 ! kg/m2/s 478 melt_eros (ji,jj) = melt_eros (ji,jj) + dMe * z1_dt_e1e2 ! erosion rate kg/m2/s 479 melt_conv (ji,jj) = melt_conv (ji,jj) + dMv * z1_dt_e1e2 ! kg/m2/s 480 net_heat_to_ocean = net_heat_to_ocean + heat * mass_scale * berg_dt ! J 467 SUBROUTINE melt_budget(ki, kj, pmnew, pheat, pmass_scale, & 468 & pdM, pdMbitsE, pdMbitsM, pdMb, pdMe, & 469 & pdMv, pz1_dt_e1e2) 470 471 INTEGER :: ki, kj 472 REAL(wp), INTENT(in) :: pmnew, pheat, pmass_scale 473 REAL(wp), INTENT(in) :: pdM, pdMbitsE, pdMbitsM, pdMb, pdMe, pdMv, pz1_dt_e1e2 474 475 IF( .NOT. ln_bergdia ) RETURN 476 477 berg_melt (ki,kj) = berg_melt (ki,kj) + pdM * pz1_dt_e1e2 ! kg/m2/s 478 bits_src (ki,kj) = bits_src (ki,kj) + pdMbitsE * pz1_dt_e1e2 ! mass flux into bergy bitskg/m2/s 479 bits_melt (ki,kj) = bits_melt (ki,kj) + pdMbitsM * pz1_dt_e1e2 ! melt rate of bergy bits kg/m2/s 480 buoy_melt (ki,kj) = buoy_melt (ki,kj) + pdMb * pz1_dt_e1e2 ! kg/m2/s 481 eros_melt (ki,kj) = eros_melt (ki,kj) + pdMe * pz1_dt_e1e2 ! erosion rate kg/m2/s 482 conv_melt (ki,kj) = conv_melt (ki,kj) + pdMv * pz1_dt_e1e2 ! kg/m2/s 483 heat_to_ocean_net = heat_to_ocean_net + pheat * pmass_scale * berg_dt ! J 481 484 IF( mnew <= 0._wp ) nbergs_melted = nbergs_melted + 1 ! Delete the berg if completely melted 482 485 483 484 486 END SUBROUTINE melt_budget 485 487 486 488 !!------------------------------------------------------------------------- 487 489 488 SUBROUTINE report_state(budgetstr,budgetunits,startstr,startval,endstr,endval,delstr,nbergs) 490 SUBROUTINE report_state( cd_budgetstr, cd_budgetunits, cd_startstr, pstartval, cd_endstr, & 491 & pendval, cd_delstr, kbergs ) 489 492 ! Arguments 490 CHARACTER*(*), INTENT(in) :: budgetstr, budgetunits, startstr, endstr,delstr491 REAL(wp), INTENT(in) :: startval,endval492 INTEGER, INTENT(in), OPTIONAL :: nbergs493 494 if (present(nbergs)) then495 WRITE(numicb,100) budgetstr//' state:',&496 startstr//' start',startval,budgetunits, &497 endstr//' end',endval,budgetunits, &498 'Delta '//delstr,endval-startval,budgetunits, &499 '# of bergs',nbergs500 else501 WRITE(numicb,100) budgetstr//' state:',&502 startstr//' start',startval,budgetunits, &503 endstr//' end',endval,budgetunits, &504 delstr//'Delta',endval-startval,budgetunits505 endif493 CHARACTER*(*), INTENT(in) :: cd_budgetstr, cd_budgetunits, cd_startstr, cd_endstr, cd_delstr 494 REAL(wp), INTENT(in) :: pstartval, pendval 495 INTEGER, INTENT(in), OPTIONAL :: kbergs 496 497 IF ( PRESENT(kbergs) ) THEN 498 WRITE(numicb,100) cd_budgetstr // ' state:', & 499 cd_startstr // ' start', pstartval, cd_budgetunits, & 500 cd_endstr // ' end', pendval, cd_budgetunits, & 501 'Delta ' // cd_delstr, pendval-pstartval, cd_budgetunits, & 502 '# of bergs', kbergs 503 ELSE 504 WRITE(numicb,100) cd_budgetstr // ' state:', & 505 cd_startstr // ' start', pstartval, cd_budgetunits, & 506 cd_endstr // ' end', pendval, cd_budgetunits, & 507 cd_delstr // 'Delta', pendval-pstartval, cd_budgetunits 508 ENDIF 506 509 100 FORMAT(a19,3(a18,"=",es14.7,x,a2,:,","),a12,i8) 507 510 END SUBROUTINE report_state … … 509 512 !!------------------------------------------------------------------------- 510 513 511 SUBROUTINE report_consistant( budgetstr,budgetunits,startstr,startval,endstr,endval)514 SUBROUTINE report_consistant( cd_budgetstr, cd_budgetunits, cd_startstr, pstartval, cd_endstr, pendval) 512 515 ! Arguments 513 CHARACTER*(*), INTENT(in) :: budgetstr, budgetunits, startstr,endstr514 REAL(wp), INTENT(in) :: startval,endval515 516 WRITE(numicb,200) budgetstr//' check:',&517 startstr,startval,budgetunits, &518 endstr,endval,budgetunits, &519 'error', (endval-startval)/((endval+startval)+1e-30),'nd'516 CHARACTER*(*), INTENT(in) :: cd_budgetstr, cd_budgetunits, cd_startstr, cd_endstr 517 REAL(wp), INTENT(in) :: pstartval, pendval 518 519 WRITE(numicb,200) cd_budgetstr // ' check:', & 520 cd_startstr, pstartval, cd_budgetunits, & 521 cd_endstr, pendval, cd_budgetunits, & 522 'error', (pendval-pstartval)/((pendval+pstartval)+1e-30), 'nd' 520 523 200 FORMAT(a19,10(a18,"=",es14.7,x,a2,:,",")) 521 524 END SUBROUTINE report_consistant … … 523 526 !!------------------------------------------------------------------------- 524 527 525 SUBROUTINE report_budget(budgetstr,budgetunits,instr,inval,outstr,outval,delstr,startval,endval) 528 SUBROUTINE report_budget( cd_budgetstr, cd_budgetunits, cd_instr, pinval, cd_outstr, & 529 & poutval, cd_delstr, pstartval, pendval) 526 530 ! Arguments 527 CHARACTER*(*), INTENT(in) :: budgetstr, budgetunits, instr, outstr, delstr 528 REAL(wp), INTENT(in) :: inval, outval, startval, endval 529 530 WRITE(numicb,200) budgetstr//' budget:', & 531 instr//' in',inval,budgetunits, & 532 outstr//' out',outval,budgetunits, & 533 'Delta '//delstr,inval-outval,budgetunits, & 534 'error',((endval-startval)-(inval-outval))/max(1.e-30,max(abs(endval-startval),abs(inval-outval))),'nd' 531 CHARACTER*(*), INTENT(in) :: cd_budgetstr, cd_budgetunits, cd_instr, cd_outstr, cd_delstr 532 REAL(wp), INTENT(in) :: pinval, poutval, pstartval, pendval 533 ! Local variables 534 REAL(wp) :: zval 535 536 zval = ( ( pendval - pstartval ) - ( pinval - poutval ) ) / & 537 MAX( 1.e-30, MAX( abs( pendval - pstartval ) , ABS( pinval - poutval ) ) ) 538 539 WRITE(numicb,200) cd_budgetstr // ' budget:', & 540 cd_instr // ' in', pinval, cd_budgetunits, & 541 cd_outstr // ' out', poutval, cd_budgetunits, & 542 'Delta ' // cd_delstr, pinval-poutval, cd_budgetunits, & 543 'error', zval, 'nd' 535 544 200 FORMAT(a19,3(a18,"=",es14.7,x,a2,:,","),a8,"=",es10.3,x,a2) 536 545 END SUBROUTINE report_budget … … 538 547 !!------------------------------------------------------------------------- 539 548 540 SUBROUTINE report_istate( budgetstr,startstr,startval,endstr,endval,delstr)549 SUBROUTINE report_istate( cd_budgetstr, cd_startstr, pstartval, cd_endstr, pendval, cd_delstr) 541 550 ! Arguments 542 CHARACTER*(*), INTENT(in) :: budgetstr, startstr, endstr,delstr543 INTEGER, INTENT(in) :: startval,endval544 545 WRITE(numicb,100) budgetstr//' state:',&546 startstr//' start',startval, &547 endstr//' end',endval,&548 delstr//'Delta',endval-startval551 CHARACTER*(*), INTENT(in) :: cd_budgetstr, cd_startstr, cd_endstr, cd_delstr 552 INTEGER, INTENT(in) :: pstartval, pendval 553 554 WRITE(numicb,100) cd_budgetstr // ' state:', & 555 cd_startstr // ' start', pstartval, & 556 cd_endstr // ' end', pendval, & 557 cd_delstr // 'Delta', pendval-pstartval 549 558 100 FORMAT(a19,3(a18,"=",i14,x,:,",")) 550 559 END SUBROUTINE report_istate … … 552 561 !!------------------------------------------------------------------------- 553 562 554 SUBROUTINE report_ibudget(budgetstr,instr,inval,outstr,outval,delstr,startval,endval) 563 SUBROUTINE report_ibudget( cd_budgetstr, cd_instr, pinval, cd_outstr, poutval, & 564 & cd_delstr, pstartval, pendval) 555 565 ! Arguments 556 CHARACTER*(*), INTENT(in) :: budgetstr, instr, outstr,delstr557 INTEGER, INTENT(in) :: inval, outval, startval,endval558 559 WRITE(numicb,200) budgetstr//' budget:', &560 instr//' in',inval, &561 outstr//' out',outval, &562 'Delta '//delstr,inval-outval, &563 'error',((endval-startval)-(inval-outval))566 CHARACTER*(*), INTENT(in) :: cd_budgetstr, cd_instr, cd_outstr, cd_delstr 567 INTEGER, INTENT(in) :: pinval, poutval, pstartval, pendval 568 569 WRITE(numicb,200) cd_budgetstr // ' budget:', & 570 cd_instr // ' in', pinval, & 571 cd_outstr // ' out', poutval, & 572 'Delta ' // cd_delstr, pinval-poutval, & 573 'error', ( ( pendval - pstartval ) - ( pinval - poutval ) ) 564 574 200 FORMAT(a19,10(a18,"=",i14,x,:,",")) 565 575 END SUBROUTINE report_ibudget -
branches/2012/dev_r3337_NOCS10_ICB/NEMOGCM/NEMO/OPA_SRC/ICB/icbdyn.F90
r3339 r3359 13 13 !! - ! the broken one) 14 14 !!---------------------------------------------------------------------- 15 !!----------------------------------------------------------------------16 !! icb_init : initialise17 !! icb_gen : generate test icebergs18 !! icb_nam : read iceberg namelist19 !!----------------------------------------------------------------------20 15 USE par_oce ! NEMO parameters 21 16 USE dom_oce ! NEMO ocean domain … … 29 24 PRIVATE 30 25 31 PUBLIC evolve_icebergs ! routine called in xxx.F90 module26 PUBLIC evolve_icebergs ! routine called in icbrun.F90 module 32 27 33 28 CONTAINS … … 42 37 !!---------------------------------------------------------------------- 43 38 ! 44 REAL(wp) :: uvel1 , vvel1 , u1, v1, ax1, ay1, xi1 ,yj145 REAL(wp) :: uvel2 , vvel2 , u2, v2, ax2, ay2, xi2 ,yj246 REAL(wp) :: uvel3 , vvel3 , u3, v3, ax3, ay3, xi3 ,yj347 REAL(wp) :: uvel4 , vvel4 , u4, v4, ax4, ay4, xi4 ,yj448 REAL(wp) :: uvel_n, vvel_n , xi_n,yj_n49 REAL(wp) :: zdt, zdt_2, zdt_6, e1,e250 LOGICAL :: bounced39 REAL(wp) :: zuvel1 , zvvel1 , zu1, zv1, zax1, zay1, zxi1 , zyj1 40 REAL(wp) :: zuvel2 , zvvel2 , zu2, zv2, zax2, zay2, zxi2 , zyj2 41 REAL(wp) :: zuvel3 , zvvel3 , zu3, zv3, zax3, zay3, zxi3 , zyj3 42 REAL(wp) :: zuvel4 , zvvel4 , zu4, zv4, zax4, zay4, zxi4 , zyj4 43 REAL(wp) :: zuvel_n, zvvel_n , zxi_n, zyj_n 44 REAL(wp) :: zdt, zdt_2, zdt_6, ze1, ze2 45 LOGICAL :: ll_bounced 51 46 TYPE(iceberg), POINTER :: berg 52 47 TYPE(point) , POINTER :: pt … … 77 72 pt => berg%current_point 78 73 79 bounced = .false.74 ll_bounced = .false. 80 75 81 76 82 77 ! STEP 1 ! 83 78 ! ====== ! 84 xi1 = pt%xi ;uvel1 = pt%uvel !** X1 in (i,j) ; V1 in m/s85 yj1 = pt%yj ;vvel1 = pt%vvel79 zxi1 = pt%xi ; zuvel1 = pt%uvel !** X1 in (i,j) ; V1 in m/s 80 zyj1 = pt%yj ; zvvel1 = pt%vvel 86 81 87 82 88 83 ! !** A1 = A(X1,V1) 89 CALL accel( xi1, e1, uvel1, uvel1,ax1, &90 & berg , yj1, e2, vvel1, vvel1,ay1, zdt_2 )91 ! 92 u1 = uvel1 /e1 !** V1 in d(i,j)/dt93 v1 = vvel1 /e284 CALL accel( zxi1, ze1, zuvel1, zuvel1, zax1, & 85 & berg , zyj1, ze2, zvvel1, zvvel1, zay1, zdt_2 ) 86 ! 87 zu1 = zuvel1 / ze1 !** V1 in d(i,j)/dt 88 zv1 = zvvel1 / ze2 94 89 95 90 ! STEP 2 ! … … 97 92 ! !** X2 = X1+dt/2*V1 ; V2 = V1+dt/2*A1 98 93 ! position using di/dt & djdt ! V2 in m/s 99 xi2 = xi1 + zdt_2 * u1 ; uvel2 = uvel1 + zdt_2 * ax1 100 yj2 = yj1 + zdt_2 * v1 ; vvel2 = vvel1 + zdt_2 * ay1 101 ! 102 CALL adjust_to_ground( xi2, xi1, u1, yj2, yj1, v1, bounced ) 94 zxi2 = zxi1 + zdt_2 * zu1 ; zuvel2 = zuvel1 + zdt_2 * zax1 95 zyj2 = zyj1 + zdt_2 * zv1 ; zvvel2 = zvvel1 + zdt_2 * zay1 96 ! 97 CALL adjust_to_ground( zxi2, zxi1, zu1, & 98 & zyj2, zyj1, zv1, ll_bounced ) 103 99 104 100 ! !** A2 = A(X2,V2) 105 CALL accel( xi2, e1, uvel2, uvel1,ax2, &106 & berg , yj2, e2, vvel2, vvel1,ay2, zdt_2 )107 ! 108 u2 = uvel2 /e1 !** V2 in d(i,j)/dt109 v2 = vvel2 /e2101 CALL accel( zxi2, ze1, zuvel2, zuvel1, zax2, & 102 & berg , zyj2, ze2, zvvel2, zvvel1, zay2, zdt_2 ) 103 ! 104 zu2 = zuvel2 / ze1 !** V2 in d(i,j)/dt 105 zv2 = zvvel2 / ze2 110 106 ! 111 107 ! STEP 3 ! 112 108 ! ====== ! 113 109 ! !** X3 = X1+dt/2*V2 ; V3 = V1+dt/2*A2; A3=A(X3) 114 xi3 = xi1 + zdt_2 * u2 ; uvel3 = uvel1 + zdt_2 * ax2 115 yj3 = yj1 + zdt_2 * v2 ; vvel3 = vvel1 + zdt_2 * ay2 116 ! 117 CALL adjust_to_ground( xi3, xi1, u3, yj3, yj1, v3, bounced ) 110 zxi3 = zxi1 + zdt_2 * zu2 ; zuvel3 = zuvel1 + zdt_2 * zax2 111 zyj3 = zyj1 + zdt_2 * zv2 ; zvvel3 = zvvel1 + zdt_2 * zay2 112 ! 113 CALL adjust_to_ground( zxi3, zxi1, zu3, & 114 & zyj3, zyj1, zv3, ll_bounced ) 118 115 119 116 ! !** A3 = A(X3,V3) 120 CALL accel( xi3, e1, uvel3, uvel1,ax3, &121 & berg , yj3, e2, vvel3, vvel1,ay3, zdt )122 ! 123 u3 = uvel3 /e1 !** V3 in d(i,j)/dt124 v3 = vvel3 /e2117 CALL accel( zxi3, ze1, zuvel3, zuvel1, zax3, & 118 & berg , zyj3, ze2, zvvel3, zvvel1, zay3, zdt ) 119 ! 120 zu3 = zuvel3 / ze1 !** V3 in d(i,j)/dt 121 zv3 = zvvel3 / ze2 125 122 126 123 ! STEP 4 ! 127 124 ! ====== ! 128 125 ! !** X4 = X1+dt*V3 ; V4 = V1+dt*A3 129 xi4 = xi1 + zdt * u3 ; uvel4 = uvel1 + zdt * ax3 130 yj4 = yj1 + zdt * v3 ; vvel4 = vvel1 + zdt * ay3 131 132 CALL adjust_to_ground( xi4, xi1, u4, yj4, yj1, v4, bounced ) 126 zxi4 = zxi1 + zdt * zu3 ; zuvel4 = zuvel1 + zdt * zax3 127 zyj4 = zyj1 + zdt * zv3 ; zvvel4 = zvvel1 + zdt * zay3 128 129 CALL adjust_to_ground( zxi4, zxi1, zu4, & 130 & zyj4, zyj1, zv4, ll_bounced ) 133 131 134 132 ! !** A4 = A(X4,V4) 135 CALL accel( xi4, e1, uvel4, uvel1,ax4, &136 & berg , yj4, e2, vvel4, vvel1,ay4, zdt )137 138 u4 = uvel4 /e1 !** V4 in d(i,j)/dt139 v4 = vvel4 /e2133 CALL accel( zxi4, ze1, zuvel4, zuvel1, zax4, & 134 & berg , zyj4, ze2, zvvel4, zvvel1, zay4, zdt ) 135 136 zu4 = zuvel4 / ze1 !** V4 in d(i,j)/dt 137 zv4 = zvvel4 / ze2 140 138 141 139 ! FINAL STEP ! … … 143 141 ! !** Xn = X1+dt*(V1+2*V2+2*V3+V4)/6 144 142 ! !** Vn = V1+dt*(A1+2*A2+2*A3+A4)/6 145 xi_n = pt%xi + zdt_6 * ( u1 + 2.*(u2 + u3 ) + u4 ) 146 yj_n = pt%yj + zdt_6 * ( v1 + 2.*(v2 + v3 ) + v4 ) 147 uvel_n = pt%uvel + zdt_6 * ( ax1 + 2.*(ax2 + ax3) + ax4 ) 148 vvel_n = pt%vvel + zdt_6 * ( ay1 + 2.*(ay2 + ay3) + ay4 ) 149 150 CALL adjust_to_ground( xi_n, xi1, uvel_n, yj_n, yj1, vvel_n, bounced ) 151 152 pt%uvel = uvel_n !** save in berg structure 153 pt%vvel = vvel_n 154 pt%xi = xi_n 155 pt%yj = yj_n 156 157 ! sga - update actual position 143 zxi_n = pt%xi + zdt_6 * ( zu1 + 2.*(zu2 + zu3 ) + zu4 ) 144 zyj_n = pt%yj + zdt_6 * ( zv1 + 2.*(zv2 + zv3 ) + zv4 ) 145 zuvel_n = pt%uvel + zdt_6 * ( zax1 + 2.*(zax2 + zax3) + zax4 ) 146 zvvel_n = pt%vvel + zdt_6 * ( zay1 + 2.*(zay2 + zay3) + zay4 ) 147 148 CALL adjust_to_ground( zxi_n, zxi1, zuvel_n, & 149 & zyj_n, zyj1, zvvel_n, ll_bounced ) 150 151 pt%uvel = zuvel_n !** save in berg structure 152 pt%vvel = zvvel_n 153 pt%xi = zxi_n 154 pt%yj = zyj_n 155 156 ! update actual position 158 157 pt%lon = bilin_x(glamt, pt%xi, pt%yj ) 159 158 pt%lat = bilin(gphit, pt%xi, pt%yj, 'T', 0, 0 ) … … 166 165 167 166 168 SUBROUTINE adjust_to_ground( pi, pi0, pu, pj, pj0, pv, bounced ) 167 SUBROUTINE adjust_to_ground( pi, pi0, pu, & 168 & pj, pj0, pv, ld_bounced ) 169 169 !!---------------------------------------------------------------------- 170 170 !! *** ROUTINE adjust_to_ground *** … … 177 177 REAL(wp), INTENT(in ) :: pi0, pj0 ! previous iceberg position 178 178 REAL(wp), INTENT(inout) :: pu , pv ! current iceberg velocities 179 LOGICAL , INTENT( out) :: bounced ! bounced indicator 180 ! 181 REAL(wp), PARAMETER :: posn_eps = 0.05_wp ! bouncing distance in (i,j) referential 179 LOGICAL , INTENT( out) :: ld_bounced ! bounced indicator 182 180 ! 183 181 INTEGER :: ii, ii0 184 182 INTEGER :: ij, ij0 185 INTEGER :: bounce_method 186 !!---------------------------------------------------------------------- 187 ! 188 bounced = .false. 189 bounce_method = 2 183 INTEGER :: ibounce_method 184 !!---------------------------------------------------------------------- 185 ! 186 ld_bounced = .FALSE. 190 187 ! 191 188 ii0 = INT( pi0+0.5 ) ; ij0 = INT( pj0+0.5 ) ! initial gridpoint position (T-cell) … … 204 201 ! From here, berg have reach land: treat grounding/bouncing 205 202 ! ------------------------------- 206 bounced = .true.203 ld_bounced = .TRUE. 207 204 208 205 !! not obvious what should happen now … … 213 210 !! of travel to zero; at a coastal boundary this has the effect of sliding the berg along the coast 214 211 215 SELECT CASE ( bounce_method ) 212 ibounce_method = 2 213 SELECT CASE ( ibounce_method ) 216 214 CASE ( 1 ) 217 215 pi = pi0 … … 233 231 234 232 235 SUBROUTINE accel( xi, e1, uvel, uvel0,ax, &236 & berg , yj, e2, vvel, vvel0, ay,dt )233 SUBROUTINE accel( pxi, pe1, puvel, puvel0, pax, & 234 & berg , pyj, pe2, pvvel, pvvel0, pay, pdt ) 237 235 !!---------------------------------------------------------------------- 238 236 !! *** ROUTINE accel *** … … 243 241 !!---------------------------------------------------------------------- 244 242 TYPE(iceberg ), POINTER :: berg 245 REAL(wp), INTENT(in ) :: xi ,yj ! berg position in (i,j) referential246 REAL(wp), INTENT(in ) :: uvel ,vvel ! berg velocity [m/s]247 REAL(wp), INTENT(in ) :: uvel0,vvel0 ! initial berg velocity [m/s]248 REAL(wp), INTENT( out) :: e1,e2 ! horizontal scale factor at (xi,yj)249 REAL(wp), INTENT(inout) :: ax,ay ! berg acceleration250 REAL(wp), INTENT(in ) :: dt ! berg time step251 ! 252 REAL(wp), PARAMETER :: alpha = 0._wp !253 REAL(wp), PARAMETER :: beta = 1._wp !254 REAL(wp), PARAMETER :: vel_lim =15._wp ! max allowed berg speed255 REAL(wp), PARAMETER :: accel_lim = 1.e-2_wp ! max allowed berg acceleration256 REAL(wp), PARAMETER :: Cr0 = 0.06_wp !243 REAL(wp), INTENT(in ) :: pxi , pyj ! berg position in (i,j) referential 244 REAL(wp), INTENT(in ) :: puvel , pvvel ! berg velocity [m/s] 245 REAL(wp), INTENT(in ) :: puvel0, pvvel0 ! initial berg velocity [m/s] 246 REAL(wp), INTENT( out) :: pe1, pe2 ! horizontal scale factor at (xi,yj) 247 REAL(wp), INTENT(inout) :: pax, pay ! berg acceleration 248 REAL(wp), INTENT(in ) :: pdt ! berg time step 249 ! 250 REAL(wp), PARAMETER :: pp_alpha = 0._wp ! 251 REAL(wp), PARAMETER :: pp_beta = 1._wp ! 252 REAL(wp), PARAMETER :: pp_vel_lim =15._wp ! max allowed berg speed 253 REAL(wp), PARAMETER :: pp_accel_lim = 1.e-2_wp ! max allowed berg acceleration 254 REAL(wp), PARAMETER :: pp_Cr0 = 0.06_wp ! 257 255 ! 258 256 INTEGER :: itloop 259 REAL(wp) :: uo, vo, ui, vi, ua, va, uwave, vwave, ssh_x, ssh_y, sst, cn,hi260 REAL(wp) :: zff, T, D, W, L, M,F261 REAL(wp) :: drag_ocn, drag_atm, drag_ice,wave_rad262 REAL(wp) :: c_ocn, c_atm, c_ice263 REAL(wp) :: ampl, wmod, Cr, Lwavelength, Lcutoff,Ltop264 REAL(wp) :: lambda, detA, A11, A12, axe, aye,D_hi265 REAL(wp) :: uveln, vveln, us, vs, speed, loc_dx, new_speed257 REAL(wp) :: zuo, zvo, zui, zvi, zua, zva, zuwave, zvwave, zssh_x, zssh_y, zsst, zcn, zhi 258 REAL(wp) :: zff, zT, zD, zW, zL, zM, zF 259 REAL(wp) :: zdrag_ocn, zdrag_atm, zdrag_ice, zwave_rad 260 REAL(wp) :: z_ocn, z_atm, z_ice 261 REAL(wp) :: zampl, zwmod, zCr, zLwavelength, zLcutoff, zLtop 262 REAL(wp) :: zlambda, zdetA, zA11, zA12, zaxe, zaye, zD_hi 263 REAL(wp) :: zuveln, zvveln, zus, zvs, zspeed, zloc_dx, zspeed_new 266 264 !!---------------------------------------------------------------------- 267 265 268 266 ! Interpolate gridded fields to berg 269 knberg = berg%number(1)270 CALL interp_flds( xi, e1, uo, ui, ua,ssh_x, &271 & yj, e2, vo, vi, va, ssh_y, sst, cn,hi, zff )272 273 M = berg%current_point%mass274 T = berg%current_point%thickness ! total thickness275 D = ( rn_rho_bergs / rho_seawater ) *T ! draught (keel depth)276 F = T -D ! freeboard277 W = berg%current_point%width278 L = berg%current_point%length279 280 hi = MIN( hi ,D )281 D_hi = MAX( 0._wp, D-hi )267 nknberg = berg%number(1) 268 CALL interp_flds( pxi, pe1, zuo, zui, zua, zssh_x, & 269 & pyj, pe2, zvo, zvi, zva, zssh_y, zsst, zcn, zhi, zff ) 270 271 zM = berg%current_point%mass 272 zT = berg%current_point%thickness ! total thickness 273 zD = ( rn_rho_bergs / pp_rho_seawater ) * zT ! draught (keel depth) 274 zF = zT - zD ! freeboard 275 zW = berg%current_point%width 276 zL = berg%current_point%length 277 278 zhi = MIN( zhi , zD ) 279 zD_hi = MAX( 0._wp, zD-zhi ) 282 280 283 281 ! Wave radiation 284 uwave = ua - uo ; vwave = va -vo ! Use wind speed rel. to ocean for wave model285 wmod = uwave*uwave + vwave*vwave ! The wave amplitude and length depend on the current;282 zuwave = zua - zuo ; zvwave = zva - zvo ! Use wind speed rel. to ocean for wave model 283 zwmod = zuwave*zuwave + zvwave*zvwave ! The wave amplitude and length depend on the current; 286 284 ! ! wind speed relative to the ocean. Actually wmod is wmod**2 here. 287 ampl = 0.5 * 0.02025 *wmod ! This is "a", the wave amplitude288 Lwavelength = 0.32 *wmod ! Surface wave length fitted to data in table at285 zampl = 0.5 * 0.02025 * zwmod ! This is "a", the wave amplitude 286 zLwavelength = 0.32 * zwmod ! Surface wave length fitted to data in table at 289 287 ! ! http://www4.ncsu.edu/eos/users/c/ceknowle/public/chapter10/part2.html 290 Lcutoff = 0.125 *Lwavelength291 Ltop = 0.25 *Lwavelength292 Cr = Cr0 * MIN( MAX( 0., (L-Lcutoff) / ((Ltop-Lcutoff)+1.e-30)) , 1.) ! Wave radiation coefficient288 zLcutoff = 0.125 * zLwavelength 289 zLtop = 0.25 * zLwavelength 290 zCr = pp_Cr0 * MIN( MAX( 0., (zL-zLcutoff) / ((zLtop-zLcutoff)+1.e-30)) , 1.) ! Wave radiation coefficient 293 291 ! ! fitted to graph from Carrieres et al., POAC Drift Model. 294 wave_rad = 0.5 * rho_seawater / M * Cr * grav * ampl * MIN( ampl,F ) * (2.*W*L) / (W+L)295 wmod = SQRT( ua*ua + va*va ) ! Wind speed296 IF( wmod /= 0._wp ) THEN297 uwave = ua/wmod ! Wave radiation force acts in wind direction ... !!gm this should be the wind rel. to ocean ?298 vwave = va/wmod292 zwave_rad = 0.5 * pp_rho_seawater / zM * zCr * grav * zampl * MIN( zampl,zF ) * (2.*zW*zL) / (zW+zL) 293 zwmod = SQRT( zua*zua + zva*zva ) ! Wind speed 294 IF( zwmod /= 0._wp ) THEN 295 zuwave = zua/zwmod ! Wave radiation force acts in wind direction ... !!gm this should be the wind rel. to ocean ? 296 zvwave = zva/zwmod 299 297 ELSE 300 uwave = 0. ; vwave=0. ;wave_rad=0. ! ... and only when wind is present. !!gm wave_rad=0. is useless298 zuwave = 0. ; zvwave=0. ; zwave_rad=0. ! ... and only when wind is present. !!gm wave_rad=0. is useless 301 299 ENDIF 302 300 303 301 ! Weighted drag coefficients 304 c_ocn = rho_seawater / M * (0.5*Cd_wv*W*(D_hi)+Cd_wh*W*L)305 c_atm = rho_air / M * (0.5*Cd_av*W*F +Cd_ah*W*L)306 c_ice = rho_ice / M * (0.5*Cd_iv*W*hi )307 IF( abs( ui) + abs(vi) == 0._wp ) c_ice = 0._wp308 309 uveln = uvel ; vveln =vvel ! Copy starting uvel, vvel302 z_ocn = pp_rho_seawater / zM * (0.5*pp_Cd_wv*zW*(zD_hi)+pp_Cd_wh*zW*zL) 303 z_atm = pp_rho_air / zM * (0.5*pp_Cd_av*zW*zF +pp_Cd_ah*zW*zL) 304 z_ice = pp_rho_ice / zM * (0.5*pp_Cd_iv*zW*zhi ) 305 IF( abs(zui) + abs(zvi) == 0._wp ) z_ice = 0._wp 306 307 zuveln = puvel ; zvveln = pvvel ! Copy starting uvel, vvel 310 308 ! 311 309 DO itloop = 1, 2 ! Iterate on drag coefficients 312 310 ! 313 us = 0.5 * ( uveln +uvel )314 vs = 0.5 * ( vveln +vvel )315 drag_ocn = c_ocn * SQRT( (us-uo)*(us-uo) + (vs-vo)*(vs-vo) )316 drag_atm = c_atm * SQRT( (us-ua)*(us-ua) + (vs-va)*(vs-va) )317 drag_ice = c_ice * SQRT( (us-ui)*(us-ui) + (vs-vi)*(vs-vi) )311 zus = 0.5 * ( zuveln + puvel ) 312 zvs = 0.5 * ( zvveln + pvvel ) 313 zdrag_ocn = z_ocn * SQRT( (zus-zuo)*(zus-zuo) + (zvs-zvo)*(zvs-zvo) ) 314 zdrag_atm = z_atm * SQRT( (zus-zua)*(zus-zua) + (zvs-zva)*(zvs-zva) ) 315 zdrag_ice = z_ice * SQRT( (zus-zui)*(zus-zui) + (zvs-zvi)*(zvs-zvi) ) 318 316 ! 319 317 ! Explicit accelerations 320 ! axe= zff*vvel -grav*ssh_x +wave_rad*uwave &321 ! - drag_ocn*(uvel-uo) -drag_atm*(uvel-ua) -drag_ice*(uvel-ui)322 ! aye=-zff*uvel -grav*ssh_y +wave_rad*vwave &323 ! - drag_ocn*(vvel-vo) -drag_atm*(vvel-va) -drag_ice*(vvel-vi)324 axe = -grav * ssh_x + wave_rad *uwave325 aye = -grav * ssh_y + wave_rad *vwave326 IF( alpha > 0._wp ) THEN ! If implicit, use time-level (n) rather than RK4 latest327 axe = axe + zff*vvel0328 aye = aye - zff*uvel0318 !zaxe= zff*pvvel -grav*zssh_x +zwave_rad*zuwave & 319 ! -zdrag_ocn*(puvel-zuo) -zdrag_atm*(puvel-zua) -zdrag_ice*(puvel-zui) 320 !zaye=-zff*puvel -grav*zssh_y +zwave_rad*zvwave & 321 ! -zdrag_ocn*(pvvel-zvo) -zdrag_atm*(pvvel-zva) -zdrag_ice*(pvvel-zvi) 322 zaxe = -grav * zssh_x + zwave_rad * zuwave 323 zaye = -grav * zssh_y + zwave_rad * zvwave 324 IF( pp_alpha > 0._wp ) THEN ! If implicit, use time-level (n) rather than RK4 latest 325 zaxe = zaxe + zff*pvvel0 326 zaye = zaye - zff*puvel0 329 327 else 330 axe=axe+zff*vvel331 aye=aye-zff*uvel328 zaxe=zaxe+zff*pvvel 329 zaye=zaye-zff*puvel 332 330 endif 333 if ( beta>0.) then ! If implicit, use time-level (n) rather than RK4 latest334 axe=axe-drag_ocn*(uvel0-uo) -drag_atm*(uvel0-ua) -drag_ice*(uvel0-ui)335 aye=aye-drag_ocn*(vvel0-vo) -drag_atm*(vvel0-va) -drag_ice*(vvel0-vi)331 if (pp_beta>0.) then ! If implicit, use time-level (n) rather than RK4 latest 332 zaxe=zaxe-zdrag_ocn*(puvel0-zuo) -zdrag_atm*(puvel0-zua) -zdrag_ice*(puvel0-zui) 333 zaye=zaye-zdrag_ocn*(pvvel0-zvo) -zdrag_atm*(pvvel0-zva) -zdrag_ice*(pvvel0-zvi) 336 334 else 337 axe=axe-drag_ocn*(uvel-uo) -drag_atm*(uvel-ua) -drag_ice*(uvel-ui)338 aye=aye-drag_ocn*(vvel-vo) -drag_atm*(vvel-va) -drag_ice*(vvel-vi)335 zaxe=zaxe-zdrag_ocn*(puvel-zuo) -zdrag_atm*(puvel-zua) -zdrag_ice*(puvel-zui) 336 zaye=zaye-zdrag_ocn*(pvvel-zvo) -zdrag_atm*(pvvel-zva) -zdrag_ice*(pvvel-zvi) 339 337 endif 340 338 341 339 ! Solve for implicit accelerations 342 IF( alpha +beta > 0._wp ) THEN343 lambda = drag_ocn + drag_atm +drag_ice344 A11 = 1.+beta*dt*lambda345 A12 = alpha*dt*zff346 detA = 1._wp / ( A11*A11 + A12*A12 )347 ax = detA * ( A11*axe+A12*aye )348 ay = detA * ( A11*aye-A12*axe )340 IF( pp_alpha + pp_beta > 0._wp ) THEN 341 zlambda = zdrag_ocn + zdrag_atm + zdrag_ice 342 zA11 = 1.+pp_beta*pdt*zlambda 343 zA12 = pp_alpha*pdt*zff 344 zdetA = 1._wp / ( zA11*zA11 + zA12*zA12 ) 345 pax = zdetA * ( zA11*zaxe+zA12*zaye ) 346 pay = zdetA * ( zA11*zaye-zA12*zaxe ) 349 347 else 350 ax=axe ; ay=aye348 pax=zaxe ; pay=zaye 351 349 endif 352 350 353 uveln = uvel0 + dt*ax354 vveln = vvel0 + dt*ay351 zuveln = puvel0 + pdt*pax 352 zvveln = pvvel0 + pdt*pay 355 353 ! 356 354 END DO ! itloop 357 355 358 356 IF( rn_speed_limit > 0.) THEN ! Limit speed of bergs based on a CFL criteria (if asked) 359 speed = SQRT( uveln*uveln + vveln*vveln ) ! Speed of berg360 IF( speed > 0.) THEN361 loc_dx = MIN( e1,e2 ) ! minimum grid spacing362 new_speed = loc_dx /dt * rn_speed_limit ! Speed limit as a factor of dx / dt363 IF( new_speed <speed ) THEN364 uveln = uveln * ( new_speed /speed ) ! Scale velocity to reduce speed365 vveln = vveln * ( new_speed /speed ) ! without changing the direction357 zspeed = SQRT( zuveln*zuveln + zvveln*zvveln ) ! Speed of berg 358 IF( zspeed > 0.) THEN 359 zloc_dx = MIN( pe1, pe2 ) ! minimum grid spacing 360 zspeed_new = zloc_dx / pdt * rn_speed_limit ! Speed limit as a factor of dx / dt 361 IF( zspeed_new < zspeed ) THEN 362 zuveln = zuveln * ( zspeed_new / zspeed ) ! Scale velocity to reduce speed 363 zvveln = zvveln * ( zspeed_new / zspeed ) ! without changing the direction 366 364 CALL speed_budget() 367 365 ENDIF … … 369 367 ENDIF 370 368 ! ! check the speed and acceleration limits 371 IF( ABS( uveln ) > vel_lim .OR. ABS( vveln ) >vel_lim ) &369 IF( ABS( zuveln ) > pp_vel_lim .OR. ABS( zvveln ) > pp_vel_lim ) & 372 370 WRITE(numicb,'("pe=",i3,x,a)') narea,'Dump triggered by excessive velocity' 373 IF( ABS( ax ) > accel_lim .OR. ABS( ay ) >accel_lim ) &371 IF( ABS( pax ) > pp_accel_lim .OR. ABS( pay ) > pp_accel_lim ) & 374 372 WRITE(numicb,'("pe=",i3,x,a)') narea,'Dump triggered by excessive acceleration' 375 373 ! -
branches/2012/dev_r3337_NOCS10_ICB/NEMOGCM/NEMO/OPA_SRC/ICB/icbini.F90
r3339 r3359 39 39 TYPE(FLD_N), PRIVATE :: sn_icb ! information about the calving file to be read 40 40 41 P RIVATE icb_alloc ! routine called in xxx.F90 module42 PUBLIC icb_ init ! routine called in xxx.F90 module43 P UBLIC icb_gen ! routine called in xxx.F90 module44 P UBLIC icb_nam ! routine called in xxx.F90 module41 PUBLIC icb_init ! routine called in nemogcm.F90 module 42 PUBLIC icb_gen ! routine called in icbclv.F90 module 43 PRIVATE icb_nam 44 PRIVATE icb_alloc 45 45 46 46 CONTAINS … … 100 100 ALLOCATE( ssh_e(0:jpi+1,0:jpj+1) , STAT=ill) 101 101 icb_alloc = icb_alloc + ill 102 ALLOCATE( initial_width (nclasses) , STAT=ill)103 icb_alloc = icb_alloc + ill 104 ALLOCATE( initial_length (nclasses) , STAT=ill)105 icb_alloc = icb_alloc + ill 106 ALLOCATE( p_calving(jpi,jpj) , STAT=ill)107 icb_alloc = icb_alloc + ill 108 ALLOCATE( p_calving_hflx(jpi,jpj) , STAT=ill)109 icb_alloc = icb_alloc + ill 110 111 ALLOCATE( icbfldpts(jpi) , STAT=ill)112 icb_alloc = icb_alloc + ill 113 ALLOCATE( icbflddest(jpi) , STAT=ill)114 icb_alloc = icb_alloc + ill 115 ALLOCATE( icbfldproc(jpni) , STAT=ill)102 ALLOCATE( first_width (nclasses) , STAT=ill) 103 icb_alloc = icb_alloc + ill 104 ALLOCATE( first_length (nclasses) , STAT=ill) 105 icb_alloc = icb_alloc + ill 106 ALLOCATE( src_calving(jpi,jpj) , STAT=ill) 107 icb_alloc = icb_alloc + ill 108 ALLOCATE( src_calving_hflx(jpi,jpj) , STAT=ill) 109 icb_alloc = icb_alloc + ill 110 111 ALLOCATE( nicbfldpts(jpi) , STAT=ill) 112 icb_alloc = icb_alloc + ill 113 ALLOCATE( nicbflddest(jpi) , STAT=ill) 114 icb_alloc = icb_alloc + ill 115 ALLOCATE( nicbfldproc(jpni) , STAT=ill) 116 116 icb_alloc = icb_alloc + ill 117 117 … … 121 121 END FUNCTION icb_alloc 122 122 123 SUBROUTINE icb_init( dt, kt )123 SUBROUTINE icb_init( pdt, kt ) 124 124 !!---------------------------------------------------------------------- 125 125 !! *** ROUTINE dom_init *** … … 129 129 !! ** Method : - blah blah 130 130 !!---------------------------------------------------------------------- 131 REAL(wp) , INTENT(in) :: dt! iceberg time-step (rdt*nn_fsbc)131 REAL(wp) , INTENT(in) :: pdt ! iceberg time-step (rdt*nn_fsbc) 132 132 INTEGER , INTENT(in) :: kt ! time step number 133 133 ! 134 INTEGER :: ji, jj, jn, i1, j1 135 INTEGER :: ii, inum, ivar, k1 134 INTEGER :: ji, jj, jn ! loop counters 135 INTEGER :: i1, i2, i3 ! dummy integers 136 INTEGER :: ii, inum, ivar 136 137 INTEGER :: istat1, istat2, istat3 137 CHARACTER(len=300) :: sdist138 CHARACTER(len=300) :: cl_sdist 138 139 !!---------------------------------------------------------------------- 139 140 … … 146 147 ! set parameters (mostly from namelist) 147 148 ! 148 berg_dt = dt149 initial_width (:) = SQRT( rn_initial_mass(:) / ( rn_LoW_ratio*rn_rho_bergs*rn_initial_thickness(:) ) )150 initial_length(:) = rn_LoW_ratio * initial_width(:)149 berg_dt = pdt 150 first_width (:) = SQRT( rn_initial_mass(:) / ( rn_LoW_ratio * rn_rho_bergs * rn_initial_thickness(:) ) ) 151 first_length(:) = rn_LoW_ratio * first_width(:) 151 152 152 153 berg_grid%calving (:,:) = 0._wp … … 154 155 berg_grid%stored_heat (:,:) = 0._wp 155 156 berg_grid%floating_melt(:,:) = 0._wp 156 berg_grid%maxclass (:,:)= nclasses157 berg_grid%maxclass (:,:) = nclasses 157 158 berg_grid%stored_ice (:,:,:) = 0._wp 158 159 berg_grid%tmp (:,:) = 0._wp 159 p_calving(:,:) = 0._wp160 p_calving_hflx(:,:) = 0._wp160 src_calving (:,:) = 0._wp 161 src_calving_hflx (:,:) = 0._wp 161 162 162 163 ! domain for icebergs … … 166 167 ! for the north fold we work out which points communicate by asking 167 168 ! lbc_lnk to pass processor number (valid even in single processor case) 168 ! borrow p_calving arrays for this169 ! borrow src_calving arrays for this 169 170 ! 170 171 ! pack i and j together using a scaling of a power of 10 171 icbpack = 10000172 IF( jpiglo >= icbpack ) CALL ctl_stop('icbini: processor index packing failure')173 icbfldproc(:) = -1172 nicbpack = 10000 173 IF( jpiglo >= nicbpack ) CALL ctl_stop('icbini: processor index packing failure') 174 nicbfldproc(:) = -1 174 175 175 176 DO jj = 1, jpj 176 177 DO ji = 1, jpi 177 p_calving_hflx(ji,jj) = narea178 p_calving(ji,jj) =icbpack*(njmpp+jj-1) + nimpp+ji-1178 src_calving_hflx(ji,jj) = narea 179 src_calving(ji,jj) = nicbpack*(njmpp+jj-1) + nimpp+ji-1 179 180 ENDDO 180 181 ENDDO 181 CALL lbc_lnk( p_calving_hflx, 'T', 1._wp )182 CALL lbc_lnk( p_calving , 'T', 1._wp )182 CALL lbc_lnk( src_calving_hflx, 'T', 1._wp ) 183 CALL lbc_lnk( src_calving , 'T', 1._wp ) 183 184 184 185 ! work out interior of processor from exchange array … … 186 187 ! last entry is right hand interior index 187 188 jj = jpj/2 188 icbdi = -1189 icbei = -1189 nicbdi = -1 190 nicbei = -1 190 191 DO ji = 1,jpi 191 k1 = INT( p_calving(ji,jj) )192 j1 = INT( k1/icbpack )193 i1 = k1 - j1*icbpack192 i3 = INT( src_calving(ji,jj) ) 193 i2 = INT( i3/nicbpack ) 194 i1 = i3 - i2*nicbpack 194 195 IF( i1 == nimpp+ji-1 ) THEN 195 IF( icbdi < 0 ) THEN196 icbdi = ji196 IF( nicbdi < 0 ) THEN 197 nicbdi = ji 197 198 ELSE 198 icbei = ji199 nicbei = ji 199 200 ENDIF 200 201 ENDIF … … 203 204 ! repeat for j direction 204 205 ji = jpi/2 205 icbdj = -1206 icbej = -1206 nicbdj = -1 207 nicbej = -1 207 208 DO jj = 1,jpj 208 k1 = INT( p_calving(ji,jj) )209 j1 = INT( k1/icbpack )210 i1 = k1 - j1*icbpack211 IF( j1== njmpp+jj-1 ) THEN212 IF( icbdj < 0 ) THEN213 icbdj = jj209 i3 = INT( src_calving(ji,jj) ) 210 i2 = INT( i3/nicbpack ) 211 i1 = i3 - i2*nicbpack 212 IF( i2 == njmpp+jj-1 ) THEN 213 IF( nicbdj < 0 ) THEN 214 nicbdj = jj 214 215 ELSE 215 icbej = jj216 nicbej = jj 216 217 ENDIF 217 218 ENDIF … … 219 220 220 221 ! special for east-west boundary exchange we save the destination index 221 i1 = MAX( icbdi-1, 1)222 k1 = INT( p_calving(i1,jpj/2) )223 jj = INT( k1/icbpack )224 icb_left = REAL( k1 -icbpack*jj, wp )225 i1 = MIN( icbei+1, jpi )226 k1 = INT( p_calving(i1,jpj/2) )227 jj = INT( k1/icbpack )228 icb_right = REAL( k1 -icbpack*jj, wp )222 i1 = MAX( nicbdi-1, 1) 223 i3 = INT( src_calving(i1,jpj/2) ) 224 jj = INT( i3/nicbpack ) 225 ricb_left = REAL( i3 - nicbpack*jj, wp ) 226 i1 = MIN( nicbei+1, jpi ) 227 i3 = INT( src_calving(i1,jpj/2) ) 228 jj = INT( i3/nicbpack ) 229 ricb_right = REAL( i3 - nicbpack*jj, wp ) 229 230 230 231 ! north fold 231 232 IF( npolj > 0 ) THEN 232 233 ! 233 ! icebergs in row icbej+1 get passed across fold234 icbfldpts(:) = INT( p_calving(:,icbej+1) )235 icbflddest(:) = INT( p_calving_hflx(:,icbej+1) )234 ! icebergs in row nicbej+1 get passed across fold 235 nicbfldpts(:) = INT( src_calving(:,nicbej+1) ) 236 nicbflddest(:) = INT( src_calving_hflx(:,nicbej+1) ) 236 237 ! 237 238 ! work out list of unique processors to talk to 238 DO ji = icbdi,icbei239 ii = icbflddest(ji)239 DO ji = nicbdi, nicbei 240 ii = nicbflddest(ji) 240 241 DO jn = 1, jpni 241 IF( icbfldproc(jn) == -1 ) THEN242 icbfldproc(jn) = ii242 IF( nicbfldproc(jn) == -1 ) THEN 243 nicbfldproc(jn) = ii 243 244 EXIT 244 245 ENDIF 245 IF( icbfldproc(jn) == ii ) EXIT246 IF( nicbfldproc(jn) == ii ) EXIT 246 247 ENDDO 247 248 ENDDO … … 253 254 WRITE(numicb,*) 'nldi, nlei ', nldi, nlei 254 255 WRITE(numicb,*) 'nldj, nlej ', nldj, nlej 255 WRITE(numicb,*) 'berg i interior ', icbdi,icbei256 WRITE(numicb,*) 'berg j interior ', icbdj,icbej257 WRITE(numicb,*) 'berg left ', icb_left258 WRITE(numicb,*) 'berg right ', icb_right256 WRITE(numicb,*) 'berg i interior ', nicbdi, nicbei 257 WRITE(numicb,*) 'berg j interior ', nicbdj, nicbej 258 WRITE(numicb,*) 'berg left ', ricb_left 259 WRITE(numicb,*) 'berg right ', ricb_right 259 260 jj = jpj/2 260 261 WRITE(numicb,*) "central j line:" 261 262 WRITE(numicb,*) "i processor" 262 WRITE(numicb,*) (INT( p_calving_hflx(ji,jj)), ji=1,jpi)263 WRITE(numicb,*) (INT(src_calving_hflx(ji,jj)), ji=1,jpi) 263 264 WRITE(numicb,*) "i point" 264 WRITE(numicb,*) (INT( p_calving(ji,jj)), ji=1,jpi)265 WRITE(numicb,*) (INT(src_calving(ji,jj)), ji=1,jpi) 265 266 ji = jpi/2 266 267 WRITE(numicb,*) "central i line:" 267 268 WRITE(numicb,*) "j processor" 268 WRITE(numicb,*) (INT( p_calving_hflx(ji,jj)), jj=1,jpj)269 WRITE(numicb,*) (INT(src_calving_hflx(ji,jj)), jj=1,jpj) 269 270 WRITE(numicb,*) "j point" 270 WRITE(numicb,*) (INT( p_calving(ji,jj)), jj=1,jpj)271 WRITE(numicb,*) (INT(src_calving(ji,jj)), jj=1,jpj) 271 272 IF( npolj > 0 ) THEN 272 273 WRITE(numicb,*) 'north fold destination points ' 273 WRITE(numicb,*) icbfldpts274 WRITE(numicb,*) nicbfldpts 274 275 WRITE(numicb,*) 'north fold destination procs ' 275 WRITE(numicb,*) icbflddest276 WRITE(numicb,*) nicbflddest 276 277 ENDIF 277 278 CALL flush(numicb) 278 279 ENDIF 279 280 280 p_calving(:,:) = 0._wp281 p_calving_hflx(:,:) = 0._wp281 src_calving(:,:) = 0._wp 282 src_calving_hflx(:,:) = 0._wp 282 283 283 284 ! assign each new iceberg with a unique number constructed from the processor number 284 285 ! and incremented by the total number of processors 285 kount_bergs(:) = 0286 kount_bergs(1) = narea - jpnij286 num_bergs(:) = 0 287 num_bergs(1) = narea - jpnij 287 288 288 289 ! when not generating test icebergs we need to setup calving file … … 290 291 ! 291 292 ! maximum distribution class array does not change in time so read it once 292 sdist = TRIM( cn_dir )//TRIM( sn_icb%clname )293 CALL iom_open ( sdist, inum ) ! open file293 cl_sdist = TRIM( cn_dir )//TRIM( sn_icb%clname ) 294 CALL iom_open ( cl_sdist, inum ) ! open file 294 295 ivar = iom_varid( inum, 'maxclass', ldstop=.FALSE. ) 295 296 IF( ivar > 0 ) THEN 296 CALL iom_get ( inum, jpdom_data, 'maxclass', p_calving ) ! read the max distribution array297 berg_grid%maxclass(:,:) = INT( p_calving )298 p_calving(:,:) = 0._wp297 CALL iom_get ( inum, jpdom_data, 'maxclass', src_calving ) ! read the max distribution array 298 berg_grid%maxclass(:,:) = INT( src_calving ) 299 src_calving(:,:) = 0._wp 299 300 ENDIF 300 301 CALL iom_close( inum ) ! close file … … 321 322 ! 322 323 CALL icebergs_read_restart() 323 restarted_bergs = .TRUE.324 l_restarted_bergs = .TRUE. 324 325 ENDIF 325 326 ENDIF … … 336 337 337 338 ! Local variables 338 INTEGER :: i, j, nbergs339 INTEGER :: ji, jj, ibergs 339 340 TYPE(iceberg) :: localberg ! NOT a pointer but an actual local variable 340 341 TYPE(point) :: localpt … … 358 359 ! use tmask here because tmask_i has been doctored on one side of the north fold line 359 360 360 DO j = icbdj,icbej361 DO i = icbdi,icbei362 IF (tmask( i,j,1) .GT. 0._wp .AND. &363 gphit( i,j) .GT. rn_test_box(3) .AND. gphit(i,j) .LT. rn_test_box(4) .AND. &364 glamt( i,j) .GT. rn_test_box(1) .AND. glamt(i,j) .LT. rn_test_box(2)) THEN361 DO jj = nicbdj,nicbej 362 DO ji = nicbdi,nicbei 363 IF (tmask(ji,jj,1) .GT. 0._wp .AND. & 364 gphit(ji,jj) .GT. rn_test_box(3) .AND. gphit(ji,jj) .LT. rn_test_box(4) .AND. & 365 glamt(ji,jj) .GT. rn_test_box(1) .AND. glamt(ji,jj) .LT. rn_test_box(2)) THEN 365 366 localberg%mass_scaling = rn_mass_scaling(iberg) 366 localpt%xi = REAL( nimpp+ i-1, wp )367 localpt%yj = REAL( njmpp+j -1, wp )367 localpt%xi = REAL( nimpp+ji-1, wp ) 368 localpt%yj = REAL( njmpp+jj-1, wp ) 368 369 localpt%lon = bilin(glamt, localpt%xi, localpt%yj, 'T', 0, 0 ) 369 370 localpt%lat = bilin(gphit, localpt%xi, localpt%yj, 'T', 0, 0 ) 370 371 localpt%mass = rn_initial_mass(iberg) 371 372 localpt%thickness = rn_initial_thickness(iberg) 372 localpt%width = initial_width(iberg)373 localpt%length = initial_length(iberg)373 localpt%width = first_width(iberg) 374 localpt%length = first_length(iberg) 374 375 localpt%year = iyr 375 376 localpt%day = FLOAT(iday)+(FLOAT(ihr)+FLOAT(imin)/60._wp)/24._wp … … 379 380 localpt%vvel = 0._wp 380 381 CALL increment_kounter() 381 localberg%number(:) = kount_bergs(:)382 localberg%number(:) = num_bergs(:) 382 383 call add_new_berg_to_list(localberg, localpt) 383 384 ENDIF … … 385 386 ENDDO 386 387 387 nbergs = count_bergs()388 IF( lk_mpp ) CALL mpp_sum( nbergs)389 WRITE(numicb,'(a,i6,a)') 'diamonds, icb_gen: ', nbergs,' were generated'388 ibergs = count_bergs() 389 IF( lk_mpp ) CALL mpp_sum(ibergs) 390 WRITE(numicb,'(a,i6,a)') 'diamonds, icb_gen: ',ibergs,' were generated' 390 391 391 392 END SUBROUTINE icb_gen … … 406 407 & ln_time_average_weight, nn_test_icebergs, rn_test_box, rn_speed_limit, & 407 408 & cn_dir, sn_icb 408 INTEGER :: k! loop counter409 REAL(wp) :: dfact409 INTEGER :: jn ! loop counter 410 REAL(wp) :: zfact 410 411 !!---------------------------------------------------------------------- 411 412 … … 431 432 ENDIF 432 433 433 dfact = SUM( rn_distribution )434 IF( dfact < 1._wp ) THEN435 IF( dfact <= 0._wp ) THEN434 zfact = SUM( rn_distribution ) 435 IF( zfact < 1._wp ) THEN 436 IF( zfact <= 0._wp ) THEN 436 437 CALL ctl_stop('icb_init: sum of berg distribution equal to zero') 437 438 ELSE 438 rn_distribution(:) = rn_distribution(:) / dfact439 rn_distribution(:) = rn_distribution(:) / zfact 439 440 CALL ctl_warn('icb_init: sum of berg input distribution not equal to one and so RESCALED') 440 441 ENDIF … … 448 449 WRITE(numout,*) ' Period between sampling of position for trajectory storage nn_sample_rate = ', nn_sample_rate 449 450 WRITE(numout,*) ' Mass thresholds between iceberg classes (kg) rn_initial_mass =' 450 DO k=1,nclasses451 WRITE(numout,'(a,f15.2)') ' ',rn_initial_mass( k)451 DO jn=1,nclasses 452 WRITE(numout,'(a,f15.2)') ' ',rn_initial_mass(jn) 452 453 ENDDO 453 454 WRITE(numout,*) ' Fraction of calving to apply to this class (non-dim) rn_distribution =' 454 DO k=1,nclasses455 WRITE(numout,'(a,f10.2)') ' ',rn_distribution( k)455 DO jn=1,nclasses 456 WRITE(numout,'(a,f10.2)') ' ',rn_distribution(jn) 456 457 ENDDO 457 458 WRITE(numout,*) ' Ratio between effective and real iceberg mass (non-dim) rn_mass_scaling = ' 458 DO k=1,nclasses459 WRITE(numout,'(a,f10.2)') ' ',rn_mass_scaling( k)459 DO jn=1,nclasses 460 WRITE(numout,'(a,f10.2)') ' ',rn_mass_scaling(jn) 460 461 ENDDO 461 462 WRITE(numout,*) ' Total thickness of newly calved bergs (m) rn_initial_thickness = ' 462 DO k=1,nclasses463 WRITE(numout,'(a,f10.2)') ' ',rn_initial_thickness( k)463 DO jn=1,nclasses 464 WRITE(numout,'(a,f10.2)') ' ',rn_initial_thickness(jn) 464 465 ENDDO 465 466 WRITE(numout,*) ' Timesteps between verbose messages nn_verbose_write = ', nn_verbose_write -
branches/2012/dev_r3337_NOCS10_ICB/NEMOGCM/NEMO/OPA_SRC/ICB/icblbc.F90
r3339 r3359 19 19 !! routines because they do not lie on regular jpi,jpj grids 20 20 !! Processor exchanges are handled as in lib_mpp whenever icebergs step 21 !! across boundary of interior domain ( icbdi-icbei, icbdj-icbej)21 !! across boundary of interior domain (nicbdi-nicbei, nicbdj-nicbej) 22 22 !! so that iceberg does not exist in more than one processor 23 23 !! North fold exchanges controlled by three arrays: 24 !! icbflddest - unique processor numbers that current one exchanges with25 !! icbfldproc - processor number that current grid point exchanges with26 !! icbfldpts - packed i,j point in exchanging processor24 !! nicbflddest - unique processor numbers that current one exchanges with 25 !! nicbfldproc - processor number that current grid point exchanges with 26 !! nicbfldpts - packed i,j point in exchanging processor 27 27 !!---------------------------------------------------------------------- 28 28 … … 56 56 TYPE(buffer), POINTER :: obuffer_f=>NULL() , ibuffer_f=>NULL() 57 57 58 INTEGER, PARAMETER, PUBLIC :: delta_buf = 25 ! Size by which to increment buffers59 INTEGER, PARAMETER, PUBLIC :: buffer_width = 15+nkounts ! items to store for each berg58 INTEGER, PARAMETER, PUBLIC :: pp_delta_buf = 25 ! Size by which to increment buffers 59 INTEGER, PARAMETER, PUBLIC :: pp_buffer_width = 15+nkounts ! items to store for each berg 60 60 61 61 #endif … … 88 88 TYPE(iceberg), POINTER :: this 89 89 TYPE(point) , POINTER :: pt 90 INTEGER :: i ne90 INTEGER :: iine 91 91 92 92 !! periodic east/west boundaries … … 98 98 DO WHILE( ASSOCIATED(this) ) 99 99 pt => this%current_point 100 i ne = INT( pt%xi + 0.5 )101 IF( i ne .GT. nimpp+icbei-1 ) THEN102 pt%xi = icb_right + MOD(pt%xi, 1._wp ) - 1._wp103 ELSE IF( i ne .LT. nimpp+icbdi-1 ) THEN104 pt%xi = icb_left + MOD(pt%xi, 1._wp )100 iine = INT( pt%xi + 0.5 ) 101 IF( iine .GT. nimpp+nicbei-1 ) THEN 102 pt%xi = ricb_right + MOD(pt%xi, 1._wp ) - 1._wp 103 ELSE IF( iine .LT. nimpp+nicbdi-1 ) THEN 104 pt%xi = ricb_left + MOD(pt%xi, 1._wp ) 105 105 ENDIF 106 106 this => this%next … … 136 136 TYPE(iceberg), POINTER :: this 137 137 TYPE(point) , POINTER :: pt 138 INTEGER :: i ne, jne, kpts139 INTEGER :: jiglo, jjglo138 INTEGER :: iine, ijne, ipts 139 INTEGER :: iiglo, ijglo 140 140 141 141 this => first_berg 142 142 DO WHILE( ASSOCIATED(this) ) 143 143 pt => this%current_point 144 jne = INT( pt%yj + 0.5 )145 IF( jne .GT. njmpp+icbej-1 ) THEN144 ijne = INT( pt%yj + 0.5 ) 145 IF( ijne .GT. njmpp+nicbej-1 ) THEN 146 146 ! 147 i ne = INT( pt%xi + 0.5 )148 kpts = icbfldpts (ine-nimpp+1)147 iine = INT( pt%xi + 0.5 ) 148 ipts = nicbfldpts (iine-nimpp+1) 149 149 ! 150 150 ! moving across the cut line means both position and 151 151 ! velocity must change 152 jjglo = INT( kpts/icbpack )153 jiglo = kpts - icbpack*jjglo154 pt%xi = jiglo - ( pt%xi - REAL(ine,wp) )155 pt%yj = jjglo - ( pt%yj - REAL(jne,wp) )152 ijglo = INT( ipts/nicbpack ) 153 iiglo = ipts - nicbpack*ijglo 154 pt%xi = iiglo - ( pt%xi - REAL(iine,wp) ) 155 pt%yj = ijglo - ( pt%yj - REAL(ijne,wp) ) 156 156 pt%uvel = -1._wp * pt%uvel 157 157 pt%vvel = -1._wp * pt%vvel … … 181 181 TYPE(iceberg) , POINTER :: tmpberg, this 182 182 TYPE(point) , POINTER :: pt 183 INTEGER :: nbergs_to_send_e, nbergs_to_send_w184 INTEGER :: nbergs_to_send_n, nbergs_to_send_s185 INTEGER :: nbergs_rcvd_from_e, nbergs_rcvd_from_w186 INTEGER :: nbergs_rcvd_from_n, nbergs_rcvd_from_s187 INTEGER :: i, nbergs_start, nbergs_end188 INTEGER :: i ne,jne189 INTEGER :: pe_N, pe_S, pe_W,pe_E190 REAL(wp), DIMENSION(2) :: ewbergs, webergs, nsbergs,snbergs191 INTEGER :: ml_req1, ml_req2, ml_req3,ml_req4192 INTEGER :: ml_req5, ml_req6, ml_req7, ml_req8,ml_err193 INTEGER, DIMENSION(MPI_STATUS_SIZE) :: ml_stat183 INTEGER :: ibergs_to_send_e, ibergs_to_send_w 184 INTEGER :: ibergs_to_send_n, ibergs_to_send_s 185 INTEGER :: ibergs_rcvd_from_e, ibergs_rcvd_from_w 186 INTEGER :: ibergs_rcvd_from_n, ibergs_rcvd_from_s 187 INTEGER :: i, ibergs_start, ibergs_end 188 INTEGER :: iine, ijne 189 INTEGER :: ipe_N, ipe_S, ipe_W, ipe_E 190 REAL(wp), DIMENSION(2) :: zewbergs, zwebergs, znsbergs, zsnbergs 191 INTEGER :: iml_req1, iml_req2, iml_req3, iml_req4 192 INTEGER :: iml_req5, iml_req6, iml_req7, iml_req8, iml_err 193 INTEGER, DIMENSION(MPI_STATUS_SIZE) :: iml_stat 194 194 195 195 ! set up indices of neighbouring processors 196 pe_N = -1197 pe_S = -1198 pe_W = -1199 pe_E = -1200 IF( nbondi .EQ. 0 .OR. nbondi .EQ. 1) pe_W = nowe201 IF( nbondi .EQ. -1 .OR. nbondi .EQ. 0) pe_E = noea202 IF( nbondj .EQ. 0 .OR. nbondj .EQ. 1) pe_S = noso203 IF( nbondj .EQ. -1 .OR. nbondj .EQ. 0) pe_N = nono196 ipe_N = -1 197 ipe_S = -1 198 ipe_W = -1 199 ipe_E = -1 200 IF( nbondi .EQ. 0 .OR. nbondi .EQ. 1) ipe_W = nowe 201 IF( nbondi .EQ. -1 .OR. nbondi .EQ. 0) ipe_E = noea 202 IF( nbondj .EQ. 0 .OR. nbondj .EQ. 1) ipe_S = noso 203 IF( nbondj .EQ. -1 .OR. nbondj .EQ. 0) ipe_N = nono 204 204 ! 205 205 ! at northern line of processors with north fold handle bergs differently 206 IF( npolj > 0 ) pe_N = -1206 IF( npolj > 0 ) ipe_N = -1 207 207 208 208 ! if there's only one processor in x direction then don't let mpp try to handle periodicity 209 209 IF( jpni == 1 ) THEN 210 pe_E = -1211 pe_W = -1210 ipe_E = -1 211 ipe_W = -1 212 212 ENDIF 213 213 214 214 IF( nn_verbose_level >= 2 ) THEN 215 WRITE(numicb,*) 'processor west : ', pe_W216 WRITE(numicb,*) 'processor east : ', pe_E217 WRITE(numicb,*) 'processor north : ', pe_N218 WRITE(numicb,*) 'processor south : ', pe_S215 WRITE(numicb,*) 'processor west : ', ipe_W 216 WRITE(numicb,*) 'processor east : ', ipe_E 217 WRITE(numicb,*) 'processor north : ', ipe_N 218 WRITE(numicb,*) 'processor south : ', ipe_S 219 219 WRITE(numicb,*) 'processor nimpp : ', nimpp 220 220 WRITE(numicb,*) 'processor njmpp : ', njmpp … … 234 234 IF( nn_verbose_level > 0 ) THEN 235 235 ! store the number of icebergs on this processor at start 236 nbergs_start = count_bergs()237 ENDIF 238 239 nbergs_to_send_e = 0240 nbergs_to_send_w = 0241 nbergs_to_send_n = 0242 nbergs_to_send_s = 0243 nbergs_rcvd_from_e = 0244 nbergs_rcvd_from_w = 0245 nbergs_rcvd_from_n = 0246 nbergs_rcvd_from_s = 0236 ibergs_start = count_bergs() 237 ENDIF 238 239 ibergs_to_send_e = 0 240 ibergs_to_send_w = 0 241 ibergs_to_send_n = 0 242 ibergs_to_send_s = 0 243 ibergs_rcvd_from_e = 0 244 ibergs_rcvd_from_w = 0 245 ibergs_rcvd_from_n = 0 246 ibergs_rcvd_from_s = 0 247 247 248 248 ! Find number of bergs that headed east/west … … 251 251 DO WHILE (ASSOCIATED(this)) 252 252 pt => this%current_point 253 i ne = INT( pt%xi + 0.5 )254 IF( pe_E >= 0 .AND. ine .GT. nimpp+icbei-1 ) THEN253 iine = INT( pt%xi + 0.5 ) 254 IF( ipe_E >= 0 .AND. iine .GT. nimpp+nicbei-1 ) THEN 255 255 tmpberg => this 256 256 this => this%next 257 nbergs_to_send_e = nbergs_to_send_e + 1257 ibergs_to_send_e = ibergs_to_send_e + 1 258 258 IF( nn_verbose_level >= 4 ) THEN 259 WRITE(numicb,*) 'bergstep ', ktberg,' packing berg ',tmpberg%number(:),' for transfer to east'259 WRITE(numicb,*) 'bergstep ',nktberg,' packing berg ',tmpberg%number(:),' for transfer to east' 260 260 CALL flush( numicb ) 261 261 ENDIF 262 262 ! deal with periodic case 263 tmpberg%current_point%xi = icb_right + MOD(tmpberg%current_point%xi, 1._wp ) - 1._wp263 tmpberg%current_point%xi = ricb_right + MOD(tmpberg%current_point%xi, 1._wp ) - 1._wp 264 264 ! now pack it into buffer and delete from list 265 CALL pack_berg_into_buffer( tmpberg, obuffer_e, nbergs_to_send_e)265 CALL pack_berg_into_buffer( tmpberg, obuffer_e, ibergs_to_send_e) 266 266 CALL delete_iceberg_from_list(first_berg, tmpberg) 267 ELSE IF( pe_W >= 0 .AND. ine .LT. nimpp+icbdi-1 ) THEN267 ELSE IF( ipe_W >= 0 .AND. iine .LT. nimpp+nicbdi-1 ) THEN 268 268 tmpberg => this 269 269 this => this%next 270 nbergs_to_send_w = nbergs_to_send_w + 1270 ibergs_to_send_w = ibergs_to_send_w + 1 271 271 IF( nn_verbose_level >= 4 ) THEN 272 WRITE(numicb,*) 'bergstep ', ktberg,' packing berg ',tmpberg%number(:),' for transfer to west'272 WRITE(numicb,*) 'bergstep ',nktberg,' packing berg ',tmpberg%number(:),' for transfer to west' 273 273 CALL flush( numicb ) 274 274 ENDIF 275 275 ! deal with periodic case 276 tmpberg%current_point%xi = icb_left + MOD(tmpberg%current_point%xi, 1._wp )276 tmpberg%current_point%xi = ricb_left + MOD(tmpberg%current_point%xi, 1._wp ) 277 277 ! now pack it into buffer and delete from list 278 CALL pack_berg_into_buffer( tmpberg, obuffer_w, nbergs_to_send_w)278 CALL pack_berg_into_buffer( tmpberg, obuffer_w, ibergs_to_send_w) 279 279 CALL delete_iceberg_from_list(first_berg, tmpberg) 280 280 ELSE … … 284 284 ENDIF 285 285 if( nn_verbose_level >= 3) then 286 write(numicb,*) 'bergstep ', ktberg,' send ew: ', nbergs_to_send_e, nbergs_to_send_w286 write(numicb,*) 'bergstep ',nktberg,' send ew: ', ibergs_to_send_e, ibergs_to_send_w 287 287 call flush(numicb) 288 288 endif … … 295 295 SELECT CASE ( nbondi ) 296 296 CASE( -1 ) 297 webergs(1) = nbergs_to_send_e298 CALL mppsend( 12, webergs(1), 1, pe_E,ml_req1)299 CALL mpprecv( 11, ewbergs(2), 1 )300 IF( l_isend ) CALL mpi_wait( ml_req1, ml_stat,ml_err )301 nbergs_rcvd_from_e = INT(ewbergs(2) )297 zwebergs(1) = ibergs_to_send_e 298 CALL mppsend( 12, zwebergs(1), 1, ipe_E, iml_req1) 299 CALL mpprecv( 11, zewbergs(2), 1 ) 300 IF( l_isend ) CALL mpi_wait( iml_req1, iml_stat, iml_err ) 301 ibergs_rcvd_from_e = INT( zewbergs(2) ) 302 302 CASE( 0 ) 303 ewbergs(1) = nbergs_to_send_w304 webergs(1) = nbergs_to_send_e305 CALL mppsend( 11, ewbergs(1), 1, pe_W,ml_req2)306 CALL mppsend( 12, webergs(1), 1, pe_E,ml_req3)307 CALL mpprecv( 11, ewbergs(2), 1 )308 CALL mpprecv( 12, webergs(2), 1 )309 IF( l_isend ) CALL mpi_wait( ml_req2, ml_stat,ml_err )310 IF( l_isend ) CALL mpi_wait( ml_req3, ml_stat,ml_err )311 nbergs_rcvd_from_e = INT(ewbergs(2) )312 nbergs_rcvd_from_w = INT(webergs(2) )303 zewbergs(1) = ibergs_to_send_w 304 zwebergs(1) = ibergs_to_send_e 305 CALL mppsend( 11, zewbergs(1), 1, ipe_W, iml_req2) 306 CALL mppsend( 12, zwebergs(1), 1, ipe_E, iml_req3) 307 CALL mpprecv( 11, zewbergs(2), 1 ) 308 CALL mpprecv( 12, zwebergs(2), 1 ) 309 IF( l_isend ) CALL mpi_wait( iml_req2, iml_stat, iml_err ) 310 IF( l_isend ) CALL mpi_wait( iml_req3, iml_stat, iml_err ) 311 ibergs_rcvd_from_e = INT( zewbergs(2) ) 312 ibergs_rcvd_from_w = INT( zwebergs(2) ) 313 313 CASE( 1 ) 314 ewbergs(1) = nbergs_to_send_w315 CALL mppsend( 11, ewbergs(1), 1, pe_W,ml_req4)316 CALL mpprecv( 12, webergs(2), 1 )317 IF( l_isend ) CALL mpi_wait( ml_req4, ml_stat,ml_err )318 nbergs_rcvd_from_w = INT(webergs(2) )314 zewbergs(1) = ibergs_to_send_w 315 CALL mppsend( 11, zewbergs(1), 1, ipe_W, iml_req4) 316 CALL mpprecv( 12, zwebergs(2), 1 ) 317 IF( l_isend ) CALL mpi_wait( iml_req4, iml_stat, iml_err ) 318 ibergs_rcvd_from_w = INT( zwebergs(2) ) 319 319 END SELECT 320 320 if( nn_verbose_level >= 3) then 321 write(numicb,*) 'bergstep ', ktberg,' recv ew: ', nbergs_rcvd_from_w, nbergs_rcvd_from_e321 write(numicb,*) 'bergstep ',nktberg,' recv ew: ', ibergs_rcvd_from_w, ibergs_rcvd_from_e 322 322 call flush(numicb) 323 323 endif … … 325 325 SELECT CASE ( nbondi ) 326 326 CASE( -1 ) 327 IF( nbergs_to_send_e > 0 ) CALL mppsend( 14, obuffer_e%data, nbergs_to_send_e*buffer_width, pe_E,ml_req1 )328 IF( nbergs_rcvd_from_e > 0 ) THEN329 CALL increase_ibuffer(ibuffer_e, nbergs_rcvd_from_e)330 CALL mpprecv( 13, ibuffer_e%data, nbergs_rcvd_from_e*buffer_width )331 ENDIF 332 IF( nbergs_to_send_e > 0 .AND. l_isend ) CALL mpi_wait( ml_req1, ml_stat,ml_err )333 DO i = 1, nbergs_rcvd_from_e327 IF( ibergs_to_send_e > 0 ) CALL mppsend( 14, obuffer_e%data, ibergs_to_send_e*pp_buffer_width, ipe_E, iml_req1 ) 328 IF( ibergs_rcvd_from_e > 0 ) THEN 329 CALL increase_ibuffer(ibuffer_e, ibergs_rcvd_from_e) 330 CALL mpprecv( 13, ibuffer_e%data, ibergs_rcvd_from_e*pp_buffer_width ) 331 ENDIF 332 IF( ibergs_to_send_e > 0 .AND. l_isend ) CALL mpi_wait( iml_req1, iml_stat, iml_err ) 333 DO i = 1, ibergs_rcvd_from_e 334 334 IF( nn_verbose_level >= 4 ) THEN 335 WRITE(numicb,*) 'bergstep ', ktberg,' unpacking berg ',INT(ibuffer_e%data(16,i)),' from east'335 WRITE(numicb,*) 'bergstep ',nktberg,' unpacking berg ',INT(ibuffer_e%data(16,i)),' from east' 336 336 CALL flush( numicb ) 337 337 ENDIF … … 339 339 ENDDO 340 340 CASE( 0 ) 341 IF( nbergs_to_send_w > 0 ) CALL mppsend( 13, obuffer_w%data, nbergs_to_send_w*buffer_width, pe_W,ml_req2 )342 IF( nbergs_to_send_e > 0 ) CALL mppsend( 14, obuffer_e%data, nbergs_to_send_e*buffer_width, pe_E,ml_req3 )343 IF( nbergs_rcvd_from_e > 0 ) THEN344 CALL increase_ibuffer(ibuffer_e, nbergs_rcvd_from_e)345 CALL mpprecv( 13, ibuffer_e%data, nbergs_rcvd_from_e*buffer_width )346 ENDIF 347 IF( nbergs_rcvd_from_w > 0 ) THEN348 CALL increase_ibuffer(ibuffer_w, nbergs_rcvd_from_w)349 CALL mpprecv( 14, ibuffer_w%data, nbergs_rcvd_from_w*buffer_width )350 ENDIF 351 IF( nbergs_to_send_w > 0 .AND. l_isend ) CALL mpi_wait( ml_req2, ml_stat,ml_err )352 IF( nbergs_to_send_e > 0 .AND. l_isend ) CALL mpi_wait( ml_req3, ml_stat,ml_err )353 DO i = 1, nbergs_rcvd_from_e341 IF( ibergs_to_send_w > 0 ) CALL mppsend( 13, obuffer_w%data, ibergs_to_send_w*pp_buffer_width, ipe_W, iml_req2 ) 342 IF( ibergs_to_send_e > 0 ) CALL mppsend( 14, obuffer_e%data, ibergs_to_send_e*pp_buffer_width, ipe_E, iml_req3 ) 343 IF( ibergs_rcvd_from_e > 0 ) THEN 344 CALL increase_ibuffer(ibuffer_e, ibergs_rcvd_from_e) 345 CALL mpprecv( 13, ibuffer_e%data, ibergs_rcvd_from_e*pp_buffer_width ) 346 ENDIF 347 IF( ibergs_rcvd_from_w > 0 ) THEN 348 CALL increase_ibuffer(ibuffer_w, ibergs_rcvd_from_w) 349 CALL mpprecv( 14, ibuffer_w%data, ibergs_rcvd_from_w*pp_buffer_width ) 350 ENDIF 351 IF( ibergs_to_send_w > 0 .AND. l_isend ) CALL mpi_wait( iml_req2, iml_stat, iml_err ) 352 IF( ibergs_to_send_e > 0 .AND. l_isend ) CALL mpi_wait( iml_req3, iml_stat, iml_err ) 353 DO i = 1, ibergs_rcvd_from_e 354 354 IF( nn_verbose_level >= 4 ) THEN 355 WRITE(numicb,*) 'bergstep ', ktberg,' unpacking berg ',INT(ibuffer_e%data(16,i)),' from east'355 WRITE(numicb,*) 'bergstep ',nktberg,' unpacking berg ',INT(ibuffer_e%data(16,i)),' from east' 356 356 CALL flush( numicb ) 357 357 ENDIF 358 358 CALL unpack_berg_from_buffer(first_berg, ibuffer_e, i) 359 359 ENDDO 360 DO i = 1, nbergs_rcvd_from_w360 DO i = 1, ibergs_rcvd_from_w 361 361 IF( nn_verbose_level >= 4 ) THEN 362 WRITE(numicb,*) 'bergstep ', ktberg,' unpacking berg ',INT(ibuffer_w%data(16,i)),' from west'362 WRITE(numicb,*) 'bergstep ',nktberg,' unpacking berg ',INT(ibuffer_w%data(16,i)),' from west' 363 363 CALL flush( numicb ) 364 364 ENDIF … … 366 366 ENDDO 367 367 CASE( 1 ) 368 IF( nbergs_to_send_w > 0 ) CALL mppsend( 13, obuffer_w%data, nbergs_to_send_w*buffer_width, pe_W,ml_req4 )369 IF( nbergs_rcvd_from_w > 0 ) THEN370 CALL increase_ibuffer(ibuffer_w, nbergs_rcvd_from_w)371 CALL mpprecv( 14, ibuffer_w%data, nbergs_rcvd_from_w*buffer_width )372 ENDIF 373 IF( nbergs_to_send_w > 0 .AND. l_isend ) CALL mpi_wait( ml_req4, ml_stat,ml_err )374 DO i = 1, nbergs_rcvd_from_w368 IF( ibergs_to_send_w > 0 ) CALL mppsend( 13, obuffer_w%data, ibergs_to_send_w*pp_buffer_width, ipe_W, iml_req4 ) 369 IF( ibergs_rcvd_from_w > 0 ) THEN 370 CALL increase_ibuffer(ibuffer_w, ibergs_rcvd_from_w) 371 CALL mpprecv( 14, ibuffer_w%data, ibergs_rcvd_from_w*pp_buffer_width ) 372 ENDIF 373 IF( ibergs_to_send_w > 0 .AND. l_isend ) CALL mpi_wait( iml_req4, iml_stat, iml_err ) 374 DO i = 1, ibergs_rcvd_from_w 375 375 IF( nn_verbose_level >= 4 ) THEN 376 WRITE(numicb,*) 'bergstep ', ktberg,' unpacking berg ',INT(ibuffer_w%data(16,i)),' from west'376 WRITE(numicb,*) 'bergstep ',nktberg,' unpacking berg ',INT(ibuffer_w%data(16,i)),' from west' 377 377 CALL flush( numicb ) 378 378 ENDIF … … 390 390 DO WHILE (ASSOCIATED(this)) 391 391 pt => this%current_point 392 jne = INT( pt%yj + 0.5 )393 IF( pe_N >= 0 .AND. jne .GT. njmpp+icbej-1 ) THEN392 ijne = INT( pt%yj + 0.5 ) 393 IF( ipe_N >= 0 .AND. ijne .GT. njmpp+nicbej-1 ) THEN 394 394 tmpberg => this 395 395 this => this%next 396 nbergs_to_send_n = nbergs_to_send_n + 1396 ibergs_to_send_n = ibergs_to_send_n + 1 397 397 IF( nn_verbose_level >= 4 ) THEN 398 WRITE(numicb,*) 'bergstep ', ktberg,' packing berg ',tmpberg%number(:),' for transfer to north'398 WRITE(numicb,*) 'bergstep ',nktberg,' packing berg ',tmpberg%number(:),' for transfer to north' 399 399 CALL flush( numicb ) 400 400 ENDIF 401 CALL pack_berg_into_buffer( tmpberg, obuffer_n, nbergs_to_send_n)401 CALL pack_berg_into_buffer( tmpberg, obuffer_n, ibergs_to_send_n) 402 402 CALL delete_iceberg_from_list(first_berg, tmpberg) 403 ELSE IF( pe_S >= 0 .AND. jne .LT. njmpp+icbdj-1 ) THEN403 ELSE IF( ipe_S >= 0 .AND. ijne .LT. njmpp+nicbdj-1 ) THEN 404 404 tmpberg => this 405 405 this => this%next 406 nbergs_to_send_s = nbergs_to_send_s + 1406 ibergs_to_send_s = ibergs_to_send_s + 1 407 407 IF( nn_verbose_level >= 4 ) THEN 408 WRITE(numicb,*) 'bergstep ', ktberg,' packing berg ',tmpberg%number(:),' for transfer to south'408 WRITE(numicb,*) 'bergstep ',nktberg,' packing berg ',tmpberg%number(:),' for transfer to south' 409 409 CALL flush( numicb ) 410 410 ENDIF 411 CALL pack_berg_into_buffer( tmpberg, obuffer_s, nbergs_to_send_s)411 CALL pack_berg_into_buffer( tmpberg, obuffer_s, ibergs_to_send_s) 412 412 CALL delete_iceberg_from_list(first_berg, tmpberg) 413 413 ELSE … … 417 417 ENDIF 418 418 if( nn_verbose_level >= 3) then 419 write(numicb,*) 'bergstep ', ktberg,' send ns: ', nbergs_to_send_n, nbergs_to_send_s419 write(numicb,*) 'bergstep ',nktberg,' send ns: ', ibergs_to_send_n, ibergs_to_send_s 420 420 call flush(numicb) 421 421 endif … … 426 426 SELECT CASE ( nbondj ) 427 427 CASE( -1 ) 428 snbergs(1) = nbergs_to_send_n429 CALL mppsend( 16, snbergs(1), 1, pe_N,ml_req1)430 CALL mpprecv( 15, nsbergs(2), 1 )431 IF( l_isend ) CALL mpi_wait( ml_req1, ml_stat,ml_err )432 nbergs_rcvd_from_n = INT(nsbergs(2) )428 zsnbergs(1) = ibergs_to_send_n 429 CALL mppsend( 16, zsnbergs(1), 1, ipe_N, iml_req1) 430 CALL mpprecv( 15, znsbergs(2), 1 ) 431 IF( l_isend ) CALL mpi_wait( iml_req1, iml_stat, iml_err ) 432 ibergs_rcvd_from_n = INT( znsbergs(2) ) 433 433 CASE( 0 ) 434 nsbergs(1) = nbergs_to_send_s435 snbergs(1) = nbergs_to_send_n436 CALL mppsend( 15, nsbergs(1), 1, pe_S,ml_req2)437 CALL mppsend( 16, snbergs(1), 1, pe_N,ml_req3)438 CALL mpprecv( 15, nsbergs(2), 1 )439 CALL mpprecv( 16, snbergs(2), 1 )440 IF( l_isend ) CALL mpi_wait( ml_req2, ml_stat,ml_err )441 IF( l_isend ) CALL mpi_wait( ml_req3, ml_stat,ml_err )442 nbergs_rcvd_from_n = INT(nsbergs(2) )443 nbergs_rcvd_from_s = INT(snbergs(2) )434 znsbergs(1) = ibergs_to_send_s 435 zsnbergs(1) = ibergs_to_send_n 436 CALL mppsend( 15, znsbergs(1), 1, ipe_S, iml_req2) 437 CALL mppsend( 16, zsnbergs(1), 1, ipe_N, iml_req3) 438 CALL mpprecv( 15, znsbergs(2), 1 ) 439 CALL mpprecv( 16, zsnbergs(2), 1 ) 440 IF( l_isend ) CALL mpi_wait( iml_req2, iml_stat, iml_err ) 441 IF( l_isend ) CALL mpi_wait( iml_req3, iml_stat, iml_err ) 442 ibergs_rcvd_from_n = INT( znsbergs(2) ) 443 ibergs_rcvd_from_s = INT( zsnbergs(2) ) 444 444 CASE( 1 ) 445 nsbergs(1) = nbergs_to_send_s446 CALL mppsend( 15, nsbergs(1), 1, pe_S,ml_req4)447 CALL mpprecv( 16, snbergs(2), 1 )448 IF( l_isend ) CALL mpi_wait( ml_req4, ml_stat,ml_err )449 nbergs_rcvd_from_s = INT(snbergs(2) )445 znsbergs(1) = ibergs_to_send_s 446 CALL mppsend( 15, znsbergs(1), 1, ipe_S, iml_req4) 447 CALL mpprecv( 16, zsnbergs(2), 1 ) 448 IF( l_isend ) CALL mpi_wait( iml_req4, iml_stat, iml_err ) 449 ibergs_rcvd_from_s = INT( zsnbergs(2) ) 450 450 END SELECT 451 451 if( nn_verbose_level >= 3) then 452 write(numicb,*) 'bergstep ', ktberg,' recv ns: ', nbergs_rcvd_from_s, nbergs_rcvd_from_n452 write(numicb,*) 'bergstep ',nktberg,' recv ns: ', ibergs_rcvd_from_s, ibergs_rcvd_from_n 453 453 call flush(numicb) 454 454 endif … … 456 456 SELECT CASE ( nbondj ) 457 457 CASE( -1 ) 458 IF( nbergs_to_send_n > 0 ) CALL mppsend( 18, obuffer_n%data, nbergs_to_send_n*buffer_width, pe_N,ml_req1 )459 IF( nbergs_rcvd_from_n > 0 ) THEN460 CALL increase_ibuffer(ibuffer_n, nbergs_rcvd_from_n)461 CALL mpprecv( 17, ibuffer_n%data, nbergs_rcvd_from_n*buffer_width )462 ENDIF 463 IF( nbergs_to_send_n > 0 .AND. l_isend ) CALL mpi_wait( ml_req1, ml_stat,ml_err )464 DO i = 1, nbergs_rcvd_from_n458 IF( ibergs_to_send_n > 0 ) CALL mppsend( 18, obuffer_n%data, ibergs_to_send_n*pp_buffer_width, ipe_N, iml_req1 ) 459 IF( ibergs_rcvd_from_n > 0 ) THEN 460 CALL increase_ibuffer(ibuffer_n, ibergs_rcvd_from_n) 461 CALL mpprecv( 17, ibuffer_n%data, ibergs_rcvd_from_n*pp_buffer_width ) 462 ENDIF 463 IF( ibergs_to_send_n > 0 .AND. l_isend ) CALL mpi_wait( iml_req1, iml_stat, iml_err ) 464 DO i = 1, ibergs_rcvd_from_n 465 465 IF( nn_verbose_level >= 4 ) THEN 466 WRITE(numicb,*) 'bergstep ', ktberg,' unpacking berg ',INT(ibuffer_n%data(16,i)),' from north'466 WRITE(numicb,*) 'bergstep ',nktberg,' unpacking berg ',INT(ibuffer_n%data(16,i)),' from north' 467 467 CALL flush( numicb ) 468 468 ENDIF … … 470 470 ENDDO 471 471 CASE( 0 ) 472 IF( nbergs_to_send_s > 0 ) CALL mppsend( 17, obuffer_s%data, nbergs_to_send_s*buffer_width, pe_S,ml_req2 )473 IF( nbergs_to_send_n > 0 ) CALL mppsend( 18, obuffer_n%data, nbergs_to_send_n*buffer_width, pe_N,ml_req3 )474 IF( nbergs_rcvd_from_n > 0 ) THEN475 CALL increase_ibuffer(ibuffer_n, nbergs_rcvd_from_n)476 CALL mpprecv( 17, ibuffer_n%data, nbergs_rcvd_from_n*buffer_width )477 ENDIF 478 IF( nbergs_rcvd_from_s > 0 ) THEN479 CALL increase_ibuffer(ibuffer_s, nbergs_rcvd_from_s)480 CALL mpprecv( 18, ibuffer_s%data, nbergs_rcvd_from_s*buffer_width )481 ENDIF 482 IF( nbergs_to_send_s > 0 .AND. l_isend ) CALL mpi_wait( ml_req2, ml_stat,ml_err )483 IF( nbergs_to_send_n > 0 .AND. l_isend ) CALL mpi_wait( ml_req3, ml_stat,ml_err )484 DO i = 1, nbergs_rcvd_from_n472 IF( ibergs_to_send_s > 0 ) CALL mppsend( 17, obuffer_s%data, ibergs_to_send_s*pp_buffer_width, ipe_S, iml_req2 ) 473 IF( ibergs_to_send_n > 0 ) CALL mppsend( 18, obuffer_n%data, ibergs_to_send_n*pp_buffer_width, ipe_N, iml_req3 ) 474 IF( ibergs_rcvd_from_n > 0 ) THEN 475 CALL increase_ibuffer(ibuffer_n, ibergs_rcvd_from_n) 476 CALL mpprecv( 17, ibuffer_n%data, ibergs_rcvd_from_n*pp_buffer_width ) 477 ENDIF 478 IF( ibergs_rcvd_from_s > 0 ) THEN 479 CALL increase_ibuffer(ibuffer_s, ibergs_rcvd_from_s) 480 CALL mpprecv( 18, ibuffer_s%data, ibergs_rcvd_from_s*pp_buffer_width ) 481 ENDIF 482 IF( ibergs_to_send_s > 0 .AND. l_isend ) CALL mpi_wait( iml_req2, iml_stat, iml_err ) 483 IF( ibergs_to_send_n > 0 .AND. l_isend ) CALL mpi_wait( iml_req3, iml_stat, iml_err ) 484 DO i = 1, ibergs_rcvd_from_n 485 485 IF( nn_verbose_level >= 4 ) THEN 486 WRITE(numicb,*) 'bergstep ', ktberg,' unpacking berg ',INT(ibuffer_n%data(16,i)),' from north'486 WRITE(numicb,*) 'bergstep ',nktberg,' unpacking berg ',INT(ibuffer_n%data(16,i)),' from north' 487 487 CALL flush( numicb ) 488 488 ENDIF 489 489 CALL unpack_berg_from_buffer(first_berg, ibuffer_n, i) 490 490 ENDDO 491 DO i = 1, nbergs_rcvd_from_s491 DO i = 1, ibergs_rcvd_from_s 492 492 IF( nn_verbose_level >= 4 ) THEN 493 WRITE(numicb,*) 'bergstep ', ktberg,' unpacking berg ',INT(ibuffer_s%data(16,i)),' from south'493 WRITE(numicb,*) 'bergstep ',nktberg,' unpacking berg ',INT(ibuffer_s%data(16,i)),' from south' 494 494 CALL flush( numicb ) 495 495 ENDIF … … 497 497 ENDDO 498 498 CASE( 1 ) 499 IF( nbergs_to_send_s > 0 ) CALL mppsend( 17, obuffer_s%data, nbergs_to_send_s*buffer_width, pe_S,ml_req4 )500 IF( nbergs_rcvd_from_s > 0 ) THEN501 CALL increase_ibuffer(ibuffer_s, nbergs_rcvd_from_s)502 CALL mpprecv( 18, ibuffer_s%data, nbergs_rcvd_from_s*buffer_width )503 ENDIF 504 IF( nbergs_to_send_s > 0 .AND. l_isend ) CALL mpi_wait( ml_req4, ml_stat,ml_err )505 DO i = 1, nbergs_rcvd_from_s499 IF( ibergs_to_send_s > 0 ) CALL mppsend( 17, obuffer_s%data, ibergs_to_send_s*pp_buffer_width, ipe_S, iml_req4 ) 500 IF( ibergs_rcvd_from_s > 0 ) THEN 501 CALL increase_ibuffer(ibuffer_s, ibergs_rcvd_from_s) 502 CALL mpprecv( 18, ibuffer_s%data, ibergs_rcvd_from_s*pp_buffer_width ) 503 ENDIF 504 IF( ibergs_to_send_s > 0 .AND. l_isend ) CALL mpi_wait( iml_req4, iml_stat, iml_err ) 505 DO i = 1, ibergs_rcvd_from_s 506 506 IF( nn_verbose_level >= 4 ) THEN 507 WRITE(numicb,*) 'bergstep ', ktberg,' unpacking berg ',INT(ibuffer_s%data(16,i)),' from south'507 WRITE(numicb,*) 'bergstep ',nktberg,' unpacking berg ',INT(ibuffer_s%data(16,i)),' from south' 508 508 CALL flush( numicb ) 509 509 ENDIF … … 514 514 IF( nn_verbose_level > 0 ) THEN 515 515 ! compare the number of icebergs on this processor from the start to the end 516 nbergs_end = count_bergs()517 i = ( nbergs_rcvd_from_n + nbergs_rcvd_from_s + nbergs_rcvd_from_e + nbergs_rcvd_from_w ) - &518 ( nbergs_to_send_n + nbergs_to_send_s + nbergs_to_send_e + nbergs_to_send_w )519 IF( nbergs_end-(nbergs_start+i) .NE. 0 ) THEN516 ibergs_end = count_bergs() 517 i = ( ibergs_rcvd_from_n + ibergs_rcvd_from_s + ibergs_rcvd_from_e + ibergs_rcvd_from_w ) - & 518 ( ibergs_to_send_n + ibergs_to_send_s + ibergs_to_send_e + ibergs_to_send_w ) 519 IF( ibergs_end-(ibergs_start+i) .NE. 0 ) THEN 520 520 WRITE( numicb,* ) 'send_bergs_to_other_pes: net change in number of icebergs' 521 WRITE( numicb,1000) 'send_bergs_to_other_pes: nbergs_end=', &522 nbergs_end,' on PE',narea523 WRITE( numicb,1000) 'send_bergs_to_other_pes: nbergs_start=', &524 nbergs_start,' on PE',narea521 WRITE( numicb,1000) 'send_bergs_to_other_pes: ibergs_end=', & 522 ibergs_end,' on PE',narea 523 WRITE( numicb,1000) 'send_bergs_to_other_pes: ibergs_start=', & 524 ibergs_start,' on PE',narea 525 525 WRITE( numicb,1000) 'send_bergs_to_other_pes: delta=', & 526 526 i,' on PE',narea 527 527 WRITE( numicb,1000) 'send_bergs_to_other_pes: error=', & 528 nbergs_end-(nbergs_start+i),' on PE',narea529 WRITE( numicb,1000) 'send_bergs_to_other_pes: nbergs_to_send_n=', &530 nbergs_to_send_n,' on PE',narea531 WRITE( numicb,1000) 'send_bergs_to_other_pes: nbergs_to_send_s=', &532 nbergs_to_send_s,' on PE',narea533 WRITE( numicb,1000) 'send_bergs_to_other_pes: nbergs_to_send_e=', &534 nbergs_to_send_e,' on PE',narea535 WRITE( numicb,1000) 'send_bergs_to_other_pes: nbergs_to_send_w=', &536 nbergs_to_send_w,' on PE',narea537 WRITE( numicb,1000) 'send_bergs_to_other_pes: nbergs_rcvd_from_n=', &538 nbergs_rcvd_from_n,' on PE',narea539 WRITE( numicb,1000) 'send_bergs_to_other_pes: nbergs_rcvd_from_s=', &540 nbergs_rcvd_from_s,' on PE',narea541 WRITE( numicb,1000) 'send_bergs_to_other_pes: nbergs_rcvd_from_e=', &542 nbergs_rcvd_from_e,' on PE',narea543 WRITE( numicb,1000) 'send_bergs_to_other_pes: nbergs_rcvd_from_w=', &544 nbergs_rcvd_from_w,' on PE',narea528 ibergs_end-(ibergs_start+i),' on PE',narea 529 WRITE( numicb,1000) 'send_bergs_to_other_pes: ibergs_to_send_n=', & 530 ibergs_to_send_n,' on PE',narea 531 WRITE( numicb,1000) 'send_bergs_to_other_pes: ibergs_to_send_s=', & 532 ibergs_to_send_s,' on PE',narea 533 WRITE( numicb,1000) 'send_bergs_to_other_pes: ibergs_to_send_e=', & 534 ibergs_to_send_e,' on PE',narea 535 WRITE( numicb,1000) 'send_bergs_to_other_pes: ibergs_to_send_w=', & 536 ibergs_to_send_w,' on PE',narea 537 WRITE( numicb,1000) 'send_bergs_to_other_pes: ibergs_rcvd_from_n=', & 538 ibergs_rcvd_from_n,' on PE',narea 539 WRITE( numicb,1000) 'send_bergs_to_other_pes: ibergs_rcvd_from_s=', & 540 ibergs_rcvd_from_s,' on PE',narea 541 WRITE( numicb,1000) 'send_bergs_to_other_pes: ibergs_rcvd_from_e=', & 542 ibergs_rcvd_from_e,' on PE',narea 543 WRITE( numicb,1000) 'send_bergs_to_other_pes: ibergs_rcvd_from_w=', & 544 ibergs_rcvd_from_w,' on PE',narea 545 545 1000 FORMAT(a,i5,a,i4) 546 546 CALL ctl_stop('send_bergs_to_other_pes: lost or gained an iceberg or two') … … 557 557 DO WHILE (ASSOCIATED(this)) 558 558 pt => this%current_point 559 i ne = INT( pt%xi + 0.5 )560 jne = INT( pt%yj + 0.5 )559 iine = INT( pt%xi + 0.5 ) 560 ijne = INT( pt%yj + 0.5 ) 561 561 ! CALL check_position(grd, this, 'exchange (bot)') 562 IF( i ne .LT. nimpp+icbdi-1 .OR. &563 i ne .GT. nimpp+icbei-1 .OR. &564 jne .LT. njmpp+icbdj-1 .OR. &565 jne .GT. njmpp+icbej-1) THEN562 IF( iine .LT. nimpp+nicbdi-1 .OR. & 563 iine .GT. nimpp+nicbei-1 .OR. & 564 ijne .LT. njmpp+nicbdj-1 .OR. & 565 ijne .GT. njmpp+nicbej-1) THEN 566 566 i = i + 1 567 WRITE(numicb,*) 'berg lost in halo: ', this%number(:),i ne,jne567 WRITE(numicb,*) 'berg lost in halo: ', this%number(:),iine,ijne 568 568 WRITE(numicb,*) ' ', nimpp, njmpp 569 WRITE(numicb,*) ' ', icbdi, icbei, icbdj,icbej569 WRITE(numicb,*) ' ', nicbdi, nicbei, nicbdj, nicbej 570 570 CALL flush( numicb ) 571 571 ENDIF … … 590 590 TYPE(iceberg) , POINTER :: tmpberg, this 591 591 TYPE(point) , POINTER :: pt 592 INTEGER :: nbergs_to_send593 INTEGER :: nbergs_to_rcv594 INTEGER :: jiglo, jjglo, jk, jn595 INTEGER :: jfldproc, kproc, kpts596 INTEGER :: i ne,jne597 REAL(wp), DIMENSION(2) :: sbergs,nbergs598 INTEGER :: ml_req1, ml_req2,ml_err599 INTEGER, DIMENSION(MPI_STATUS_SIZE) :: ml_stat592 INTEGER :: ibergs_to_send 593 INTEGER :: ibergs_to_rcv 594 INTEGER :: iiglo, ijglo, jk, jn 595 INTEGER :: ifldproc, iproc, ipts 596 INTEGER :: iine, ijne 597 REAL(wp), DIMENSION(2) :: zsbergs, znbergs 598 INTEGER :: iml_req1, iml_req2, iml_err 599 INTEGER, DIMENSION(MPI_STATUS_SIZE) :: iml_stat 600 600 601 601 ! set up indices of neighbouring processors 602 602 603 ! icbfldproc is a list of unique processor numbers that this processor603 ! nicbfldproc is a list of unique processor numbers that this processor 604 604 ! exchanges with (including itself), so we loop over this array; since 605 605 ! its of fixed size, the first -1 marks end of list of processors 606 606 ! 607 607 DO jn = 1, jpni 608 IF( icbfldproc(jn) == -1 ) EXIT609 jfldproc =icbfldproc(jn)610 nbergs_to_send = 0608 IF( nicbfldproc(jn) == -1 ) EXIT 609 ifldproc = nicbfldproc(jn) 610 ibergs_to_send = 0 611 611 612 612 ! Find number of bergs that need to be exchanged 613 ! Pick out exchanges with processor jfldproc614 ! if jfldproc is this processor then don't send613 ! Pick out exchanges with processor ifldproc 614 ! if ifldproc is this processor then don't send 615 615 ! 616 616 IF( ASSOCIATED(first_berg) ) THEN … … 618 618 DO WHILE (ASSOCIATED(this)) 619 619 pt => this%current_point 620 i ne = INT( pt%xi + 0.5 )621 jne = INT( pt%yj + 0.5 )622 kpts = icbfldpts (ine-nimpp+1)623 kproc = icbflddest(ine-nimpp+1)624 IF( jne .GT. njmpp+icbej-1 ) THEN625 IF( kproc == jfldproc ) THEN620 iine = INT( pt%xi + 0.5 ) 621 ijne = INT( pt%yj + 0.5 ) 622 ipts = nicbfldpts (iine-nimpp+1) 623 iproc = nicbflddest(iine-nimpp+1) 624 IF( ijne .GT. njmpp+nicbej-1 ) THEN 625 IF( iproc == ifldproc ) THEN 626 626 ! 627 627 ! moving across the cut line means both position and 628 628 ! velocity must change 629 jjglo = INT( kpts/icbpack )630 jiglo = kpts - icbpack*jjglo631 pt%xi = jiglo - ( pt%xi - REAL(ine,wp) )632 pt%yj = jjglo - ( pt%yj - REAL(jne,wp) )629 ijglo = INT( ipts/nicbpack ) 630 iiglo = ipts - nicbpack*ijglo 631 pt%xi = iiglo - ( pt%xi - REAL(iine,wp) ) 632 pt%yj = ijglo - ( pt%yj - REAL(ijne,wp) ) 633 633 pt%uvel = -1._wp * pt%uvel 634 634 pt%vvel = -1._wp * pt%vvel 635 635 ! 636 636 ! now remove berg from list and pack it into a buffer 637 IF( kproc /= narea ) THEN637 IF( iproc /= narea ) THEN 638 638 tmpberg => this 639 nbergs_to_send = nbergs_to_send + 1639 ibergs_to_send = ibergs_to_send + 1 640 640 IF( nn_verbose_level >= 4 ) THEN 641 WRITE(numicb,*) 'bergstep ', ktberg,' packing berg ',tmpberg%number(:),' for north fold'641 WRITE(numicb,*) 'bergstep ',nktberg,' packing berg ',tmpberg%number(:),' for north fold' 642 642 CALL flush( numicb ) 643 643 ENDIF 644 CALL pack_berg_into_buffer( tmpberg, obuffer_f, nbergs_to_send)644 CALL pack_berg_into_buffer( tmpberg, obuffer_f, ibergs_to_send) 645 645 CALL delete_iceberg_from_list(first_berg, tmpberg) 646 646 ENDIF … … 652 652 ENDIF 653 653 if( nn_verbose_level >= 3) then 654 write(numicb,*) 'bergstep ', ktberg,' send nfld: ', nbergs_to_send654 write(numicb,*) 'bergstep ',nktberg,' send nfld: ', ibergs_to_send 655 655 call flush(numicb) 656 656 endif … … 658 658 ! if we're in this processor, then we've done everything we need to 659 659 ! so go on to next element of loop 660 IF( jfldproc == narea ) CYCLE661 662 sbergs(1) = nbergs_to_send663 CALL mppsend( 21, sbergs(1), 1, jfldproc-1,ml_req1)664 CALL mpprecv( 21, nbergs(2), 1 )665 IF( l_isend ) CALL mpi_wait( ml_req1, ml_stat,ml_err )666 nbergs_to_rcv = INT(nbergs(2) )660 IF( ifldproc == narea ) CYCLE 661 662 zsbergs(1) = ibergs_to_send 663 CALL mppsend( 21, zsbergs(1), 1, ifldproc-1, iml_req1) 664 CALL mpprecv( 21, znbergs(2), 1 ) 665 IF( l_isend ) CALL mpi_wait( iml_req1, iml_stat, iml_err ) 666 ibergs_to_rcv = INT( znbergs(2) ) 667 667 668 668 ! send bergs 669 669 670 IF( nbergs_to_send > 0 ) &671 CALL mppsend( 12, obuffer_f%data, nbergs_to_send*buffer_width, jfldproc-1,ml_req2 )672 IF( nbergs_to_rcv > 0 ) THEN673 CALL increase_ibuffer(ibuffer_f, nbergs_to_rcv)674 CALL mpprecv( 12, ibuffer_f%data, nbergs_to_rcv*buffer_width )675 ENDIF 676 IF( nbergs_to_send > 0 .AND. l_isend ) CALL mpi_wait( ml_req2, ml_stat,ml_err )677 DO jk = 1, nbergs_to_rcv670 IF( ibergs_to_send > 0 ) & 671 CALL mppsend( 12, obuffer_f%data, ibergs_to_send*pp_buffer_width, ifldproc-1, iml_req2 ) 672 IF( ibergs_to_rcv > 0 ) THEN 673 CALL increase_ibuffer(ibuffer_f, ibergs_to_rcv) 674 CALL mpprecv( 12, ibuffer_f%data, ibergs_to_rcv*pp_buffer_width ) 675 ENDIF 676 IF( ibergs_to_send > 0 .AND. l_isend ) CALL mpi_wait( iml_req2, iml_stat, iml_err ) 677 DO jk = 1, ibergs_to_rcv 678 678 IF( nn_verbose_level >= 4 ) THEN 679 WRITE(numicb,*) 'bergstep ', ktberg,' unpacking berg ',INT(ibuffer_f%data(16,jk)),' from north fold'679 WRITE(numicb,*) 'bergstep ',nktberg,' unpacking berg ',INT(ibuffer_f%data(16,jk)),' from north fold' 680 680 CALL flush( numicb ) 681 681 ENDIF … … 704 704 !!------------------------------------------------------------------------- 705 705 706 SUBROUTINE dealloc_buffer( buff)706 SUBROUTINE dealloc_buffer(pbuff) 707 707 ! Arguments 708 TYPE(buffer), POINTER :: buff709 710 IF( ASSOCIATED( buff) ) THEN711 IF( ASSOCIATED( buff%data)) DEALLOCATE(buff%data)712 DEALLOCATE( buff)708 TYPE(buffer), POINTER :: pbuff 709 710 IF( ASSOCIATED(pbuff) ) THEN 711 IF( ASSOCIATED(pbuff%data)) DEALLOCATE(pbuff%data) 712 DEALLOCATE(pbuff) 713 713 ENDIF 714 714 END SUBROUTINE dealloc_buffer … … 716 716 !!------------------------------------------------------------------------- 717 717 718 SUBROUTINE pack_berg_into_buffer(berg, buff, n)718 SUBROUTINE pack_berg_into_buffer(berg, pbuff, kb) 719 719 ! Arguments 720 720 TYPE(iceberg), POINTER :: berg 721 TYPE(buffer) , POINTER :: buff722 INTEGER , INTENT(in) :: n721 TYPE(buffer) , POINTER :: pbuff 722 INTEGER , INTENT(in) :: kb 723 723 ! Local 724 724 INTEGER :: k 725 725 726 IF( .NOT. ASSOCIATED( buff) ) CALL increase_buffer( buff,delta_buf )727 IF( n .GT. buff%size ) CALL increase_buffer( buff,delta_buf )726 IF( .NOT. ASSOCIATED(pbuff) ) CALL increase_buffer( pbuff, pp_delta_buf ) 727 IF( kb .GT. pbuff%size ) CALL increase_buffer( pbuff, pp_delta_buf ) 728 728 729 729 !! pack points into buffer 730 730 731 buff%data( 1,n) = berg%current_point%lon732 buff%data( 2,n) = berg%current_point%lat733 buff%data( 3,n) = berg%current_point%uvel734 buff%data( 4,n) = berg%current_point%vvel735 buff%data( 5,n) = berg%current_point%xi736 buff%data( 6,n) = berg%current_point%yj737 buff%data( 7,n) = float(berg%current_point%year)738 buff%data( 8,n) = berg%current_point%day739 buff%data( 9,n) = berg%current_point%mass740 buff%data(10,n) = berg%current_point%thickness741 buff%data(11,n) = berg%current_point%width742 buff%data(12,n) = berg%current_point%length743 buff%data(13,n) = berg%current_point%mass_of_bits744 buff%data(14,n) = berg%current_point%heat_density745 746 buff%data(15,n) = berg%mass_scaling731 pbuff%data( 1,kb) = berg%current_point%lon 732 pbuff%data( 2,kb) = berg%current_point%lat 733 pbuff%data( 3,kb) = berg%current_point%uvel 734 pbuff%data( 4,kb) = berg%current_point%vvel 735 pbuff%data( 5,kb) = berg%current_point%xi 736 pbuff%data( 6,kb) = berg%current_point%yj 737 pbuff%data( 7,kb) = float(berg%current_point%year) 738 pbuff%data( 8,kb) = berg%current_point%day 739 pbuff%data( 9,kb) = berg%current_point%mass 740 pbuff%data(10,kb) = berg%current_point%thickness 741 pbuff%data(11,kb) = berg%current_point%width 742 pbuff%data(12,kb) = berg%current_point%length 743 pbuff%data(13,kb) = berg%current_point%mass_of_bits 744 pbuff%data(14,kb) = berg%current_point%heat_density 745 746 pbuff%data(15,kb) = berg%mass_scaling 747 747 DO k=1,nkounts 748 buff%data(15+k,n) = REAL( berg%number(k), wp )748 pbuff%data(15+k,kb) = REAL( berg%number(k), wp ) 749 749 END DO 750 750 … … 753 753 !!------------------------------------------------------------------------- 754 754 755 SUBROUTINE unpack_berg_from_buffer(first, buff, n)755 SUBROUTINE unpack_berg_from_buffer(first, pbuff, kb) 756 756 ! Arguments 757 757 TYPE(iceberg), POINTER :: first 758 TYPE(buffer) , POINTER :: buff759 INTEGER , INTENT(in) :: n758 TYPE(buffer) , POINTER :: pbuff 759 INTEGER , INTENT(in) :: kb 760 760 ! Local variables 761 LOGICAL :: lres762 761 TYPE(iceberg) :: currentberg 763 762 TYPE(point) :: pt 764 INTEGER ::k765 766 pt%lon = buff%data( 1,n)767 pt%lat = buff%data( 2,n)768 pt%uvel = buff%data( 3,n)769 pt%vvel = buff%data( 4,n)770 pt%xi = buff%data( 5,n)771 pt%yj = buff%data( 6,n)772 pt%year = INT( buff%data( 7,n) )773 pt%day = buff%data( 8,n)774 pt%mass = buff%data( 9,n)775 pt%thickness = buff%data(10,n)776 pt%width = buff%data(11,n)777 pt%length = buff%data(12,n)778 pt%mass_of_bits = buff%data(13,n)779 pt%heat_density = buff%data(14,n)780 781 currentberg%mass_scaling = buff%data(15,n)782 DO k=1,nkounts783 currentberg%number( k) = INT( buff%data(15+k,n) )763 INTEGER :: ik 764 765 pt%lon = pbuff%data( 1,kb) 766 pt%lat = pbuff%data( 2,kb) 767 pt%uvel = pbuff%data( 3,kb) 768 pt%vvel = pbuff%data( 4,kb) 769 pt%xi = pbuff%data( 5,kb) 770 pt%yj = pbuff%data( 6,kb) 771 pt%year = INT( pbuff%data( 7,kb) ) 772 pt%day = pbuff%data( 8,kb) 773 pt%mass = pbuff%data( 9,kb) 774 pt%thickness = pbuff%data(10,kb) 775 pt%width = pbuff%data(11,kb) 776 pt%length = pbuff%data(12,kb) 777 pt%mass_of_bits = pbuff%data(13,kb) 778 pt%heat_density = pbuff%data(14,kb) 779 780 currentberg%mass_scaling = pbuff%data(15,kb) 781 DO ik=1,nkounts 782 currentberg%number(ik) = INT( pbuff%data(15+ik,kb) ) 784 783 END DO 785 784 … … 790 789 !!------------------------------------------------------------------------- 791 790 792 SUBROUTINE increase_buffer(old, delta)791 SUBROUTINE increase_buffer(old,kdelta) 793 792 ! Arguments 794 793 TYPE(buffer), POINTER :: old 795 INTEGER , INTENT(in) :: delta794 INTEGER , INTENT(in) :: kdelta 796 795 ! Local variables 797 796 TYPE(buffer), POINTER :: new 798 INTEGER :: new_size797 INTEGER :: inew_size 799 798 800 799 IF( .NOT. ASSOCIATED(old) ) THEN 801 new_size =delta800 inew_size = kdelta 802 801 ELSE 803 new_size = old%size +delta802 inew_size = old%size + kdelta 804 803 ENDIF 805 804 ALLOCATE( new ) 806 ALLOCATE( new%data( buffer_width,new_size) )807 new%size = new_size805 ALLOCATE( new%data( pp_buffer_width, inew_size) ) 806 new%size = inew_size 808 807 IF( ASSOCIATED(old) ) THEN 809 808 new%data(:,1:old%size) = old%data(:,1:old%size) … … 817 816 !!------------------------------------------------------------------------- 818 817 819 SUBROUTINE increase_ibuffer(old, delta)818 SUBROUTINE increase_ibuffer(old,kdelta) 820 819 ! Arguments 821 820 TYPE(buffer), POINTER :: old 822 INTEGER , INTENT(in) :: delta821 INTEGER , INTENT(in) :: kdelta 823 822 ! Local variables 824 823 TYPE(buffer), POINTER :: new 825 INTEGER :: new_size,old_size824 INTEGER :: inew_size, iold_size 826 825 827 826 IF( .NOT. ASSOCIATED(old) ) THEN 828 new_size = delta +delta_buf829 old_size = 0827 inew_size = kdelta + pp_delta_buf 828 iold_size = 0 830 829 ELSE 831 old_size = old%size832 IF( delta .LT. old%size ) THEN833 new_size = old%size +delta830 iold_size = old%size 831 IF( kdelta .LT. old%size ) THEN 832 inew_size = old%size + kdelta 834 833 ELSE 835 new_size = delta +delta_buf836 ENDIF 837 ENDIF 838 839 IF( old_size .NE.new_size ) THEN834 inew_size = kdelta + pp_delta_buf 835 ENDIF 836 ENDIF 837 838 IF( iold_size .NE. inew_size ) THEN 840 839 ALLOCATE( new ) 841 ALLOCATE( new%data( buffer_width,new_size) )842 new%size = new_size840 ALLOCATE( new%data( pp_buffer_width, inew_size) ) 841 new%size = inew_size 843 842 IF( ASSOCIATED(old) ) THEN 844 843 new%data(:,1:old%size) = old%data(:,1:old%size) … … 847 846 ENDIF 848 847 old => new 849 !WRITE( numicb,*) 'increase_ibuffer',narea,' increased to', new_size848 !WRITE( numicb,*) 'increase_ibuffer',narea,' increased to',inew_size 850 849 ENDIF 851 850 -
branches/2012/dev_r3337_NOCS10_ICB/NEMOGCM/NEMO/OPA_SRC/ICB/icbrst.F90
r3339 r3359 29 29 30 30 REAL(wp), DIMENSION(:,:,:), PRIVATE, POINTER :: griddata => NULL() ! need 2d array to read in with 31 INTEGER, DIMENSION(3), PRIVATE :: strt3, lngth3 32 INTEGER, PRIVATE :: lonid, latid, xid, yid, uvelid, vvelid 33 INTEGER, PRIVATE :: massid, thicknessid, widthid, lengthid 34 INTEGER, PRIVATE :: yearid, dayid 35 INTEGER, PRIVATE :: scaling_id, mass_of_bits_id, heat_density_id, numberid 36 INTEGER, PRIVATE :: siceid, sheatid, calvid, calvhid, kountid 37 38 PUBLIC icebergs_read_restart ! routine called in xxx.F90 module 39 PUBLIC icebergs_write_restart ! routine called in xxx.F90 module 31 INTEGER, DIMENSION(3), PRIVATE :: nstrt3, nlngth3 32 INTEGER, PRIVATE :: nlonid, nlatid, nxid, nyid, nuvelid, nvvelid 33 INTEGER, PRIVATE :: nmassid, nthicknessid, nwidthid, nlengthid 34 INTEGER, PRIVATE :: nyearid, ndayid 35 INTEGER, PRIVATE :: nscaling_id, nmass_of_bits_id, nheat_density_id, numberid 36 INTEGER, PRIVATE :: nsiceid, nsheatid, ncalvid, ncalvhid, nkountid 37 INTEGER, PRIVATE :: nret, ncid, nc_dim 38 39 PUBLIC icebergs_read_restart ! routine called in icbini.F90 module 40 PUBLIC icebergs_write_restart ! routine called in icbrun.F90 module 40 41 41 42 CONTAINS … … 47 48 48 49 ! Local variables 49 INTEGER :: ndim, nvar, natt50 INTEGER :: k, iret, ncid, unlim_dim, x_dim, y_dim, c_dim, nbergs_in_file51 INTEGER :: mclass52 INTEGER, DIMENSION(1) :: strt,lngth, idata53 INTEGER, DIMENSION(2) :: strt2,lngth254 INTEGER, DIMENSION(nkounts) :: idata255 REAL(wp), DIMENSION(1) :: data ! need 1d array to read in with50 INTEGER :: idim, ivar, iatt 51 INTEGER :: jn, iunlim_dim, ibergs_in_file 52 INTEGER :: iclass 53 INTEGER, DIMENSION(1) :: istrt, ilngth, idata 54 INTEGER, DIMENSION(2) :: istrt2, ilngth2 55 INTEGER, DIMENSION(nkounts) :: idata2 56 REAL(wp), DIMENSION(1) :: zdata ! need 1d array to read in with 56 57 ! start and count arrays 57 LOGICAL :: lres, found_restart,multiPErestart=.FALSE.58 CHARACTER(len=80) :: filename59 CHARACTER(len=NF90_MAX_NAME) :: dname60 TYPE(iceberg) :: localberg ! NOT a pointer but an actual local variable61 TYPE(point) :: localpt ! NOT a pointer but an actual local variable58 LOGICAL :: ll_found_restart, ll_multiPErestart=.FALSE. 59 CHARACTER(len=80) :: cl_filename 60 CHARACTER(len=NF90_MAX_NAME) :: cl_dname 61 TYPE(iceberg) :: localberg ! NOT a pointer but an actual local variable 62 TYPE(point) :: localpt ! NOT a pointer but an actual local variable 62 63 63 64 IF( .NOT. ASSOCIATED(griddata) ) ALLOCATE( griddata(jpi,jpj,1) ) 64 65 65 66 ! Find a restart file 66 multiPErestart=.FALSE.67 ll_multiPErestart=.FALSE. 67 68 DO 68 filename = ' '69 filename = 'restart_icebergs.nc'70 INQUIRE( file=TRIM(filename),exist=found_restart)71 IF ( found_restart) EXIT72 filename = ' '73 WRITE( filename,'("restart_icebergs_",I4.4,".nc")') narea-174 INQUIRE( file=TRIM(filename),exist=found_restart)75 IF ( found_restart) THEN76 multiPErestart=.TRUE.69 cl_filename = ' ' 70 cl_filename = 'restart_icebergs.nc' 71 INQUIRE( file=TRIM(cl_filename), exist=ll_found_restart ) 72 IF ( ll_found_restart ) EXIT 73 cl_filename = ' ' 74 WRITE( cl_filename, '("restart_icebergs_",I4.4,".nc")' ) narea-1 75 INQUIRE( file=TRIM(cl_filename), exist=ll_found_restart ) 76 IF ( ll_found_restart ) THEN 77 ll_multiPErestart = .TRUE. 77 78 EXIT 78 79 ENDIF 79 80 IF (nn_verbose_level >= 0 .AND. lwp) & 80 WRITE( numout,'(a)') 'read_restart_bergs: no restart file found'81 multiPErestart=.TRUE. ! force sanity checking in a multi-PE modeif no file found on this PE81 WRITE( numout, '(a)' ) 'read_restart_bergs: no restart file found' 82 ll_multiPErestart = .TRUE. ! force checking in a MPP if no file found on this PE 82 83 EXIT 83 84 ENDDO 84 85 85 IF ( .NOT. found_restart) THEN! only do the following if a file was found86 IF ( .NOT. ll_found_restart) THEN ! only do the following if a file was found 86 87 CALL ctl_stop('icebergs: no restart file found') 87 88 ENDIF 88 89 89 90 IF (nn_verbose_level >= 0 .AND. lwp) & 90 WRITE(numout,'(2a)') 'icebergs, read_restart_bergs: found restart file = ',TRIM( filename)91 92 iret = NF90_OPEN(TRIM(filename), NF90_NOWRITE, ncid)93 IF ( iret .ne. NF90_NOERR) CALL ctl_stop('icebergs, read_restart_bergs: nf_open failed')94 95 iret = nf90_inquire(ncid, ndim, nvar, natt,unlim_dim)96 IF ( iret .ne. NF90_NOERR) CALL ctl_stop('icebergs, read_restart_bergs: nf_inquire failed')97 98 IF( unlim_dim .NE. -1) THEN99 100 iret = nf90_inquire_dimension(ncid, unlim_dim, dname, nbergs_in_file)101 IF ( iret .ne. NF90_NOERR) CALL ctl_stop('icebergs, read_restart_bergs: nf_inq_dimlen failed')102 103 iret = NF90_INQ_VARID(ncid, 'number', numberid)104 iret = NF90_INQ_VARID(ncid, 'mass_scaling',scaling_id)105 iret = NF90_INQ_VARID(ncid, 'xi',xid)106 iret = NF90_INQ_VARID(ncid, 'yj',yid)107 iret = NF90_INQ_VARID(ncid, 'lon',lonid)108 iret = NF90_INQ_VARID(ncid, 'lat',latid)109 iret = NF90_INQ_VARID(ncid, 'uvel',uvelid)110 iret = NF90_INQ_VARID(ncid, 'vvel',vvelid)111 iret = NF90_INQ_VARID(ncid, 'mass',massid)112 iret = NF90_INQ_VARID(ncid, 'thickness',thicknessid)113 iret = NF90_INQ_VARID(ncid, 'width',widthid)114 iret = NF90_INQ_VARID(ncid, 'length',lengthid)115 iret = NF90_INQ_VARID(ncid, 'year',yearid)116 iret = NF90_INQ_VARID(ncid, 'day',dayid)117 iret = NF90_INQ_VARID(ncid, 'mass_of_bits',mass_of_bits_id)118 iret = NF90_INQ_VARID(ncid, 'heat_density',heat_density_id)119 120 lngth(1) = 1121 strt2(1) = 1122 lngth2(1) = nkounts123 lngth2(2) = 1124 DO k=1, nbergs_in_file125 126 strt(1) = k127 strt2(2) = k128 129 iret = NF90_GET_VAR(ncid, numberid, idata2, strt2,lngth2 )91 WRITE(numout,'(2a)') 'icebergs, read_restart_bergs: found restart file = ',TRIM(cl_filename) 92 93 nret = NF90_OPEN(TRIM(cl_filename), NF90_NOWRITE, ncid) 94 IF (nret .ne. NF90_NOERR) CALL ctl_stop('icebergs, read_restart_bergs: nf_open failed') 95 96 nret = nf90_inquire(ncid, idim, ivar, iatt, iunlim_dim) 97 IF (nret .ne. NF90_NOERR) CALL ctl_stop('icebergs, read_restart_bergs: nf_inquire failed') 98 99 IF( iunlim_dim .NE. -1) THEN 100 101 nret = nf90_inquire_dimension(ncid, iunlim_dim, cl_dname, ibergs_in_file) 102 IF (nret .ne. NF90_NOERR) CALL ctl_stop('icebergs, read_restart_bergs: nf_inq_dimlen failed') 103 104 nret = NF90_INQ_VARID(ncid, 'number', numberid) 105 nret = NF90_INQ_VARID(ncid, 'mass_scaling', nscaling_id) 106 nret = NF90_INQ_VARID(ncid, 'xi', nxid) 107 nret = NF90_INQ_VARID(ncid, 'yj', nyid) 108 nret = NF90_INQ_VARID(ncid, 'lon', nlonid) 109 nret = NF90_INQ_VARID(ncid, 'lat', nlatid) 110 nret = NF90_INQ_VARID(ncid, 'uvel', nuvelid) 111 nret = NF90_INQ_VARID(ncid, 'vvel', nvvelid) 112 nret = NF90_INQ_VARID(ncid, 'mass', nmassid) 113 nret = NF90_INQ_VARID(ncid, 'thickness', nthicknessid) 114 nret = NF90_INQ_VARID(ncid, 'width', nwidthid) 115 nret = NF90_INQ_VARID(ncid, 'length', nlengthid) 116 nret = NF90_INQ_VARID(ncid, 'year', nyearid) 117 nret = NF90_INQ_VARID(ncid, 'day', ndayid) 118 nret = NF90_INQ_VARID(ncid, 'mass_of_bits', nmass_of_bits_id) 119 nret = NF90_INQ_VARID(ncid, 'heat_density', nheat_density_id) 120 121 ilngth(1) = 1 122 istrt2(1) = 1 123 ilngth2(1) = nkounts 124 ilngth2(2) = 1 125 DO jn=1, ibergs_in_file 126 127 istrt(1) = jn 128 istrt2(2) = jn 129 130 nret = NF90_GET_VAR(ncid, numberid, idata2, istrt2, ilngth2 ) 130 131 localberg%number(:) = idata2(:) 131 132 132 iret = NF90_GET_VAR(ncid, scaling_id, data, strt,lngth )133 localberg%mass_scaling = data(1)134 135 iret = NF90_GET_VAR(ncid, lonid, data, strt,lngth)136 localpt%lon = data(1)137 iret = NF90_GET_VAR(ncid, latid, data, strt,lngth)138 localpt%lat = data(1)133 nret = NF90_GET_VAR(ncid, nscaling_id, zdata, istrt, ilngth ) 134 localberg%mass_scaling = zdata(1) 135 136 nret = NF90_GET_VAR(ncid, nlonid, zdata, istrt, ilngth) 137 localpt%lon = zdata(1) 138 nret = NF90_GET_VAR(ncid, nlatid, zdata, istrt, ilngth) 139 localpt%lat = zdata(1) 139 140 IF (nn_verbose_level >= 2 .AND. lwp) THEN 140 WRITE(numout,'(a,i5,a,2f10.4,a,i5)') 'icebergs, read_restart_bergs: berg ', k,' is at ', &141 WRITE(numout,'(a,i5,a,2f10.4,a,i5)') 'icebergs, read_restart_bergs: berg ',jn,' is at ', & 141 142 localpt%lon,localpt%lat,' on PE ',narea-1 142 143 ENDIF 143 iret = NF90_GET_VAR(ncid, xid, data, strt,lngth)144 localpt%xi = data(1)145 iret = NF90_GET_VAR(ncid, yid, data, strt,lngth)146 localpt%yj = data(1)147 iret = NF90_GET_VAR(ncid, uvelid, data, strt,lngth )148 localpt%uvel = data(1)149 iret = NF90_GET_VAR(ncid, vvelid, data, strt,lngth )150 localpt%vvel = data(1)151 iret = NF90_GET_VAR(ncid, massid, data, strt,lngth )152 localpt%mass = data(1)153 iret = NF90_GET_VAR(ncid, thicknessid, data, strt,lngth )154 localpt%thickness = data(1)155 iret = NF90_GET_VAR(ncid, widthid, data, strt,lngth )156 localpt%width = data(1)157 iret = NF90_GET_VAR(ncid, lengthid, data, strt,lngth )158 localpt%length = data(1)159 iret = NF90_GET_VAR(ncid, yearid, idata, strt,lngth )144 nret = NF90_GET_VAR(ncid, nxid, zdata, istrt, ilngth) 145 localpt%xi = zdata(1) 146 nret = NF90_GET_VAR(ncid, nyid, zdata, istrt, ilngth) 147 localpt%yj = zdata(1) 148 nret = NF90_GET_VAR(ncid, nuvelid, zdata, istrt, ilngth ) 149 localpt%uvel = zdata(1) 150 nret = NF90_GET_VAR(ncid, nvvelid, zdata, istrt, ilngth ) 151 localpt%vvel = zdata(1) 152 nret = NF90_GET_VAR(ncid, nmassid, zdata, istrt, ilngth ) 153 localpt%mass = zdata(1) 154 nret = NF90_GET_VAR(ncid, nthicknessid, zdata, istrt, ilngth ) 155 localpt%thickness = zdata(1) 156 nret = NF90_GET_VAR(ncid, nwidthid, zdata, istrt, ilngth ) 157 localpt%width = zdata(1) 158 nret = NF90_GET_VAR(ncid, nlengthid, zdata, istrt, ilngth ) 159 localpt%length = zdata(1) 160 nret = NF90_GET_VAR(ncid, nyearid, idata, istrt, ilngth ) 160 161 localpt%year = idata(1) 161 iret = NF90_GET_VAR(ncid, dayid, data, strt,lngth )162 localpt%day = data(1)163 iret = NF90_GET_VAR(ncid, mass_of_bits_id, data, strt,lngth )164 localpt%mass_of_bits = data(1)165 iret = NF90_GET_VAR(ncid, heat_density_id, data, strt,lngth )166 localpt%heat_density = data(1)162 nret = NF90_GET_VAR(ncid, ndayid, zdata, istrt, ilngth ) 163 localpt%day = zdata(1) 164 nret = NF90_GET_VAR(ncid, nmass_of_bits_id, zdata, istrt, ilngth ) 165 localpt%mass_of_bits = zdata(1) 166 nret = NF90_GET_VAR(ncid, nheat_density_id, zdata, istrt, ilngth ) 167 localpt%heat_density = zdata(1) 167 168 ! 168 169 CALL add_new_berg_to_list( localberg, localpt ) … … 171 172 ENDIF 172 173 173 iret = NF90_INQ_DIMID( ncid, 'c',c_dim )174 IF ( iret .ne. NF90_NOERR) CALL ctl_stop('icebergs, read_restart: nf_inq_dimid c failed')175 176 iret = NF90_INQUIRE_DIMENSION( ncid, c_dim, dname, mclass )177 IF ( iret .ne. NF90_NOERR) CALL ctl_stop('icebergs, read_restart: nf_inquire_dimension failed')178 179 iret = NF90_INQ_VARID(ncid, 'kount' ,kountid)180 iret = NF90_INQ_VARID(ncid, 'calving' ,calvid)181 iret = NF90_INQ_VARID(ncid, 'calving_hflx',calvhid)182 iret = NF90_INQ_VARID(ncid, 'stored_ice' ,siceid)183 iret = NF90_INQ_VARID(ncid, 'stored_heat' ,sheatid)184 185 strt3(1) = 1186 strt3(2) = 1187 lngth3(1) = jpi188 lngth3(2) = jpj189 lngth3(3) = 1190 191 DO k=1,mclass192 193 strt3(3) = k194 195 iret = NF90_GET_VAR( ncid, siceid , griddata, strt3,lngth3 )196 berg_grid%stored_ice(:,:, k) = griddata(:,:,1)174 nret = NF90_INQ_DIMID( ncid, 'c', nc_dim ) 175 IF (nret .ne. NF90_NOERR) CALL ctl_stop('icebergs, read_restart: nf_inq_dimid c failed') 176 177 nret = NF90_INQUIRE_DIMENSION( ncid, nc_dim, cl_dname, iclass ) 178 IF (nret .ne. NF90_NOERR) CALL ctl_stop('icebergs, read_restart: nf_inquire_dimension failed') 179 180 nret = NF90_INQ_VARID(ncid, 'kount' , nkountid) 181 nret = NF90_INQ_VARID(ncid, 'calving' , ncalvid) 182 nret = NF90_INQ_VARID(ncid, 'calving_hflx', ncalvhid) 183 nret = NF90_INQ_VARID(ncid, 'stored_ice' , nsiceid) 184 nret = NF90_INQ_VARID(ncid, 'stored_heat' , nsheatid) 185 186 nstrt3(1) = 1 187 nstrt3(2) = 1 188 nlngth3(1) = jpi 189 nlngth3(2) = jpj 190 nlngth3(3) = 1 191 192 DO jn=1,iclass 193 194 nstrt3(3) = jn 195 196 nret = NF90_GET_VAR( ncid, nsiceid , griddata, nstrt3, nlngth3 ) 197 berg_grid%stored_ice(:,:,jn) = griddata(:,:,1) 197 198 198 199 ENDDO 199 200 200 iret = NF90_GET_VAR( ncid, calvid , p_calving(:,:) )201 iret = NF90_GET_VAR( ncid, calvhid, p_calving_hflx(:,:) )202 iret = NF90_GET_VAR( ncid,sheatid, berg_grid%stored_heat(:,:) )203 iret = NF90_GET_VAR( ncid,kountid, idata2(:) )204 kount_bergs(:) = idata2(:)201 nret = NF90_GET_VAR( ncid, ncalvid , src_calving(:,:) ) 202 nret = NF90_GET_VAR( ncid, ncalvhid, src_calving_hflx(:,:) ) 203 nret = NF90_GET_VAR( ncid, nsheatid, berg_grid%stored_heat(:,:) ) 204 nret = NF90_GET_VAR( ncid, nkountid, idata2(:) ) 205 num_bergs(:) = idata2(:) 205 206 206 207 ! Finish up 207 iret = NF90_CLOSE(ncid)208 IF ( iret .ne. NF90_NOERR) CALL ctl_stop('icebergs, read_restart: nf_close failed')208 nret = NF90_CLOSE(ncid) 209 IF (nret .ne. NF90_NOERR) CALL ctl_stop('icebergs, read_restart: nf_close failed') 209 210 210 211 ! Sanity check 211 k= count_bergs()212 jn = count_bergs() 212 213 IF (nn_verbose_level >= 0) & 213 WRITE(numout,'(2(a,i5))') 'icebergs, read_restart_bergs: # bergs =', k,' on PE',narea-1214 WRITE(numout,'(2(a,i5))') 'icebergs, read_restart_bergs: # bergs =',jn,' on PE',narea-1 214 215 IF( lk_mpp ) THEN 215 IF ( multiPErestart) CALL mpp_sum(nbergs_in_file) ! In case PE 0 didn't open a file216 CALL mpp_sum( k)216 IF (ll_multiPErestart) CALL mpp_sum(ibergs_in_file) ! In case PE 0 didn't open a file 217 CALL mpp_sum(jn) 217 218 ENDIF 218 219 IF ( lwp ) THEN 219 WRITE(numout,'(a,i5,a,i5,a)') 'icebergs, read_restart_bergs: there were', nbergs_in_file, &220 ' bergs in the restart file and', k,' bergs have been read'220 WRITE(numout,'(a,i5,a,i5,a)') 'icebergs, read_restart_bergs: there were',ibergs_in_file, & 221 ' bergs in the restart file and', jn,' bergs have been read' 221 222 ENDIF 222 223 … … 232 233 INTEGER, INTENT( in ) :: kt 233 234 ! Local variables 234 INTEGER :: iret, ncid, k_dim, n_dim, x_dim, y_dim, c_dim, k 235 CHARACTER(len=80) :: filename 235 INTEGER, :: ix_dim, iy_dim, ik_dim, in_dim 236 INTEGER :: jn 237 CHARACTER(len=80) :: cl_filename 236 238 TYPE(iceberg), POINTER :: this 237 239 TYPE(point) , POINTER :: pt … … 239 241 IF( .NOT. ASSOCIATED(griddata) ) ALLOCATE( griddata(jpi,jpj,1) ) 240 242 241 WRITE( filename,'("icebergs_",I8.8,"_restart_",I4.4,".nc")') kt, narea-1242 IF (nn_verbose_level >= 0) WRITE(numout,'(2a)') 'icebergs, write_restart: creating ',TRIM( filename)243 244 iret = NF90_CREATE(TRIM(filename), NF90_CLOBBER, ncid)245 IF ( iret .ne. NF90_NOERR) CALL ctl_stop('icebergs, write_restart: nf_create failed')243 WRITE(cl_filename,'("icebergs_",I8.8,"_restart_",I4.4,".nc")') kt, narea-1 244 IF (nn_verbose_level >= 0) WRITE(numout,'(2a)') 'icebergs, write_restart: creating ',TRIM(cl_filename) 245 246 nret = NF90_CREATE(TRIM(cl_filename), NF90_CLOBBER, ncid) 247 IF (nret .ne. NF90_NOERR) CALL ctl_stop('icebergs, write_restart: nf_create failed') 246 248 247 249 ! Dimensions 248 iret = NF90_DEF_DIM(ncid, 'x', jpi,x_dim)249 IF ( iret .ne. NF90_NOERR) CALL ctl_stop('icebergs, write_restart: nf_def_dim x failed')250 251 iret = NF90_DEF_DIM(ncid, 'y', jpj,y_dim)252 IF ( iret .ne. NF90_NOERR) CALL ctl_stop('icebergs, write_restart: nf_def_dim y failed')253 254 iret = NF90_DEF_DIM(ncid, 'c', nclasses,c_dim)255 IF ( iret .ne. NF90_NOERR) CALL ctl_stop('icebergs, write_restart: nf_def_dim c failed')256 257 iret = NF90_DEF_DIM(ncid, 'k', nkounts,k_dim)258 IF ( iret .ne. NF90_NOERR) CALL ctl_stop('icebergs, write_restart: nf_def_dim k failed')250 nret = NF90_DEF_DIM(ncid, 'x', jpi, ix_dim) 251 IF (nret .ne. NF90_NOERR) CALL ctl_stop('icebergs, write_restart: nf_def_dim x failed') 252 253 nret = NF90_DEF_DIM(ncid, 'y', jpj, iy_dim) 254 IF (nret .ne. NF90_NOERR) CALL ctl_stop('icebergs, write_restart: nf_def_dim y failed') 255 256 nret = NF90_DEF_DIM(ncid, 'c', nclasses, nc_dim) 257 IF (nret .ne. NF90_NOERR) CALL ctl_stop('icebergs, write_restart: nf_def_dim c failed') 258 259 nret = NF90_DEF_DIM(ncid, 'k', nkounts, ik_dim) 260 IF (nret .ne. NF90_NOERR) CALL ctl_stop('icebergs, write_restart: nf_def_dim k failed') 259 261 260 262 IF (associated(first_berg)) then 261 iret = NF90_DEF_DIM(ncid, 'n', NF90_UNLIMITED,n_dim)262 IF ( iret .ne. NF90_NOERR) CALL ctl_stop('icebergs, write_restart: nf_def_dim n failed')263 nret = NF90_DEF_DIM(ncid, 'n', NF90_UNLIMITED, in_dim) 264 IF (nret .ne. NF90_NOERR) CALL ctl_stop('icebergs, write_restart: nf_def_dim n failed') 263 265 ENDIF 264 266 265 267 ! Variables 266 iret = NF90_DEF_VAR(ncid, 'kount' , NF90_INT , (/ k_dim /),kountid)267 iret = NF90_DEF_VAR(ncid, 'calving' , NF90_DOUBLE, (/ x_dim, y_dim /),calvid)268 iret = NF90_DEF_VAR(ncid, 'calving_hflx', NF90_DOUBLE, (/ x_dim, y_dim /),calvhid)269 iret = NF90_DEF_VAR(ncid, 'stored_ice' , NF90_DOUBLE, (/ x_dim, y_dim, c_dim /),siceid)270 iret = NF90_DEF_VAR(ncid, 'stored_heat' , NF90_DOUBLE, (/ x_dim, y_dim /),sheatid)268 nret = NF90_DEF_VAR(ncid, 'kount' , NF90_INT , (/ ik_dim /), nkountid) 269 nret = NF90_DEF_VAR(ncid, 'calving' , NF90_DOUBLE, (/ ix_dim, iy_dim /), ncalvid) 270 nret = NF90_DEF_VAR(ncid, 'calving_hflx', NF90_DOUBLE, (/ ix_dim, iy_dim /), ncalvhid) 271 nret = NF90_DEF_VAR(ncid, 'stored_ice' , NF90_DOUBLE, (/ ix_dim, iy_dim, nc_dim /), nsiceid) 272 nret = NF90_DEF_VAR(ncid, 'stored_heat' , NF90_DOUBLE, (/ ix_dim, iy_dim /), nsheatid) 271 273 272 274 ! Attributes 273 iret = NF90_PUT_ATT(ncid,calvid , 'long_name', 'iceberg calving')274 iret = NF90_PUT_ATT(ncid,calvid , 'units', 'some')275 iret = NF90_PUT_ATT(ncid,calvhid, 'long_name', 'heat flux associated with iceberg calving')276 iret = NF90_PUT_ATT(ncid,calvhid, 'units', 'some')277 iret = NF90_PUT_ATT(ncid,siceid , 'long_name', 'stored ice used to calve icebergs')278 iret = NF90_PUT_ATT(ncid,siceid , 'units', 'kg/s')279 iret = NF90_PUT_ATT(ncid,sheatid, 'long_name', 'heat in stored ice used to calve icebergs')280 iret = NF90_PUT_ATT(ncid,sheatid, 'units', 'J/kg/s')275 nret = NF90_PUT_ATT(ncid, ncalvid , 'long_name', 'iceberg calving') 276 nret = NF90_PUT_ATT(ncid, ncalvid , 'units', 'some') 277 nret = NF90_PUT_ATT(ncid, ncalvhid, 'long_name', 'heat flux associated with iceberg calving') 278 nret = NF90_PUT_ATT(ncid, ncalvhid, 'units', 'some') 279 nret = NF90_PUT_ATT(ncid, nsiceid , 'long_name', 'stored ice used to calve icebergs') 280 nret = NF90_PUT_ATT(ncid, nsiceid , 'units', 'kg/s') 281 nret = NF90_PUT_ATT(ncid, nsheatid, 'long_name', 'heat in stored ice used to calve icebergs') 282 nret = NF90_PUT_ATT(ncid, nsheatid, 'units', 'J/kg/s') 281 283 282 284 IF ( ASSOCIATED(first_berg) ) THEN … … 285 287 286 288 ! Variables 287 iret = NF90_DEF_VAR(ncid, 'lon', NF90_DOUBLE, n_dim,lonid)288 iret = NF90_DEF_VAR(ncid, 'lat', NF90_DOUBLE, n_dim,latid)289 iret = NF90_DEF_VAR(ncid, 'xi', NF90_DOUBLE, n_dim,xid)290 iret = NF90_DEF_VAR(ncid, 'yj', NF90_DOUBLE, n_dim,yid)291 iret = NF90_DEF_VAR(ncid, 'uvel', NF90_DOUBLE, n_dim,uvelid)292 iret = NF90_DEF_VAR(ncid, 'vvel', NF90_DOUBLE, n_dim,vvelid)293 iret = NF90_DEF_VAR(ncid, 'mass', NF90_DOUBLE, n_dim,massid)294 iret = NF90_DEF_VAR(ncid, 'thickness', NF90_DOUBLE, n_dim,thicknessid)295 iret = NF90_DEF_VAR(ncid, 'width', NF90_DOUBLE, n_dim,widthid)296 iret = NF90_DEF_VAR(ncid, 'length', NF90_DOUBLE, n_dim,lengthid)297 iret = NF90_DEF_VAR(ncid, 'number', NF90_INT, (/k_dim,n_dim/), numberid)298 iret = NF90_DEF_VAR(ncid, 'year', NF90_INT, n_dim,yearid)299 iret = NF90_DEF_VAR(ncid, 'day', NF90_DOUBLE, n_dim,dayid)300 iret = NF90_DEF_VAR(ncid, 'mass_scaling', NF90_DOUBLE, n_dim,scaling_id)301 iret = NF90_DEF_VAR(ncid, 'mass_of_bits', NF90_DOUBLE, n_dim,mass_of_bits_id)302 iret = NF90_DEF_VAR(ncid, 'heat_density', NF90_DOUBLE, n_dim,heat_density_id)289 nret = NF90_DEF_VAR(ncid, 'lon', NF90_DOUBLE, in_dim, nlonid) 290 nret = NF90_DEF_VAR(ncid, 'lat', NF90_DOUBLE, in_dim, nlatid) 291 nret = NF90_DEF_VAR(ncid, 'xi', NF90_DOUBLE, in_dim, nxid) 292 nret = NF90_DEF_VAR(ncid, 'yj', NF90_DOUBLE, in_dim, nyid) 293 nret = NF90_DEF_VAR(ncid, 'uvel', NF90_DOUBLE, in_dim, nuvelid) 294 nret = NF90_DEF_VAR(ncid, 'vvel', NF90_DOUBLE, in_dim, nvvelid) 295 nret = NF90_DEF_VAR(ncid, 'mass', NF90_DOUBLE, in_dim, nmassid) 296 nret = NF90_DEF_VAR(ncid, 'thickness', NF90_DOUBLE, in_dim, nthicknessid) 297 nret = NF90_DEF_VAR(ncid, 'width', NF90_DOUBLE, in_dim, nwidthid) 298 nret = NF90_DEF_VAR(ncid, 'length', NF90_DOUBLE, in_dim, nlengthid) 299 nret = NF90_DEF_VAR(ncid, 'number', NF90_INT, (/ik_dim,in_dim/), numberid) 300 nret = NF90_DEF_VAR(ncid, 'year', NF90_INT, in_dim, nyearid) 301 nret = NF90_DEF_VAR(ncid, 'day', NF90_DOUBLE, in_dim, ndayid) 302 nret = NF90_DEF_VAR(ncid, 'mass_scaling', NF90_DOUBLE, in_dim, nscaling_id) 303 nret = NF90_DEF_VAR(ncid, 'mass_of_bits', NF90_DOUBLE, in_dim, nmass_of_bits_id) 304 nret = NF90_DEF_VAR(ncid, 'heat_density', NF90_DOUBLE, in_dim, nheat_density_id) 303 305 304 306 ! Attributes 305 iret = NF90_PUT_ATT(ncid,lonid, 'long_name', 'longitude')306 iret = NF90_PUT_ATT(ncid,lonid, 'units', 'degrees_E')307 iret = NF90_PUT_ATT(ncid,latid, 'long_name', 'latitude')308 iret = NF90_PUT_ATT(ncid,latid, 'units', 'degrees_N')309 iret = NF90_PUT_ATT(ncid,xid, 'long_name', 'x grid box position')310 iret = NF90_PUT_ATT(ncid,xid, 'units', 'fractional')311 iret = NF90_PUT_ATT(ncid,yid, 'long_name', 'y grid box position')312 iret = NF90_PUT_ATT(ncid,yid, 'units', 'fractional')313 iret = NF90_PUT_ATT(ncid,uvelid, 'long_name', 'zonal velocity')314 iret = NF90_PUT_ATT(ncid,uvelid, 'units', 'm/s')315 iret = NF90_PUT_ATT(ncid,vvelid, 'long_name', 'meridional velocity')316 iret = NF90_PUT_ATT(ncid,vvelid, 'units', 'm/s')317 iret = NF90_PUT_ATT(ncid,massid, 'long_name', 'mass')318 iret = NF90_PUT_ATT(ncid,massid, 'units', 'kg')319 iret = NF90_PUT_ATT(ncid,thicknessid, 'long_name', 'thickness')320 iret = NF90_PUT_ATT(ncid,thicknessid, 'units', 'm')321 iret = NF90_PUT_ATT(ncid,widthid, 'long_name', 'width')322 iret = NF90_PUT_ATT(ncid,widthid, 'units', 'm')323 iret = NF90_PUT_ATT(ncid,lengthid, 'long_name', 'length')324 iret = NF90_PUT_ATT(ncid,lengthid, 'units', 'm')325 iret = NF90_PUT_ATT(ncid, numberid, 'long_name', 'iceberg number on this processor')326 iret = NF90_PUT_ATT(ncid, numberid, 'units', 'count')327 iret = NF90_PUT_ATT(ncid,yearid, 'long_name', 'calendar year of calving event')328 iret = NF90_PUT_ATT(ncid,yearid, 'units', 'years')329 iret = NF90_PUT_ATT(ncid,dayid, 'long_name', 'year day of calving event')330 iret = NF90_PUT_ATT(ncid,dayid, 'units', 'days')331 iret = NF90_PUT_ATT(ncid,scaling_id, 'long_name', 'scaling factor for mass of calving berg')332 iret = NF90_PUT_ATT(ncid,scaling_id, 'units', 'none')333 iret = NF90_PUT_ATT(ncid,mass_of_bits_id, 'long_name', 'mass of bergy bits')334 iret = NF90_PUT_ATT(ncid,mass_of_bits_id, 'units', 'kg')335 iret = NF90_PUT_ATT(ncid,heat_density_id, 'long_name', 'heat density')336 iret = NF90_PUT_ATT(ncid,heat_density_id, 'units', 'J/kg')307 nret = NF90_PUT_ATT(ncid, nlonid, 'long_name', 'longitude') 308 nret = NF90_PUT_ATT(ncid, nlonid, 'units', 'degrees_E') 309 nret = NF90_PUT_ATT(ncid, nlatid, 'long_name', 'latitude') 310 nret = NF90_PUT_ATT(ncid, nlatid, 'units', 'degrees_N') 311 nret = NF90_PUT_ATT(ncid, nxid, 'long_name', 'x grid box position') 312 nret = NF90_PUT_ATT(ncid, nxid, 'units', 'fractional') 313 nret = NF90_PUT_ATT(ncid, nyid, 'long_name', 'y grid box position') 314 nret = NF90_PUT_ATT(ncid, nyid, 'units', 'fractional') 315 nret = NF90_PUT_ATT(ncid, nuvelid, 'long_name', 'zonal velocity') 316 nret = NF90_PUT_ATT(ncid, nuvelid, 'units', 'm/s') 317 nret = NF90_PUT_ATT(ncid, nvvelid, 'long_name', 'meridional velocity') 318 nret = NF90_PUT_ATT(ncid, nvvelid, 'units', 'm/s') 319 nret = NF90_PUT_ATT(ncid, nmassid, 'long_name', 'mass') 320 nret = NF90_PUT_ATT(ncid, nmassid, 'units', 'kg') 321 nret = NF90_PUT_ATT(ncid, nthicknessid, 'long_name', 'thickness') 322 nret = NF90_PUT_ATT(ncid, nthicknessid, 'units', 'm') 323 nret = NF90_PUT_ATT(ncid, nwidthid, 'long_name', 'width') 324 nret = NF90_PUT_ATT(ncid, nwidthid, 'units', 'm') 325 nret = NF90_PUT_ATT(ncid, nlengthid, 'long_name', 'length') 326 nret = NF90_PUT_ATT(ncid, nlengthid, 'units', 'm') 327 nret = NF90_PUT_ATT(ncid, numberid, 'long_name', 'iceberg number on this processor') 328 nret = NF90_PUT_ATT(ncid, numberid, 'units', 'count') 329 nret = NF90_PUT_ATT(ncid, nyearid, 'long_name', 'calendar year of calving event') 330 nret = NF90_PUT_ATT(ncid, nyearid, 'units', 'years') 331 nret = NF90_PUT_ATT(ncid, ndayid, 'long_name', 'year day of calving event') 332 nret = NF90_PUT_ATT(ncid, ndayid, 'units', 'days') 333 nret = NF90_PUT_ATT(ncid, nscaling_id, 'long_name', 'scaling factor for mass of calving berg') 334 nret = NF90_PUT_ATT(ncid, nscaling_id, 'units', 'none') 335 nret = NF90_PUT_ATT(ncid, nmass_of_bits_id, 'long_name', 'mass of bergy bits') 336 nret = NF90_PUT_ATT(ncid, nmass_of_bits_id, 'units', 'kg') 337 nret = NF90_PUT_ATT(ncid, nheat_density_id, 'long_name', 'heat density') 338 nret = NF90_PUT_ATT(ncid, nheat_density_id, 'units', 'J/kg') 337 339 338 340 ENDIF ! associated(first_berg) 339 341 340 342 ! End define mode 341 iret = NF90_ENDDEF(ncid)343 nret = NF90_ENDDEF(ncid) 342 344 343 345 ! -------------------------------- 344 346 ! now write some data 345 347 346 strt3(1) = 1347 strt3(2) = 1348 lngth3(1) = jpi349 lngth3(2) = jpj350 lngth3(3) = 1351 352 DO k=1,nclasses353 griddata(:,:,1) = berg_grid%stored_ice(:,:, k)354 strt3(3) = k355 iret = NF90_PUT_VAR( ncid, siceid, griddata, strt3,lngth3 )356 IF ( iret .ne. NF90_NOERR) THEN357 IF( lwp ) WRITE(numout,*) TRIM(NF90_STRERROR( iret ))348 nstrt3(1) = 1 349 nstrt3(2) = 1 350 nlngth3(1) = jpi 351 nlngth3(2) = jpj 352 nlngth3(3) = 1 353 354 DO jn=1,nclasses 355 griddata(:,:,1) = berg_grid%stored_ice(:,:,jn) 356 nstrt3(3) = jn 357 nret = NF90_PUT_VAR( ncid, nsiceid, griddata, nstrt3, nlngth3 ) 358 IF (nret .ne. NF90_NOERR) THEN 359 IF( lwp ) WRITE(numout,*) TRIM(NF90_STRERROR( nret )) 358 360 CALL ctl_stop('icebergs, write_restart: nf_put_var stored_ice failed') 359 361 ENDIF 360 362 ENDDO 361 IF( lwp ) WRITE(numout,*) 'file: ',TRIM( filename),' var: stored_ice written'362 363 iret = NF90_PUT_VAR( ncid, kountid, kount_bergs(:) )364 IF ( iret .ne. NF90_NOERR) CALL ctl_stop('icebergs, write_restart: nf_put_var kount failed')365 366 iret = NF90_PUT_VAR( ncid,sheatid, berg_grid%stored_heat(:,:) )367 IF ( iret .ne. NF90_NOERR) CALL ctl_stop('icebergs, write_restart: nf_put_var stored_heat failed')368 IF( lwp ) WRITE(numout,*) 'file: ',TRIM( filename),' var: stored_heat written'369 370 iret = NF90_PUT_VAR( ncid, calvid , p_calving(:,:) )371 IF ( iret .ne. NF90_NOERR) CALL ctl_stop('icebergs, write_restart: nf_put_var calving failed')372 iret = NF90_PUT_VAR( ncid, calvhid, p_calving_hflx(:,:) )373 IF ( iret .ne. NF90_NOERR) CALL ctl_stop('icebergs, write_restart: nf_put_var calving_hflx failed')374 IF( lwp ) WRITE(numout,*) 'file: ',TRIM( filename),' var: calving written'363 IF( lwp ) WRITE(numout,*) 'file: ',TRIM(cl_filename),' var: stored_ice written' 364 365 nret = NF90_PUT_VAR( ncid, nkountid, num_bergs(:) ) 366 IF (nret .ne. NF90_NOERR) CALL ctl_stop('icebergs, write_restart: nf_put_var kount failed') 367 368 nret = NF90_PUT_VAR( ncid, nsheatid, berg_grid%stored_heat(:,:) ) 369 IF (nret .ne. NF90_NOERR) CALL ctl_stop('icebergs, write_restart: nf_put_var stored_heat failed') 370 IF( lwp ) WRITE(numout,*) 'file: ',TRIM(cl_filename),' var: stored_heat written' 371 372 nret = NF90_PUT_VAR( ncid, ncalvid , src_calving(:,:) ) 373 IF (nret .ne. NF90_NOERR) CALL ctl_stop('icebergs, write_restart: nf_put_var calving failed') 374 nret = NF90_PUT_VAR( ncid, ncalvhid, src_calving_hflx(:,:) ) 375 IF (nret .ne. NF90_NOERR) CALL ctl_stop('icebergs, write_restart: nf_put_var calving_hflx failed') 376 IF( lwp ) WRITE(numout,*) 'file: ',TRIM(cl_filename),' var: calving written' 375 377 376 378 IF ( ASSOCIATED(first_berg) ) THEN 377 379 378 380 ! Write variables 379 ! sga -just write out the current point of the trajectory381 ! just write out the current point of the trajectory 380 382 381 383 this => first_berg 382 k= 0384 jn = 0 383 385 DO WHILE (ASSOCIATED(this)) 384 386 pt => this%current_point 385 k=k+1386 387 iret = NF90_PUT_VAR(ncid, numberid, this%number, (/1,k/), (/nkounts,1/) )388 iret = NF90_PUT_VAR(ncid, scaling_id, this%mass_scaling, (/ k/) )389 390 iret = NF90_PUT_VAR(ncid, lonid, pt%lon, (/ k/) )391 iret = NF90_PUT_VAR(ncid, latid, pt%lat, (/ k/) )392 iret = NF90_PUT_VAR(ncid, xid, pt%xi, (/ k/) )393 iret = NF90_PUT_VAR(ncid, yid, pt%yj, (/ k/) )394 iret = NF90_PUT_VAR(ncid, uvelid, pt%uvel, (/ k/) )395 iret = NF90_PUT_VAR(ncid, vvelid, pt%vvel, (/ k/) )396 iret = NF90_PUT_VAR(ncid, massid, pt%mass, (/ k/) )397 iret = NF90_PUT_VAR(ncid, thicknessid, pt%thickness, (/ k/) )398 iret = NF90_PUT_VAR(ncid, widthid, pt%width, (/ k/) )399 iret = NF90_PUT_VAR(ncid, lengthid, pt%length, (/ k/) )400 iret = NF90_PUT_VAR(ncid, yearid, pt%year, (/ k/) )401 iret = NF90_PUT_VAR(ncid, dayid, pt%day, (/ k/) )402 iret = NF90_PUT_VAR(ncid, mass_of_bits_id, pt%mass_of_bits, (/ k/) )403 iret = NF90_PUT_VAR(ncid, heat_density_id, pt%heat_density, (/ k/) )387 jn=jn+1 388 389 nret = NF90_PUT_VAR(ncid, numberid, this%number, (/1,jn/), (/nkounts,1/) ) 390 nret = NF90_PUT_VAR(ncid, nscaling_id, this%mass_scaling, (/ jn /) ) 391 392 nret = NF90_PUT_VAR(ncid, nlonid, pt%lon, (/ jn /) ) 393 nret = NF90_PUT_VAR(ncid, nlatid, pt%lat, (/ jn /) ) 394 nret = NF90_PUT_VAR(ncid, nxid, pt%xi, (/ jn /) ) 395 nret = NF90_PUT_VAR(ncid, nyid, pt%yj, (/ jn /) ) 396 nret = NF90_PUT_VAR(ncid, nuvelid, pt%uvel, (/ jn /) ) 397 nret = NF90_PUT_VAR(ncid, nvvelid, pt%vvel, (/ jn /) ) 398 nret = NF90_PUT_VAR(ncid, nmassid, pt%mass, (/ jn /) ) 399 nret = NF90_PUT_VAR(ncid, nthicknessid, pt%thickness, (/ jn /) ) 400 nret = NF90_PUT_VAR(ncid, nwidthid, pt%width, (/ jn /) ) 401 nret = NF90_PUT_VAR(ncid, nlengthid, pt%length, (/ jn /) ) 402 nret = NF90_PUT_VAR(ncid, nyearid, pt%year, (/ jn /) ) 403 nret = NF90_PUT_VAR(ncid, ndayid, pt%day, (/ jn /) ) 404 nret = NF90_PUT_VAR(ncid, nmass_of_bits_id, pt%mass_of_bits, (/ jn /) ) 405 nret = NF90_PUT_VAR(ncid, nheat_density_id, pt%heat_density, (/ jn /) ) 404 406 405 407 this=>this%next … … 409 411 410 412 ! Finish up 411 iret = NF90_CLOSE(ncid)412 IF ( iret .ne. NF90_NOERR) CALL ctl_stop('icebergs, write_restart: nf_close failed')413 nret = NF90_CLOSE(ncid) 414 IF (nret .ne. NF90_NOERR) CALL ctl_stop('icebergs, write_restart: nf_close failed') 413 415 414 416 END SUBROUTINE icebergs_write_restart -
branches/2012/dev_r3337_NOCS10_ICB/NEMOGCM/NEMO/OPA_SRC/ICB/icbrun.F90
r3339 r3359 39 39 PRIVATE 40 40 41 PUBLIC icb_stp ! routine called in xxx.F90 module42 PUBLIC icb_end ! routine called in xxx.F90 module41 PUBLIC icb_stp ! routine called in sbcmod.F90 module 42 PUBLIC icb_end ! routine called in nemogcm.F90 module 43 43 44 44 CONTAINS … … 52 52 !! ** Method : - blah blah 53 53 !!---------------------------------------------------------------------- 54 INTEGER, INTENT( in ) :: kt54 INTEGER, INTENT( in ) :: kt 55 55 ! 56 56 INTEGER :: iyr, imon, iday, ihr, imin, isec 57 LOGICAL :: l err, sample_traj, l_budget,l_verbose57 LOGICAL :: ll_sample_traj, ll_budget, ll_verbose 58 58 !!---------------------------------------------------------------------- 59 59 !! start of timestep housekeeping 60 60 61 ktberg = kt61 nktberg = kt 62 62 63 63 ! read calving data … … 65 65 ! 66 66 CALL fld_read ( kt, 1, sf_icb ) 67 p_calving(:,:) = sf_icb(1)%fnow(:,:,1) ! calving in km^3/year (water equivalent)68 p_calving_hflx(:,:) = 0._wp ! NO heat flux for now67 src_calving(:,:) = sf_icb(1)%fnow(:,:,1) ! calving in km^3/year (water equivalent) 68 src_calving_hflx(:,:) = 0._wp ! NO heat flux for now 69 69 ! 70 70 ENDIF … … 87 87 current_yearday = yearday(imon, iday, ihr, imin, isec) 88 88 89 l _verbose = .FALSE.89 ll_verbose = .FALSE. 90 90 IF( nn_verbose_write .GT. 0 .AND. & 91 MOD(kt-1,nn_verbose_write ) == 0 ) l _verbose = nn_verbose_level >= 092 IF( l _verbose ) WRITE(numicb,9100) iyr, imon, iday, ihr, imin, isec, &91 MOD(kt-1,nn_verbose_write ) == 0 ) ll_verbose = nn_verbose_level >= 0 92 IF( ll_verbose ) WRITE(numicb,9100) iyr, imon, iday, ihr, imin, isec, & 93 93 current_year, current_yearday 94 94 9100 FORMAT('y,m,d=',3i5,' h,m,s=',3i5,' yr,yrdy=',i5,f8.3) … … 112 112 113 113 IF( lk_mpp ) THEN 114 114 CALL mpp_send_bergs () ! Send bergs to other PEs 115 115 ELSE 116 116 CALL lbc_send_bergs() ! Deal with any cyclic boundaries in non-mpp case 117 117 ENDIF 118 118 … … 122 122 !! end of timestep housekeeping 123 123 124 sample_traj = .FALSE.125 IF( nn_sample_rate .GT. 0 .AND. MOD(kt-1,nn_sample_rate) == 0 ) sample_traj = .TRUE.126 IF( sample_traj .AND. &124 ll_sample_traj = .FALSE. 125 IF( nn_sample_rate .GT. 0 .AND. MOD(kt-1,nn_sample_rate) == 0 ) ll_sample_traj = .TRUE. 126 IF( ll_sample_traj .AND. & 127 127 ASSOCIATED(first_berg) ) CALL traj_write ( kt ) ! For each berg, record trajectory 128 128 … … 142 142 143 143 ! Diagnose budgets 144 l _budget = .FALSE.145 IF( nn_verbose_write .GT. 0 .AND. MOD(kt-1,nn_verbose_write) == 0 ) l _budget = ln_bergdia146 CALL icb_budget( l _budget )144 ll_budget = .FALSE. 145 IF( nn_verbose_write .GT. 0 .AND. MOD(kt-1,nn_verbose_write) == 0 ) ll_budget = ln_bergdia 146 CALL icb_budget( ll_budget ) 147 147 148 148 IF( MOD(kt,nn_stock) == 0 ) THEN … … 175 175 DEALLOCATE( ssh_e ) 176 176 177 DEALLOCATE( icbfldpts ) 177 DEALLOCATE( nicbfldpts ) 178 DEALLOCATE( nicbflddest ) 179 DEALLOCATE( nicbfldproc ) 178 180 179 181 IF( lk_mpp ) CALL dealloc_buffers() … … 206 208 DEALLOCATE(berg_grid) 207 209 208 DEALLOCATE( initial_width)209 DEALLOCATE( initial_length)210 DEALLOCATE(first_width) 211 DEALLOCATE(first_length) 210 212 211 213 IF (lwp) WRITE(numout,'(a,i6)') 'icebergs: icb_end complete',narea -
branches/2012/dev_r3337_NOCS10_ICB/NEMOGCM/NEMO/OPA_SRC/ICB/icbthm.F90
r3339 r3359 29 29 PRIVATE 30 30 31 PUBLIC thermodynamics ! routine called in xxx.F90 module31 PUBLIC thermodynamics ! routine called in icbrun.F90 module 32 32 33 33 CONTAINS … … 43 43 INTEGER :: kt ! timestep number, just passed to print_berg 44 44 ! 45 REAL(wp) :: M, T, W, L, SST, Vol, Ln, Wn, Tn, nVol, IC,Dn46 REAL(wp) :: Mv, Me, Mb, melt, dvo, dva, dM, Ss, dMe, dMb,dMv47 REAL(wp) :: Mnew, Mnew1, Mnew2,heat48 REAL(wp) :: Mbits, nMbits, dMbitsE, dMbitsM, Lbits, Abits,Mbb49 REAL(wp) :: xi,yj, zff, z1_rday, z1_e1e2, zdt, z1_dt, z1_dt_e1e245 REAL(wp) :: zM, zT, zW, zL, zSST, zVol, zLn, zWn, zTn, znVol, zIC, zDn 46 REAL(wp) :: zMv, zMe, zMb, zmelt, zdvo, zdva, zdM, zSs, zdMe, zdMb, zdMv 47 REAL(wp) :: zMnew, zMnew1, zMnew2, zheat 48 REAL(wp) :: zMbits, znMbits, zdMbitsE, zdMbitsM, zLbits, zAbits, zMbb 49 REAL(wp) :: zxi, zyj, zff, z1_rday, z1_e1e2, zdt, z1_dt, z1_dt_e1e2 50 50 INTEGER :: ii, ij 51 51 TYPE(iceberg) , POINTER :: this, next … … 67 67 ! 68 68 pt => this%current_point 69 knberg = this%number(1)69 nknberg = this%number(1) 70 70 CALL interp_flds( pt%xi, pt%e1, pt%uo, pt%ui, pt%ua, pt%ssh_x, & 71 & pt%yj, pt%e2, pt%vo, pt%vi, pt%va, pt%ssh_y, pt%sst, pt%cn, pt%hi, zff ) 72 ! 73 SST = pt%sst 74 IC = MIN( 1._wp, pt%cn + rn_sicn_shift ) ! Shift sea-ice concentration !!gm ??? 75 M = pt%mass 76 T = pt%thickness ! total thickness 77 ! D = (rn_rho_bergs/rho_seawater)*T ! draught (keel depth) 78 ! F = T - D ! freeboard 79 W = pt%width 80 L = pt%length 81 xi = pt%xi ! position in (i,j) referential 82 yj = pt%yj 83 ii = INT( xi + 0.5 ) - nimpp + 1 ! t-cell of the berg 84 ij = INT( yj + 0.5 ) - njmpp + 1 85 Vol = T * W * L 71 & pt%yj, pt%e2, pt%vo, pt%vi, pt%va, pt%ssh_y, & 72 & pt%sst, pt%cn, pt%hi, zff ) 73 ! 74 zSST = pt%sst 75 zIC = MIN( 1._wp, pt%cn + rn_sicn_shift ) ! Shift sea-ice concentration !!gm ??? 76 zM = pt%mass 77 zT = pt%thickness ! total thickness 78 ! D = (rn_rho_bergs/pp_rho_seawater)*zT ! draught (keel depth) 79 ! F = zT - D ! freeboard 80 zW = pt%width 81 zL = pt%length 82 zxi = pt%xi ! position in (i,j) referential 83 zyj = pt%yj 84 ii = INT( zxi + 0.5 ) - nimpp + 1 ! t-cell of the berg 85 ij = INT( zyj + 0.5 ) - njmpp + 1 86 zVol = zT * zW * zL 86 87 zdt = berg_dt ; z1_dt = 1._wp / zdt 87 88 88 89 ! Environment 89 dvo = SQRT( (pt%uvel-pt%uo)**2 + (pt%vvel-pt%vo)**2 )90 dva = SQRT( (pt%ua -pt%uo)**2 + (pt%va -pt%vo)**2 )91 Ss = 1.5 * SQRT( dva ) + 0.1 *dva ! Sea state (eqn M.A9)90 zdvo = SQRT( (pt%uvel-pt%uo)**2 + (pt%vvel-pt%vo)**2 ) 91 zdva = SQRT( (pt%ua -pt%uo)**2 + (pt%va -pt%vo)**2 ) 92 zSs = 1.5 * SQRT( zdva ) + 0.1 * zdva ! Sea state (eqn M.A9) 92 93 93 94 ! Melt rates in m/s (i.e. division by rday) 94 Mv = MAX( 7.62e-3*SST+1.29e-3*(SST**2) , 0._wp ) * z1_rday ! Buoyant convection at sides (eqn M.A10)95 Mb = MAX( 0.58*(dvo**0.8)*(SST+4.0)/(L**0.2) , 0._wp ) * z1_rday ! Basal turbulent melting (eqn M.A7 )96 Me = MAX( 1./12.*(SST+2.)*Ss*(1+cos(rpi*(IC**3))) , 0._wp ) * z1_rday ! Wave erosion (eqn M.A8 )95 zMv = MAX( 7.62e-3*zSST+1.29e-3*(zSST**2) , 0._wp ) * z1_rday ! Buoyant convection at sides (eqn M.A10) 96 zMb = MAX( 0.58*(zdvo**0.8)*(zSST+4.0)/(zL**0.2) , 0._wp ) * z1_rday ! Basal turbulent melting (eqn M.A7 ) 97 zMe = MAX( 1./12.*(zSST+2.)*zSs*(1+cos(rpi*(zIC**3))) , 0._wp ) * z1_rday ! Wave erosion (eqn M.A8 ) 97 98 98 99 IF( ln_operator_splitting ) THEN ! Operator split update of volume/mass 99 Tn = MAX( T -Mb*zdt , 0._wp ) ! new total thickness (m)100 nVol = Tn * W *L ! new volume (m^3)101 Mnew1 = (nVol/Vol) *M ! new mass (kg)102 dMb = M -Mnew1 ! mass lost to basal melting (>0) (kg)103 ! 104 Ln = MAX( L -Mv*zdt , 0._wp ) ! new length (m)105 Wn = MAX( W -Mv*zdt , 0._wp ) ! new width (m)106 nVol = Tn * Wn *Ln ! new volume (m^3)107 Mnew2 = (nVol/Vol) *M ! new mass (kg)108 dMv = Mnew1 -Mnew2 ! mass lost to buoyant convection (>0) (kg)109 ! 110 Ln = MAX( Ln -Me*zdt , 0._wp ) ! new length (m)111 Wn = MAX( Wn -Me*zdt , 0._wp ) ! new width (m)112 nVol = Tn * Wn *Ln ! new volume (m^3)113 Mnew = ( nVol / Vol ) *M ! new mass (kg)114 dMe = Mnew2 -Mnew ! mass lost to erosion (>0) (kg)115 dM = M -Mnew ! mass lost to all erosion and melting (>0) (kg)100 zTn = MAX( zT - zMb*zdt , 0._wp ) ! new total thickness (m) 101 znVol = zTn * zW * zL ! new volume (m^3) 102 zMnew1 = (znVol/zVol) * zM ! new mass (kg) 103 zdMb = zM - zMnew1 ! mass lost to basal melting (>0) (kg) 104 ! 105 zLn = MAX( zL - zMv*zdt , 0._wp ) ! new length (m) 106 zWn = MAX( zW - zMv*zdt , 0._wp ) ! new width (m) 107 znVol = zTn * zWn * zLn ! new volume (m^3) 108 zMnew2 = (znVol/zVol) * zM ! new mass (kg) 109 zdMv = zMnew1 - zMnew2 ! mass lost to buoyant convection (>0) (kg) 110 ! 111 zLn = MAX( zLn - zMe*zdt , 0._wp ) ! new length (m) 112 zWn = MAX( zWn - zMe*zdt , 0._wp ) ! new width (m) 113 znVol = zTn * zWn * zLn ! new volume (m^3) 114 zMnew = ( znVol / zVol ) * zM ! new mass (kg) 115 zdMe = zMnew2 - zMnew ! mass lost to erosion (>0) (kg) 116 zdM = zM - zMnew ! mass lost to all erosion and melting (>0) (kg) 116 117 ! 117 118 ELSE ! Update dimensions of berg 118 Ln = MAX( L -(Mv+Me)*zdt ,0._wp ) ! (m)119 Wn = MAX( W -(Mv+Me)*zdt ,0._wp ) ! (m)120 Tn = MAX( T -Mb *zdt ,0._wp ) ! (m)119 zLn = MAX( zL -(zMv+zMe)*zdt ,0._wp ) ! (m) 120 zWn = MAX( zW -(zMv+zMe)*zdt ,0._wp ) ! (m) 121 zTn = MAX( zT - zMb *zdt ,0._wp ) ! (m) 121 122 ! Update volume and mass of berg 122 nVol = Tn*Wn*Ln ! (m^3)123 Mnew = (nVol/Vol)*M ! (kg)124 dM = M -Mnew ! (kg)125 dMb = (M/Vol) * (W* L ) *Mb*zdt ! approx. mass loss to basal melting (kg)126 dMe = (M/Vol) * (T*(W+L)) *Me*zdt ! approx. mass lost to erosion (kg)127 dMv = (M/Vol) * (T*(W+L)) *Mv*zdt ! approx. mass loss to buoyant convection (kg)123 znVol = zTn*zWn*zLn ! (m^3) 124 zMnew = (znVol/zVol)*zM ! (kg) 125 zdM = zM - zMnew ! (kg) 126 zdMb = (zM/zVol) * (zW* zL ) *zMb*zdt ! approx. mass loss to basal melting (kg) 127 zdMe = (zM/zVol) * (zT*(zW+zL)) *zMe*zdt ! approx. mass lost to erosion (kg) 128 zdMv = (zM/zVol) * (zT*(zW+zL)) *zMv*zdt ! approx. mass loss to buoyant convection (kg) 128 129 ENDIF 129 130 130 131 IF( rn_bits_erosion_fraction > 0._wp ) THEN ! Bergy bits 131 132 ! 132 Mbits = pt%mass_of_bits ! mass of bergy bits (kg)133 dMbitsE = rn_bits_erosion_fraction *dMe ! change in mass of bits (kg)134 nMbits = Mbits +dMbitsE ! add new bergy bits to mass (kg)135 Lbits = MIN( L, W,T, 40._wp ) ! assume bergy bits are smallest dimension or 40 meters136 Abits = ( Mbits / rn_rho_bergs ) /Lbits ! Effective bottom area (assuming T=Lbits)137 Mbb = MAX( 0.58*(dvo**0.8)*(SST+2.0)/(Lbits**0.2), 0.) * z1_rday ! Basal turbulent melting (for bits)138 Mbb = rn_rho_bergs * Abits *Mbb ! in kg/s139 dMbitsM = MIN( Mbb*zdt ,nMbits ) ! bergy bits mass lost to melting (kg)140 nMbits = nMbits-dMbitsM ! remove mass lost to bergy bits melt141 IF( Mnew == 0._wp ) THEN ! if parent berg has completely melted then142 dMbitsM = dMbitsM +nMbits ! instantly melt all the bergy bits143 nMbits = 0._wp133 zMbits = pt%mass_of_bits ! mass of bergy bits (kg) 134 zdMbitsE = rn_bits_erosion_fraction * zdMe ! change in mass of bits (kg) 135 znMbits = zMbits + zdMbitsE ! add new bergy bits to mass (kg) 136 zLbits = MIN( zL, zW, zT, 40._wp ) ! assume bergy bits are smallest dimension or 40 meters 137 zAbits = ( zMbits / rn_rho_bergs ) / zLbits ! Effective bottom area (assuming T=Lbits) 138 zMbb = MAX( 0.58*(zdvo**0.8)*(zSST+2.0)/(zLbits**0.2), 0.) * z1_rday ! Basal turbulent melting (for bits) 139 zMbb = rn_rho_bergs * zAbits * zMbb ! in kg/s 140 zdMbitsM = MIN( zMbb*zdt , znMbits ) ! bergy bits mass lost to melting (kg) 141 znMbits = znMbits-zdMbitsM ! remove mass lost to bergy bits melt 142 IF( zMnew == 0._wp ) THEN ! if parent berg has completely melted then 143 zdMbitsM = zdMbitsM + znMbits ! instantly melt all the bergy bits 144 znMbits = 0._wp 144 145 ENDIF 145 146 ELSE ! No bergy bits 146 Abits = 0._wp147 dMbitsE = 0._wp148 dMbitsM = 0._wp149 nMbits = pt%mass_of_bits ! retain previous value incase non-zero147 zAbits = 0._wp 148 zdMbitsE = 0._wp 149 zdMbitsM = 0._wp 150 znMbits = pt%mass_of_bits ! retain previous value incase non-zero 150 151 ENDIF 151 152 … … 154 155 z1_e1e2 = 1._wp / e1e2t(ii,ij) * this%mass_scaling 155 156 z1_dt_e1e2 = z1_dt * z1_e1e2 156 melt = ( dM - ( dMbitsE - dMbitsM ) ) * z1_dt ! kg/s 157 berg_grid%floating_melt(ii,ij) = berg_grid%floating_melt(ii,ij) + melt * z1_e1e2 ! kg/m2/s 158 heat = melt * pt%heat_density ! kg/s x J/kg = J/s 159 berg_grid%calving_hflx (ii,ij) = berg_grid%calving_hflx (ii,ij) + heat * z1_e1e2 ! W/m2 160 CALL melt_budget(ii, ij, Mnew, heat, this%mass_scaling, dM, dMbitsE, dMbitsM, dMb, dMe, dMv, z1_dt_e1e2 ) 157 zmelt = ( zdM - ( zdMbitsE - zdMbitsM ) ) * z1_dt ! kg/s 158 berg_grid%floating_melt(ii,ij) = berg_grid%floating_melt(ii,ij) + zmelt * z1_e1e2 ! kg/m2/s 159 zheat = zmelt * pt%heat_density ! kg/s x J/kg = J/s 160 berg_grid%calving_hflx (ii,ij) = berg_grid%calving_hflx (ii,ij) + zheat * z1_e1e2 ! W/m2 161 CALL melt_budget( ii, ij, zMnew, zheat, this%mass_scaling, & 162 & zdM, zdMbitsE, zdMbitsM, zdMb, zdMe, & 163 & zdMv, z1_dt_e1e2 ) 161 164 ELSE 162 165 WRITE(numout,*) 'thermodynamics: berg ',this%number(:),' appears to have grounded at ',narea,ii,ij … … 167 170 168 171 ! Rolling 169 Dn = ( rn_rho_bergs / rho_seawater ) *Tn ! draught (keel depth)170 IF( Dn > 0._wp .AND. MAX(Wn,Ln) < SQRT( 0.92*(Dn**2) + 58.32*Dn ) ) THEN171 T =Tn172 Tn =Wn173 Wn =T172 zDn = ( rn_rho_bergs / pp_rho_seawater ) * zTn ! draught (keel depth) 173 IF( zDn > 0._wp .AND. MAX(zWn,zLn) < SQRT( 0.92*(zDn**2) + 58.32*zDn ) ) THEN 174 zT = zTn 175 zTn = zWn 176 zWn = zT 174 177 endif 175 178 176 179 ! Store the new state of iceberg (with L>W) 177 pt%mass = Mnew178 pt%mass_of_bits = nMbits179 pt%thickness = Tn180 pt%width = min( Wn,Ln)181 pt%length = max( Wn,Ln)180 pt%mass = zMnew 181 pt%mass_of_bits = znMbits 182 pt%thickness = zTn 183 pt%width = min(zWn,zLn) 184 pt%length = max(zWn,zLn) 182 185 183 186 next=>this%next … … 185 188 !!gm add a test to avoid over melting ? 186 189 187 IF( Mnew <= 0._wp ) THEN ! Delete the berg if completely melted190 IF( zMnew <= 0._wp ) THEN ! Delete the berg if completely melted 188 191 CALL delete_iceberg_from_list( first_berg, this ) 189 192 ! 190 193 ELSE ! Diagnose mass distribution on grid 191 194 z1_e1e2 = 1._wp / e1e2t(ii,ij) * this%mass_scaling 192 CALL size_budget(ii, ij, Wn, Ln, Abits, this%mass_scaling, Mnew, nMbits, z1_e1e2) 195 CALL size_budget( ii, ij, zWn, zLn, zAbits, & 196 & this%mass_scaling, zMnew, znMbits, z1_e1e2) 193 197 ENDIF 194 198 ! -
branches/2012/dev_r3337_NOCS10_ICB/NEMOGCM/NEMO/OPA_SRC/ICB/icbtrj.F90
r3339 r3359 26 26 PRIVATE 27 27 28 INTEGER, PRIVATE :: kount_traj29 INTEGER, PRIVATE :: n_dim, k_dim30 INTEGER, PRIVATE :: trajid31 INTEGER, PRIVATE :: numberid, stepid,scaling_id32 INTEGER, PRIVATE :: lonid, latid, xid, yid, uvelid, vvelid,massid33 INTEGER, PRIVATE :: uoid, void, uaid, vaid, uiid,viid34 INTEGER, PRIVATE :: sshxid, sshyid, sstid, cntid,thkid35 INTEGER, PRIVATE :: thicknessid, widthid,lengthid36 INTEGER, PRIVATE :: yearid,dayid37 INTEGER, PRIVATE :: mass_of_bits_id,heat_density_id38 39 PUBLIC traj_init ! routine called in xxx.F90 module40 PUBLIC traj_write ! routine called in xxx.F90 module41 PUBLIC traj_sync ! routine called in xxx.F90 module42 PUBLIC traj_end ! routine called in xxx.F90 module28 INTEGER, PRIVATE :: num_traj 29 INTEGER, PRIVATE :: n_dim, m_dim 30 INTEGER, PRIVATE :: ntrajid 31 INTEGER, PRIVATE :: numberid, nstepid, nscaling_id 32 INTEGER, PRIVATE :: nlonid, nlatid, nxid, nyid, nuvelid, nvvelid, nmassid 33 INTEGER, PRIVATE :: nuoid, nvoid, nuaid, nvaid, nuiid, nviid 34 INTEGER, PRIVATE :: nsshxid, nsshyid, nsstid, ncntid, nthkid 35 INTEGER, PRIVATE :: nthicknessid, nwidthid, nlengthid 36 INTEGER, PRIVATE :: nyearid, ndayid 37 INTEGER, PRIVATE :: nmass_of_bits_id, nheat_density_id 38 39 PUBLIC traj_init ! routine called in icbini.F90 module 40 PUBLIC traj_write ! routine called in icbrun.F90 module 41 PUBLIC traj_sync ! routine called in icbrun.F90 module 42 PUBLIC traj_end ! routine called in icbrun.F90 module 43 43 44 44 CONTAINS … … 46 46 !!------------------------------------------------------------------------- 47 47 48 SUBROUTINE traj_init( nitend )48 SUBROUTINE traj_init( ktend ) 49 49 50 50 ! local variables 51 INTEGER, INTENT( in ) :: nitend52 ! Local variables 53 INTEGER :: iret , k54 CHARACTER(len=80) :: filename51 INTEGER, INTENT( in ) :: ktend 52 ! Local variables 53 INTEGER :: iret 54 CHARACTER(len=80) :: cl_filename 55 55 TYPE(iceberg), POINTER :: this 56 56 TYPE(point) , POINTER :: pt 57 57 58 58 IF( lk_mpp ) THEN 59 WRITE( filename,'("trajectory_icebergs_",I6.6,"_",I4.4,".nc")') nitend, narea-159 WRITE(cl_filename,'("trajectory_icebergs_",I6.6,"_",I4.4,".nc")') ktend, narea-1 60 60 ELSE 61 WRITE( filename,'("trajectory_icebergs_",I6.6,".nc")') nitend61 WRITE(cl_filename,'("trajectory_icebergs_",I6.6,".nc")') ktend 62 62 ENDIF 63 IF ( lwp .AND. nn_verbose_level >= 0) WRITE(numout,'(2a)') 'icebergs, traj_init: creating ',TRIM( filename)64 65 iret = NF90_CREATE(TRIM( filename), NF90_CLOBBER,trajid)63 IF ( lwp .AND. nn_verbose_level >= 0) WRITE(numout,'(2a)') 'icebergs, traj_init: creating ',TRIM(cl_filename) 64 65 iret = NF90_CREATE(TRIM(cl_filename), NF90_CLOBBER, ntrajid) 66 66 IF (iret .NE. NF90_NOERR) CALL ctl_stop('icebergs, traj_init: nf_create failed') 67 67 68 68 ! Dimensions 69 iret = NF90_DEF_DIM( trajid, 'n', NF90_UNLIMITED, n_dim)69 iret = NF90_DEF_DIM(ntrajid, 'n', NF90_UNLIMITED, n_dim) 70 70 IF (iret .NE. NF90_NOERR) CALL ctl_stop('icebergs, traj_init: nf_def_dim n failed') 71 iret = NF90_DEF_DIM( trajid, 'k', nkounts, k_dim)71 iret = NF90_DEF_DIM(ntrajid, 'k', nkounts, m_dim) 72 72 IF (iret .NE. NF90_NOERR) CALL ctl_stop('icebergs, traj_init: nf_def_dim k failed') 73 73 74 74 ! Variables 75 iret = NF90_DEF_VAR( trajid, 'iceberg_number', NF90_INT, (/k_dim,n_dim/), numberid)76 iret = NF90_DEF_VAR( trajid, 'timestep', NF90_INT, n_dim,stepid)77 iret = NF90_DEF_VAR( trajid, 'mass_scaling', NF90_DOUBLE, n_dim,scaling_id)78 iret = NF90_DEF_VAR( trajid, 'lon', NF90_DOUBLE, n_dim,lonid)79 iret = NF90_DEF_VAR( trajid, 'lat', NF90_DOUBLE, n_dim,latid)80 iret = NF90_DEF_VAR( trajid, 'xi', NF90_DOUBLE, n_dim,xid)81 iret = NF90_DEF_VAR( trajid, 'yj', NF90_DOUBLE, n_dim,yid)82 iret = NF90_DEF_VAR( trajid, 'uvel', NF90_DOUBLE, n_dim,uvelid)83 iret = NF90_DEF_VAR( trajid, 'vvel', NF90_DOUBLE, n_dim,vvelid)84 iret = NF90_DEF_VAR( trajid, 'uto', NF90_DOUBLE, n_dim,uoid)85 iret = NF90_DEF_VAR( trajid, 'vto', NF90_DOUBLE, n_dim,void)86 iret = NF90_DEF_VAR( trajid, 'uta', NF90_DOUBLE, n_dim,uaid)87 iret = NF90_DEF_VAR( trajid, 'vta', NF90_DOUBLE, n_dim,vaid)88 iret = NF90_DEF_VAR( trajid, 'uti', NF90_DOUBLE, n_dim,uiid)89 iret = NF90_DEF_VAR( trajid, 'vti', NF90_DOUBLE, n_dim,viid)90 iret = NF90_DEF_VAR( trajid, 'ssh_x', NF90_DOUBLE, n_dim,sshyid)91 iret = NF90_DEF_VAR( trajid, 'ssh_y', NF90_DOUBLE, n_dim,sshxid)92 iret = NF90_DEF_VAR( trajid, 'sst', NF90_DOUBLE, n_dim,sstid)93 iret = NF90_DEF_VAR( trajid, 'icnt', NF90_DOUBLE, n_dim,cntid)94 iret = NF90_DEF_VAR( trajid, 'ithk', NF90_DOUBLE, n_dim,thkid)95 iret = NF90_DEF_VAR( trajid, 'mass', NF90_DOUBLE, n_dim,massid)96 iret = NF90_DEF_VAR( trajid, 'thickness', NF90_DOUBLE, n_dim,thicknessid)97 iret = NF90_DEF_VAR( trajid, 'width', NF90_DOUBLE, n_dim,widthid)98 iret = NF90_DEF_VAR( trajid, 'length', NF90_DOUBLE, n_dim,lengthid)99 iret = NF90_DEF_VAR( trajid, 'year', NF90_INT, n_dim,yearid)100 iret = NF90_DEF_VAR( trajid, 'day', NF90_DOUBLE, n_dim,dayid)101 iret = NF90_DEF_VAR( trajid, 'mass_of_bits', NF90_DOUBLE, n_dim,mass_of_bits_id)102 iret = NF90_DEF_VAR( trajid, 'heat_density', NF90_DOUBLE, n_dim,heat_density_id)75 iret = NF90_DEF_VAR(ntrajid, 'iceberg_number', NF90_INT, (/m_dim,n_dim/), numberid) 76 iret = NF90_DEF_VAR(ntrajid, 'timestep', NF90_INT, n_dim, nstepid) 77 iret = NF90_DEF_VAR(ntrajid, 'mass_scaling', NF90_DOUBLE, n_dim, nscaling_id) 78 iret = NF90_DEF_VAR(ntrajid, 'lon', NF90_DOUBLE, n_dim, nlonid) 79 iret = NF90_DEF_VAR(ntrajid, 'lat', NF90_DOUBLE, n_dim, nlatid) 80 iret = NF90_DEF_VAR(ntrajid, 'xi', NF90_DOUBLE, n_dim, nxid) 81 iret = NF90_DEF_VAR(ntrajid, 'yj', NF90_DOUBLE, n_dim, nyid) 82 iret = NF90_DEF_VAR(ntrajid, 'uvel', NF90_DOUBLE, n_dim, nuvelid) 83 iret = NF90_DEF_VAR(ntrajid, 'vvel', NF90_DOUBLE, n_dim, nvvelid) 84 iret = NF90_DEF_VAR(ntrajid, 'uto', NF90_DOUBLE, n_dim, nuoid) 85 iret = NF90_DEF_VAR(ntrajid, 'vto', NF90_DOUBLE, n_dim, nvoid) 86 iret = NF90_DEF_VAR(ntrajid, 'uta', NF90_DOUBLE, n_dim, nuaid) 87 iret = NF90_DEF_VAR(ntrajid, 'vta', NF90_DOUBLE, n_dim, nvaid) 88 iret = NF90_DEF_VAR(ntrajid, 'uti', NF90_DOUBLE, n_dim, nuiid) 89 iret = NF90_DEF_VAR(ntrajid, 'vti', NF90_DOUBLE, n_dim, nviid) 90 iret = NF90_DEF_VAR(ntrajid, 'ssh_x', NF90_DOUBLE, n_dim, nsshyid) 91 iret = NF90_DEF_VAR(ntrajid, 'ssh_y', NF90_DOUBLE, n_dim, nsshxid) 92 iret = NF90_DEF_VAR(ntrajid, 'sst', NF90_DOUBLE, n_dim, nsstid) 93 iret = NF90_DEF_VAR(ntrajid, 'icnt', NF90_DOUBLE, n_dim, ncntid) 94 iret = NF90_DEF_VAR(ntrajid, 'ithk', NF90_DOUBLE, n_dim, nthkid) 95 iret = NF90_DEF_VAR(ntrajid, 'mass', NF90_DOUBLE, n_dim, nmassid) 96 iret = NF90_DEF_VAR(ntrajid, 'thickness', NF90_DOUBLE, n_dim, nthicknessid) 97 iret = NF90_DEF_VAR(ntrajid, 'width', NF90_DOUBLE, n_dim, nwidthid) 98 iret = NF90_DEF_VAR(ntrajid, 'length', NF90_DOUBLE, n_dim, nlengthid) 99 iret = NF90_DEF_VAR(ntrajid, 'year', NF90_INT, n_dim, nyearid) 100 iret = NF90_DEF_VAR(ntrajid, 'day', NF90_DOUBLE, n_dim, ndayid) 101 iret = NF90_DEF_VAR(ntrajid, 'mass_of_bits', NF90_DOUBLE, n_dim, nmass_of_bits_id) 102 iret = NF90_DEF_VAR(ntrajid, 'heat_density', NF90_DOUBLE, n_dim, nheat_density_id) 103 103 104 104 ! Attributes 105 iret = NF90_PUT_ATT( trajid, numberid, 'long_name', 'iceberg number on this processor')106 iret = NF90_PUT_ATT( trajid, numberid, 'units', 'count')107 iret = NF90_PUT_ATT( trajid,stepid, 'long_name', 'timestep number kt')108 iret = NF90_PUT_ATT( trajid,stepid, 'units', 'count')109 iret = NF90_PUT_ATT( trajid,lonid, 'long_name', 'longitude')110 iret = NF90_PUT_ATT( trajid,lonid, 'units', 'degrees_E')111 iret = NF90_PUT_ATT( trajid,latid, 'long_name', 'latitude')112 iret = NF90_PUT_ATT( trajid,latid, 'units', 'degrees_N')113 iret = NF90_PUT_ATT( trajid,xid, 'long_name', 'x grid box position')114 iret = NF90_PUT_ATT( trajid,xid, 'units', 'fractional')115 iret = NF90_PUT_ATT( trajid,yid, 'long_name', 'y grid box position')116 iret = NF90_PUT_ATT( trajid,yid, 'units', 'fractional')117 iret = NF90_PUT_ATT( trajid,uvelid, 'long_name', 'zonal velocity')118 iret = NF90_PUT_ATT( trajid,uvelid, 'units', 'm/s')119 iret = NF90_PUT_ATT( trajid,vvelid, 'long_name', 'meridional velocity')120 iret = NF90_PUT_ATT( trajid,vvelid, 'units', 'm/s')121 iret = NF90_PUT_ATT( trajid,uoid, 'long_name', 'ocean u component')122 iret = NF90_PUT_ATT( trajid,uoid, 'units', 'm/s')123 iret = NF90_PUT_ATT( trajid,void, 'long_name', 'ocean v component')124 iret = NF90_PUT_ATT( trajid,void, 'units', 'm/s')125 iret = NF90_PUT_ATT( trajid,uaid, 'long_name', 'atmosphere u component')126 iret = NF90_PUT_ATT( trajid,uaid, 'units', 'm/s')127 iret = NF90_PUT_ATT( trajid,vaid, 'long_name', 'atmosphere v component')128 iret = NF90_PUT_ATT( trajid,vaid, 'units', 'm/s')129 iret = NF90_PUT_ATT( trajid,uiid, 'long_name', 'sea ice u component')130 iret = NF90_PUT_ATT( trajid,uiid, 'units', 'm/s')131 iret = NF90_PUT_ATT( trajid,viid, 'long_name', 'sea ice v component')132 iret = NF90_PUT_ATT( trajid,viid, 'units', 'm/s')133 iret = NF90_PUT_ATT( trajid,sshxid, 'long_name', 'sea surface height gradient from x points')134 iret = NF90_PUT_ATT( trajid,sshxid, 'units', 'm/m')135 iret = NF90_PUT_ATT( trajid,sshyid, 'long_name', 'sea surface height gradient from y points')136 iret = NF90_PUT_ATT( trajid,sshyid, 'units', 'm/m')137 iret = NF90_PUT_ATT( trajid,sstid, 'long_name', 'sea surface temperature')138 iret = NF90_PUT_ATT( trajid,sstid, 'units', 'degC')139 iret = NF90_PUT_ATT( trajid,cntid, 'long_name', 'sea ice concentration')140 iret = NF90_PUT_ATT( trajid,cntid, 'units', 'degC')141 iret = NF90_PUT_ATT( trajid,thkid, 'long_name', 'sea ice thickness')142 iret = NF90_PUT_ATT( trajid,thkid, 'units', 'm')143 iret = NF90_PUT_ATT( trajid,massid, 'long_name', 'mass')144 iret = NF90_PUT_ATT( trajid,massid, 'units', 'kg')145 iret = NF90_PUT_ATT( trajid,thicknessid, 'long_name', 'thickness')146 iret = NF90_PUT_ATT( trajid,thicknessid, 'units', 'm')147 iret = NF90_PUT_ATT( trajid,widthid, 'long_name', 'width')148 iret = NF90_PUT_ATT( trajid,widthid, 'units', 'm')149 iret = NF90_PUT_ATT( trajid,lengthid, 'long_name', 'length')150 iret = NF90_PUT_ATT( trajid,lengthid, 'units', 'm')151 iret = NF90_PUT_ATT( trajid,yearid, 'long_name', 'calendar year')152 iret = NF90_PUT_ATT( trajid,yearid, 'units', 'years')153 iret = NF90_PUT_ATT( trajid,dayid, 'long_name', 'day of year')154 iret = NF90_PUT_ATT( trajid,dayid, 'units', 'days')155 iret = NF90_PUT_ATT( trajid,scaling_id, 'long_name', 'scaling factor for mass of berg')156 iret = NF90_PUT_ATT( trajid,scaling_id, 'units', 'none')157 iret = NF90_PUT_ATT( trajid,mass_of_bits_id, 'long_name', 'mass of bergy bits')158 iret = NF90_PUT_ATT( trajid,mass_of_bits_id, 'units', 'kg')159 iret = NF90_PUT_ATT( trajid,heat_density_id, 'long_name', 'heat density')160 iret = NF90_PUT_ATT( trajid,heat_density_id, 'units', 'J/kg')105 iret = NF90_PUT_ATT(ntrajid, numberid, 'long_name', 'iceberg number on this processor') 106 iret = NF90_PUT_ATT(ntrajid, numberid, 'units', 'count') 107 iret = NF90_PUT_ATT(ntrajid, nstepid, 'long_name', 'timestep number kt') 108 iret = NF90_PUT_ATT(ntrajid, nstepid, 'units', 'count') 109 iret = NF90_PUT_ATT(ntrajid, nlonid, 'long_name', 'longitude') 110 iret = NF90_PUT_ATT(ntrajid, nlonid, 'units', 'degrees_E') 111 iret = NF90_PUT_ATT(ntrajid, nlatid, 'long_name', 'latitude') 112 iret = NF90_PUT_ATT(ntrajid, nlatid, 'units', 'degrees_N') 113 iret = NF90_PUT_ATT(ntrajid, nxid, 'long_name', 'x grid box position') 114 iret = NF90_PUT_ATT(ntrajid, nxid, 'units', 'fractional') 115 iret = NF90_PUT_ATT(ntrajid, nyid, 'long_name', 'y grid box position') 116 iret = NF90_PUT_ATT(ntrajid, nyid, 'units', 'fractional') 117 iret = NF90_PUT_ATT(ntrajid, nuvelid, 'long_name', 'zonal velocity') 118 iret = NF90_PUT_ATT(ntrajid, nuvelid, 'units', 'm/s') 119 iret = NF90_PUT_ATT(ntrajid, nvvelid, 'long_name', 'meridional velocity') 120 iret = NF90_PUT_ATT(ntrajid, nvvelid, 'units', 'm/s') 121 iret = NF90_PUT_ATT(ntrajid, nuoid, 'long_name', 'ocean u component') 122 iret = NF90_PUT_ATT(ntrajid, nuoid, 'units', 'm/s') 123 iret = NF90_PUT_ATT(ntrajid, nvoid, 'long_name', 'ocean v component') 124 iret = NF90_PUT_ATT(ntrajid, nvoid, 'units', 'm/s') 125 iret = NF90_PUT_ATT(ntrajid, nuaid, 'long_name', 'atmosphere u component') 126 iret = NF90_PUT_ATT(ntrajid, nuaid, 'units', 'm/s') 127 iret = NF90_PUT_ATT(ntrajid, nvaid, 'long_name', 'atmosphere v component') 128 iret = NF90_PUT_ATT(ntrajid, nvaid, 'units', 'm/s') 129 iret = NF90_PUT_ATT(ntrajid, nuiid, 'long_name', 'sea ice u component') 130 iret = NF90_PUT_ATT(ntrajid, nuiid, 'units', 'm/s') 131 iret = NF90_PUT_ATT(ntrajid, nviid, 'long_name', 'sea ice v component') 132 iret = NF90_PUT_ATT(ntrajid, nviid, 'units', 'm/s') 133 iret = NF90_PUT_ATT(ntrajid, nsshxid, 'long_name', 'sea surface height gradient from x points') 134 iret = NF90_PUT_ATT(ntrajid, nsshxid, 'units', 'm/m') 135 iret = NF90_PUT_ATT(ntrajid, nsshyid, 'long_name', 'sea surface height gradient from y points') 136 iret = NF90_PUT_ATT(ntrajid, nsshyid, 'units', 'm/m') 137 iret = NF90_PUT_ATT(ntrajid, nsstid, 'long_name', 'sea surface temperature') 138 iret = NF90_PUT_ATT(ntrajid, nsstid, 'units', 'degC') 139 iret = NF90_PUT_ATT(ntrajid, ncntid, 'long_name', 'sea ice concentration') 140 iret = NF90_PUT_ATT(ntrajid, ncntid, 'units', 'degC') 141 iret = NF90_PUT_ATT(ntrajid, nthkid, 'long_name', 'sea ice thickness') 142 iret = NF90_PUT_ATT(ntrajid, nthkid, 'units', 'm') 143 iret = NF90_PUT_ATT(ntrajid, nmassid, 'long_name', 'mass') 144 iret = NF90_PUT_ATT(ntrajid, nmassid, 'units', 'kg') 145 iret = NF90_PUT_ATT(ntrajid, nthicknessid, 'long_name', 'thickness') 146 iret = NF90_PUT_ATT(ntrajid, nthicknessid, 'units', 'm') 147 iret = NF90_PUT_ATT(ntrajid, nwidthid, 'long_name', 'width') 148 iret = NF90_PUT_ATT(ntrajid, nwidthid, 'units', 'm') 149 iret = NF90_PUT_ATT(ntrajid, nlengthid, 'long_name', 'length') 150 iret = NF90_PUT_ATT(ntrajid, nlengthid, 'units', 'm') 151 iret = NF90_PUT_ATT(ntrajid, nyearid, 'long_name', 'calendar year') 152 iret = NF90_PUT_ATT(ntrajid, nyearid, 'units', 'years') 153 iret = NF90_PUT_ATT(ntrajid, ndayid, 'long_name', 'day of year') 154 iret = NF90_PUT_ATT(ntrajid, ndayid, 'units', 'days') 155 iret = NF90_PUT_ATT(ntrajid, nscaling_id, 'long_name', 'scaling factor for mass of berg') 156 iret = NF90_PUT_ATT(ntrajid, nscaling_id, 'units', 'none') 157 iret = NF90_PUT_ATT(ntrajid, nmass_of_bits_id, 'long_name', 'mass of bergy bits') 158 iret = NF90_PUT_ATT(ntrajid, nmass_of_bits_id, 'units', 'kg') 159 iret = NF90_PUT_ATT(ntrajid, nheat_density_id, 'long_name', 'heat density') 160 iret = NF90_PUT_ATT(ntrajid, nheat_density_id, 'units', 'J/kg') 161 161 162 162 ! End define mode 163 iret = NF90_ENDDEF( trajid)163 iret = NF90_ENDDEF(ntrajid) 164 164 165 165 END SUBROUTINE traj_init … … 173 173 INTEGER, INTENT( in ) :: kt 174 174 ! Local variables 175 INTEGER :: iret, k176 CHARACTER(len=80) :: filename175 INTEGER :: iret, jn 176 CHARACTER(len=80) :: cl_filename 177 177 TYPE(iceberg), POINTER :: this 178 178 TYPE(point) , POINTER :: pt … … 182 182 183 183 this => first_berg 184 k = kount_traj184 jn = num_traj 185 185 DO WHILE (ASSOCIATED(this)) 186 186 pt => this%current_point 187 k=k+1188 189 iret = NF90_PUT_VAR( trajid, numberid, this%number, (/1,k/), (/nkounts,1/) )190 iret = NF90_PUT_VAR( trajid, stepid, kt, (/ k/) )191 iret = NF90_PUT_VAR( trajid, scaling_id, this%mass_scaling, (/ k/) )192 193 iret = NF90_PUT_VAR( trajid, lonid, pt%lon, (/ k/) )194 iret = NF90_PUT_VAR( trajid, latid, pt%lat, (/ k/) )195 iret = NF90_PUT_VAR( trajid, xid, pt%xi, (/ k/) )196 iret = NF90_PUT_VAR( trajid, yid, pt%yj, (/ k/) )197 iret = NF90_PUT_VAR( trajid, uvelid, pt%uvel, (/ k/) )198 iret = NF90_PUT_VAR( trajid, vvelid, pt%vvel, (/ k/) )199 iret = NF90_PUT_VAR( trajid, uoid, pt%uo, (/ k/) )200 iret = NF90_PUT_VAR( trajid, void, pt%vo, (/ k/) )201 iret = NF90_PUT_VAR( trajid, uaid, pt%ua, (/ k/) )202 iret = NF90_PUT_VAR( trajid, vaid, pt%va, (/ k/) )203 iret = NF90_PUT_VAR( trajid, uiid, pt%ui, (/ k/) )204 iret = NF90_PUT_VAR( trajid, viid, pt%vi, (/ k/) )205 iret = NF90_PUT_VAR( trajid, sshxid, pt%ssh_x, (/ k/) )206 iret = NF90_PUT_VAR( trajid, sshyid, pt%ssh_y, (/ k/) )207 iret = NF90_PUT_VAR( trajid, sstid, pt%sst, (/ k/) )208 iret = NF90_PUT_VAR( trajid, cntid, pt%cn, (/ k/) )209 iret = NF90_PUT_VAR( trajid, thkid, pt%hi, (/ k/) )210 iret = NF90_PUT_VAR( trajid, massid, pt%mass, (/ k/) )211 iret = NF90_PUT_VAR( trajid, thicknessid, pt%thickness, (/ k/) )212 iret = NF90_PUT_VAR( trajid, widthid, pt%width, (/ k/) )213 iret = NF90_PUT_VAR( trajid, lengthid, pt%length, (/ k/) )214 iret = NF90_PUT_VAR( trajid, yearid, pt%year, (/ k/) )215 iret = NF90_PUT_VAR( trajid, dayid, pt%day, (/ k/) )216 iret = NF90_PUT_VAR( trajid, mass_of_bits_id, pt%mass_of_bits, (/ k/) )217 iret = NF90_PUT_VAR( trajid, heat_density_id, pt%heat_density, (/ k/) )187 jn=jn+1 188 189 iret = NF90_PUT_VAR(ntrajid, numberid, this%number, (/1,jn/), (/nkounts,1/) ) 190 iret = NF90_PUT_VAR(ntrajid, nstepid, kt, (/ jn /) ) 191 iret = NF90_PUT_VAR(ntrajid, nscaling_id, this%mass_scaling, (/ jn /) ) 192 193 iret = NF90_PUT_VAR(ntrajid, nlonid, pt%lon, (/ jn /) ) 194 iret = NF90_PUT_VAR(ntrajid, nlatid, pt%lat, (/ jn /) ) 195 iret = NF90_PUT_VAR(ntrajid, nxid, pt%xi, (/ jn /) ) 196 iret = NF90_PUT_VAR(ntrajid, nyid, pt%yj, (/ jn /) ) 197 iret = NF90_PUT_VAR(ntrajid, nuvelid, pt%uvel, (/ jn /) ) 198 iret = NF90_PUT_VAR(ntrajid, nvvelid, pt%vvel, (/ jn /) ) 199 iret = NF90_PUT_VAR(ntrajid, nuoid, pt%uo, (/ jn /) ) 200 iret = NF90_PUT_VAR(ntrajid, nvoid, pt%vo, (/ jn /) ) 201 iret = NF90_PUT_VAR(ntrajid, nuaid, pt%ua, (/ jn /) ) 202 iret = NF90_PUT_VAR(ntrajid, nvaid, pt%va, (/ jn /) ) 203 iret = NF90_PUT_VAR(ntrajid, nuiid, pt%ui, (/ jn /) ) 204 iret = NF90_PUT_VAR(ntrajid, nviid, pt%vi, (/ jn /) ) 205 iret = NF90_PUT_VAR(ntrajid, nsshxid, pt%ssh_x, (/ jn /) ) 206 iret = NF90_PUT_VAR(ntrajid, nsshyid, pt%ssh_y, (/ jn /) ) 207 iret = NF90_PUT_VAR(ntrajid, nsstid, pt%sst, (/ jn /) ) 208 iret = NF90_PUT_VAR(ntrajid, ncntid, pt%cn, (/ jn /) ) 209 iret = NF90_PUT_VAR(ntrajid, nthkid, pt%hi, (/ jn /) ) 210 iret = NF90_PUT_VAR(ntrajid, nmassid, pt%mass, (/ jn /) ) 211 iret = NF90_PUT_VAR(ntrajid, nthicknessid, pt%thickness, (/ jn /) ) 212 iret = NF90_PUT_VAR(ntrajid, nwidthid, pt%width, (/ jn /) ) 213 iret = NF90_PUT_VAR(ntrajid, nlengthid, pt%length, (/ jn /) ) 214 iret = NF90_PUT_VAR(ntrajid, nyearid, pt%year, (/ jn /) ) 215 iret = NF90_PUT_VAR(ntrajid, ndayid, pt%day, (/ jn /) ) 216 iret = NF90_PUT_VAR(ntrajid, nmass_of_bits_id, pt%mass_of_bits, (/ jn /) ) 217 iret = NF90_PUT_VAR(ntrajid, nheat_density_id, pt%heat_density, (/ jn /) ) 218 218 219 219 this=>this%next 220 220 ENDDO 221 IF( lwp .and. nn_verbose_level > 0 ) WRITE(numout,*) 'trajectory write to frame ', k222 kount_traj = k221 IF( lwp .and. nn_verbose_level > 0 ) WRITE(numout,*) 'trajectory write to frame ', jn 222 num_traj = jn 223 223 224 224 END SUBROUTINE traj_write … … 231 231 232 232 ! flush to file 233 iret = NF90_SYNC( trajid)233 iret = NF90_SYNC(ntrajid) 234 234 IF (iret .NE. NF90_NOERR) CALL ctl_stop('icebergs, traj_sync: nf_sync failed') 235 235 … … 243 243 244 244 ! Finish up 245 iret = NF90_CLOSE( trajid)245 iret = NF90_CLOSE(ntrajid) 246 246 IF (iret .NE. NF90_NOERR) CALL ctl_stop('icebergs, traj_end: nf_close failed') 247 247 -
branches/2012/dev_r3337_NOCS10_ICB/NEMOGCM/NEMO/OPA_SRC/ICB/icbutl.F90
r3339 r3359 34 34 PRIVATE 35 35 36 PUBLIC copy_flds ! routine called in xxx.F90module37 PUBLIC interp_flds ! routine called in xxx.F90 module38 PUBLIC bilin 39 PUBLIC bilin_x 40 P UBLICbilin_e41 PUBLIC add_new_berg_to_list 42 P UBLICinsert_berg_into_list43 PUBLIC delete_iceberg_from_list 44 PUBLIC destroy_iceberg 45 PUBLIC track_berg 46 PUBLIC print_berg 47 PUBLIC print_bergs 48 PUBLIC count_bergs 49 PUBLIC increment_kounter 50 PUBLIC yearday 51 PUBLIC sum_mass 52 PUBLIC sum_heat 36 PUBLIC copy_flds ! routine called in icbrun module 37 PUBLIC interp_flds ! routine called in icbdyn, icbthm modules 38 PUBLIC bilin ! routine called in icbini, icbdyn modules 39 PUBLIC bilin_x ! routine called in icbdyn module 40 PRIVATE bilin_e 41 PUBLIC add_new_berg_to_list ! routine called in icbini.F90, icbclv, icblbc and icbrst modules 42 PRIVATE insert_berg_into_list 43 PUBLIC delete_iceberg_from_list ! routine called in icblbc, icbthm modules 44 PUBLIC destroy_iceberg ! routine called in icbrun module 45 PUBLIC track_berg ! routine not currently used, retain just in case 46 PUBLIC print_berg ! routine called in icbthm module 47 PUBLIC print_bergs ! routine called in icbini, icbrun module 48 PUBLIC count_bergs ! routine called in icbdia, icbini, icblbc, icbrst modules 49 PUBLIC increment_kounter ! routine called in icbini, icbclv modules 50 PUBLIC yearday ! routine called in icbclv, icbrun module 51 PUBLIC sum_mass ! routine called in icbdia module 52 PUBLIC sum_heat ! routine called in icbdia module 53 53 54 54 PRIVATE create_iceberg … … 109 109 !!------------------------------------------------------------------------- 110 110 111 SUBROUTINE interp_flds( pi, pe1, puo, pui, pua, pssh_i, & 112 & pj, pe2, pvo, pvi, pva, pssh_j, psst, pcn, phi, pff ) 111 SUBROUTINE interp_flds( pi, pe1, puo, pui, pua, pssh_i, & 112 & pj, pe2, pvo, pvi, pva, pssh_j, & 113 & psst, pcn, phi, pff ) 113 114 !!---------------------------------------------------------------------- 114 115 !! *** ROUTINE interp_flds *** … … 168 169 !!------------------------------------------------------------------------- 169 170 170 REAL(wp) FUNCTION bilin( pfld, pi, pj, cd_type, jdi, jdj )171 REAL(wp) FUNCTION bilin( pfld, pi, pj, cd_type, kdi, kdj ) 171 172 !!---------------------------------------------------------------------- 172 173 !! *** FUNCTION bilin *** … … 178 179 !! 179 180 !!---------------------------------------------------------------------- 180 INTEGER , INTENT(in) :: jdi, jdj ! extra halo on grid181 REAL(wp), DIMENSION(1- jdi:jpi+jdi,1-jdj:jpj+jdj), INTENT(in) :: pfld ! field to be interpolated181 INTEGER , INTENT(in) :: kdi, kdj ! extra halo on grid 182 REAL(wp), DIMENSION(1-kdi:jpi+kdi,1-kdj:jpj+kdj), INTENT(in) :: pfld ! field to be interpolated 182 183 REAL(wp) , INTENT(in) :: pi, pj ! targeted coordinates in (i,j) referential 183 184 CHARACTER(len=1) , INTENT(in) :: cd_type ! type of pfld array grid-points: = T , U , V or F points … … 424 425 !!------------------------------------------------------------------------- 425 426 426 REAL(wp) FUNCTION yearday( imon, iday, ihr, imin, isec)427 REAL(wp) FUNCTION yearday(kmon, kday, khr, kmin, ksec) 427 428 ! sga - improved but still only applies to 365 day year, need to do this properly 428 429 ! Arguments 429 INTEGER, intent(in) :: imon, iday, ihr, imin, isec430 INTEGER, intent(in) :: kmon, kday, khr, kmin, ksec 430 431 ! Local variables 431 INTEGER :: i 432 INTEGER, DIMENSION(12) :: months = (/ 0,31,28,31,30,31,30,31,31,30,31,30 /) 433 434 yearday = FLOAT( SUM( months(1:imon) ) ) 435 yearday = yearday + FLOAT(iday-1) + (FLOAT(ihr) + (FLOAT(imin) + FLOAT(isec)/60.)/60.)/24. 432 INTEGER, DIMENSION(12) :: imonths = (/ 0,31,28,31,30,31,30,31,31,30,31,30 /) 433 434 yearday = FLOAT( SUM( imonths(1:kmon) ) ) 435 yearday = yearday + FLOAT(kday-1) + (FLOAT(khr) + (FLOAT(kmin) + FLOAT(ksec)/60.)/60.)/24. 436 436 437 437 END FUNCTION yearday … … 458 458 !!------------------------------------------------------------------------- 459 459 460 SUBROUTINE destroy_iceberg( berg)460 SUBROUTINE destroy_iceberg( berg ) 461 461 ! Arguments 462 462 TYPE(iceberg), POINTER :: berg … … 472 472 !!------------------------------------------------------------------------- 473 473 474 SUBROUTINE track_berg( num,label, kt )474 SUBROUTINE track_berg( knum, cd_label, kt ) 475 475 ! Arguments 476 INTEGER, DIMENSION(nkounts) :: num ! iceberg number477 CHARACTER(len=*) :: label476 INTEGER, DIMENSION(nkounts) :: knum ! iceberg number 477 CHARACTER(len=*) :: cd_label 478 478 INTEGER :: kt ! timestep number 479 479 ! Local variables … … 486 486 match = .TRUE. 487 487 DO k=1,nkounts 488 IF( this%number(k) /= num(k) ) match = .FALSE.488 IF( this%number(k) /= knum(k) ) match = .FALSE. 489 489 END DO 490 490 IF( match ) CALL print_berg(this, kt) … … 513 513 !!------------------------------------------------------------------------- 514 514 515 SUBROUTINE print_bergs( label, kt )515 SUBROUTINE print_bergs( cd_label, kt ) 516 516 ! Arguments 517 CHARACTER(len=*) :: label517 CHARACTER(len=*) :: cd_label 518 518 INTEGER :: kt ! timestep number 519 519 ! Local variables 520 INTEGER :: nbergs, nnbergs520 INTEGER :: ibergs, inbergs 521 521 TYPE(iceberg), POINTER :: this 522 522 523 523 this => first_berg 524 524 IF( ASSOCIATED(this) ) THEN 525 WRITE(numicb,'(a," pe=(",i3,")")' ) label, narea525 WRITE(numicb,'(a," pe=(",i3,")")' ) cd_label, narea 526 526 WRITE(numicb,'(a8,4x,a6,12x,a5,15x,a7,19x,a3,17x,a5,17x,a5,17x,a5)' ) & 527 527 'timestep', 'number', 'xi,yj','lon,lat','u,v','uo,vo','ua,va','ui,vi' … … 531 531 this => this%next 532 532 ENDDO 533 nbergs = count_bergs()534 nnbergs = nbergs535 IF( lk_mpp ) CALL mpp_sum( nnbergs)536 IF ( nbergs .GT. 0 ) WRITE(numicb,'(a," there are",i5," bergs out of",i6," on PE ",i4)') &537 label, nbergs, nnbergs, narea533 ibergs = count_bergs() 534 inbergs = ibergs 535 IF( lk_mpp ) CALL mpp_sum(inbergs) 536 IF ( ibergs .GT. 0 ) WRITE(numicb,'(a," there are",i5," bergs out of",i6," on PE ",i4)') & 537 cd_label, ibergs, inbergs, narea 538 538 539 539 END SUBROUTINE print_bergs … … 543 543 SUBROUTINE increment_kounter() 544 544 ! Small routine for coping with very large integer values labelling icebergs 545 ! kount_bergs is a array of integers545 ! num_bergs is a array of integers 546 546 ! the first member is incremented in steps of jpnij starting from narea 547 547 ! this means each iceberg is labelled with a unique number … … 550 550 551 551 ! Local variables 552 INTEGER :: i , ibig553 554 ibig = HUGE( kount_bergs(1))555 IF( ibig-jpnij < kount_bergs(1) ) THEN556 kount_bergs(1) = narea557 DO i = 2,nkounts558 IF( kount_bergs(i) == ibig ) THEN559 kount_bergs(i) = 0560 IF( i == nkounts ) CALL ctl_stop('Sorry, run out of iceberg number space')552 INTEGER :: ii, ibig 553 554 ibig = HUGE(num_bergs(1)) 555 IF( ibig-jpnij < num_bergs(1) ) THEN 556 num_bergs(1) = narea 557 DO ii = 2,nkounts 558 IF( num_bergs(ii) == ibig ) THEN 559 num_bergs(ii) = 0 560 IF( ii == nkounts ) CALL ctl_stop('Sorry, run out of iceberg number space') 561 561 ELSE 562 kount_bergs(i) = kount_bergs(i) + 1562 num_bergs(ii) = num_bergs(ii) + 1 563 563 EXIT 564 564 ENDIF 565 565 END DO 566 566 ELSE 567 kount_bergs(1) = kount_bergs(1) + jpnij567 num_bergs(1) = num_bergs(1) + jpnij 568 568 ENDIF 569 569
Note: See TracChangeset
for help on using the changeset viewer.
