Changeset 1922 for trunk/NEMO/LIM_SRC_2
- Timestamp:
- 2010-06-09T10:24:58+02:00 (14 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/NEMO/LIM_SRC_2/limtrp_2.F90
r1715 r1922 4 4 !! LIM 2.0 transport ice model : sea-ice advection/diffusion 5 5 !!====================================================================== 6 !! History : LIM ! 2000-01 (UCL) Original code 7 !! 2.0 ! 2001-05 (G. Madec, R. Hordoir) opa norm 8 !! - ! 2004-01 (G. Madec, C. Ethe) F90, mpp 9 !!---------------------------------------------------------------------- 6 10 #if defined key_lim2 7 11 !!---------------------------------------------------------------------- … … 11 15 !! lim_trp_init_2 : initialization and namelist read 12 16 !!---------------------------------------------------------------------- 13 !! * Modules used14 USE phycst15 USE dom_oce 17 USE phycst ! physical constant 18 USE sbc_oce ! ocean surface boundary condition 19 USE dom_oce ! ocean domain 16 20 USE in_out_manager ! I/O manager 17 USE dom_ice_2 18 USE ice_2 19 USE limistate_2 20 USE limadv_2 21 USE limhdf_2 22 USE lbclnk 23 USE lib_mpp 21 USE dom_ice_2 ! LIM-2 domain 22 USE ice_2 ! LIM-2 variables 23 USE limistate_2 ! LIM-2 initial state 24 USE limadv_2 ! LIM-2 advection 25 USE limhdf_2 ! LIM-2 horizontal diffusion 26 USE lbclnk ! lateral boundary conditions -- MPP exchanges 27 USE lib_mpp ! MPP library 24 28 25 29 IMPLICIT NONE 26 30 PRIVATE 27 31 28 !! * Routine accessibility 29 PUBLIC lim_trp_2 ! called by sbc_ice_lim_2 30 31 !! * Shared module variables 32 REAL(wp), PUBLIC :: & !: 33 bound = 0.e0 !: boundary condit. (0.0 no-slip, 1.0 free-slip) 34 35 !! * Module variables 32 PUBLIC lim_trp_2 ! called by sbc_ice_lim_2 33 34 REAL(wp), PUBLIC :: bound = 0.e0 !: boundary condit. (0.0 no-slip, 1.0 free-slip) 35 36 36 REAL(wp) :: & ! constant values 37 37 epsi06 = 1.e-06 , & … … 44 44 # include "vectopt_loop_substitute.h90" 45 45 !!---------------------------------------------------------------------- 46 !! LIM 2.0, UCL-LOCEAN-IPSL (2005)46 !! NEMO/LIM 3.2, UCL-LOCEAN-IPSL (2010) 47 47 !! $Id$ 48 !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt48 !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 49 49 !!---------------------------------------------------------------------- 50 50 … … 62 62 !! 63 63 !! ** action : 64 !!65 !! History :66 !! 1.0 ! 00-01 (LIM) Original code67 !! ! 01-05 (G. Madec, R. Hordoir) opa norm68 !! 2.0 ! 04-01 (G. Madec, C. Ethe) F90, mpp69 64 !!--------------------------------------------------------------------- 70 65 INTEGER, INTENT(in) :: kt ! number of iteration 71 72 INTEGER :: ji, jj, jk, & ! dummy loop indices 73 & initad ! number of sub-timestep for the advection 74 75 REAL(wp) :: & 76 zindb , & 77 zacrith, & 78 zindsn , & 79 zindic , & 80 zusvosn, & 81 zusvoic, & 82 zignm , & 83 zindhe , & 84 zvbord , & 85 zcfl , & 86 zusnit , & 87 zrtt, ztsn, ztic1, ztic2 88 89 REAL(wp), DIMENSION(jpi,jpj) :: & ! temporary workspace 90 zui_u , zvi_v , zsm , & 91 zs0ice, zs0sn , zs0a , & 92 zs0c0 , zs0c1 , zs0c2 , & 93 zs0st 66 !! 67 INTEGER :: ji, jj, jk ! dummy loop indices 68 INTEGER :: initad ! number of sub-timestep for the advection 69 REAL(wp) :: zindb , zindsn , zindic, zacrith ! local scalars 70 REAL(wp) :: zusvosn, zusvoic, zignm , zindhe ! - - 71 REAL(wp) :: zvbord , zcfl , zusnit ! - - 72 REAL(wp) :: zrtt , ztsn , ztic1 , ztic2 ! - - 73 REAL(wp), DIMENSION(jpi,jpj) :: zui_u , zvi_v , zsm ! 2D workspace 74 REAL(wp), DIMENSION(jpi,jpj) :: zs0ice, zs0sn , zs0a ! - - 75 REAL(wp), DIMENSION(jpi,jpj) :: zs0c0 , zs0c1 , zs0c2 , zs0st ! - - 94 76 !--------------------------------------------------------------------- 95 77 … … 105 87 ! ice velocities at ocean U- and V-points (zui_u,zvi_v) 106 88 ! --------------------------------------- 107 ! zvbord factor between 1 and 2 to take into account slip or no-slip boundary conditions. 108 zvbord = 1.0 + ( 1.0 - bound ) 89 zvbord = 1.0 + ( 1.0 - bound ) ! zvbord=2 no-slip, =0 free slip boundary conditions 109 90 DO jj = 1, jpjm1 110 91 DO ji = 1, jpim1 ! NO vector opt. … … 113 94 END DO 114 95 END DO 115 ! Lateral boundary conditions on zui_u, zvi_v 116 CALL lbc_lnk( zui_u, 'U', -1. ) 117 CALL lbc_lnk( zvi_v, 'V', -1. ) 96 CALL lbc_lnk( zui_u, 'U', -1. ) ; CALL lbc_lnk( zvi_v, 'V', -1. ) ! Lateral boundary conditions 97 118 98 119 99 ! CFL test for stability … … 122 102 zcfl = MAX( zcfl, MAXVAL( ABS( zui_u(1:jpim1, : ) ) * rdt_ice / e1u(1:jpim1, : ) ) ) 123 103 zcfl = MAX( zcfl, MAXVAL( ABS( zvi_v( : ,1:jpjm1) ) * rdt_ice / e2v( : ,1:jpjm1) ) ) 124 125 IF (lk_mpp ) CALL mpp_max(zcfl)126 127 IF ( zcfl > 0.5 .AND. lwp ) WRITE(numout,*) 'lim_trp_2 : violation of cfl criterion the ',nday,'th day, cfl = ',zcfl104 ! 105 IF(lk_mpp) CALL mpp_max( zcfl ) 106 ! 107 IF( zcfl > 0.5 .AND. lwp ) WRITE(numout,*) 'lim_trp_2 : violation of cfl criterion the ',nday,'th day, cfl = ', zcfl 128 108 129 109 ! content of properties 130 110 ! --------------------- 131 111 zs0sn (:,:) = hsnm(:,:) * area(:,:) ! Snow volume. 132 zs0ice(:,:) = hicm (:,:) * area(:,:)! Ice volume.133 zs0a (:,:) = ( 1.0 - frld(:,:) ) * area(:,:)! Surface covered by ice.134 zs0c0 (:,:) = tbif(:,:,1) / rt0_snow * zs0sn (:,:)! Heat content of the snow layer.112 zs0ice(:,:) = hicm(:,:) * area(:,:) ! Ice volume. 113 zs0a (:,:) = ( 1.0 - frld(:,:) ) * area (:,:) ! Surface covered by ice. 114 zs0c0 (:,:) = tbif(:,:,1) / rt0_snow * zs0sn (:,:) ! Heat content of the snow layer. 135 115 zs0c1 (:,:) = tbif(:,:,2) / rt0_ice * zs0ice(:,:) ! Heat content of the first ice layer. 136 116 zs0c2 (:,:) = tbif(:,:,3) / rt0_ice * zs0ice(:,:) ! Heat content of the second ice layer. 137 zs0st (:,:) = qstoif(:,:) / xlic * zs0a (:,:)! Heat reservoir for brine pockets.117 zs0st (:,:) = qstoif(:,:) / xlic * zs0a (:,:) ! Heat reservoir for brine pockets. 138 118 139 119 140 ! Advection 120 ! Advection (Prather scheme) 141 121 ! --------- 142 ! If ice drift field is too fast, use an appropriate time step for advection. 143 initad = 1 + INT( MAX( rzero, SIGN( rone, zcfl-0.5 ) ) ) 144 zusnit = 1.0 / REAL( initad ) 145 146 IF ( MOD( nday , 2 ) == 0) THEN 147 DO jk = 1,initad 122 initad = 1 + INT( MAX( rzero, SIGN( rone, zcfl-0.5 ) ) ) ! If ice drift field is too fast, 123 zusnit = 1.0 / REAL( initad ) ! split the ice time step in two 124 ! 125 IF( MOD( ( kt - 1) / nn_fsbc , 2 ) == 0) THEN !== odd ice time step: adv_x then adv_y ==! 126 DO jk = 1, initad 148 127 CALL lim_adv_x_2( zusnit, zui_u, rone , zsm, zs0ice, sxice, sxxice, syice, syyice, sxyice ) 149 128 CALL lim_adv_y_2( zusnit, zvi_v, rzero, zsm, zs0ice, sxice, sxxice, syice, syyice, sxyice ) … … 161 140 CALL lim_adv_y_2( zusnit, zvi_v, rzero, zsm, zs0st , sxst , sxxst , syst , syyst , sxyst ) 162 141 END DO 163 ELSE 142 ELSE !== even ice time step: adv_x then adv_y ==! 164 143 DO jk = 1, initad 165 144 CALL lim_adv_y_2( zusnit, zvi_v, rone , zsm, zs0ice, sxice, sxxice, syice, syyice, sxyice ) … … 182 161 ! recover the properties from their contents 183 162 ! ------------------------------------------ 163 !!gm Define in limmsh one for all area = 1 /area (CPU time saved !) 184 164 zs0ice(:,:) = zs0ice(:,:) / area(:,:) 185 165 zs0sn (:,:) = zs0sn (:,:) / area(:,:) … … 205 185 END DO 206 186 END DO 187 !!gm more readable coding: (and avoid an error in F90 with sign of zero) 188 ! DO jj = 1, jpjm1 ! NB: has not to be defined on jpj line and jpi row 189 ! DO ji = 1 , fs_jpim1 ! vector opt. 190 ! IF( MIN( zs0a(ji,jj) , zs0a(ji+1,jj) ) == 0.e0 ) pahu(ji,jj) = 0.e0 191 ! IF( MIN( zs0a(ji,jj) , zs0a(ji,jj+1) ) == 0.e0 ) pahv(ji,jj) = 0.e0 192 ! END DO 193 ! END DO 194 !!gm end 207 195 208 196 ! diffusion … … 216 204 CALL lim_hdf_2( zs0st ) 217 205 218 zs0ice(:,:) = MAX( rzero, zs0ice(:,:) * area(:,:) ) !!bug: est-ce utile 219 zs0sn (:,:) = MAX( rzero, zs0sn (:,:) * area(:,:) ) !!bug: cf /area juste apres 220 zs0a (:,:) = MAX( rzero, zs0a (:,:) * area(:,:) ) !! suppression des 2 change le resultat... 221 zs0c0 (:,:) = MAX( rzero, zs0c0 (:,:) * area(:,:) ) 206 !!gm see comment this can be skipped 207 zs0ice(:,:) = MAX( rzero, zs0ice(:,:) * area(:,:) ) !!bug: useless 208 zs0sn (:,:) = MAX( rzero, zs0sn (:,:) * area(:,:) ) !!bug: cf /area just below 209 zs0a (:,:) = MAX( rzero, zs0a (:,:) * area(:,:) ) !! caution: the suppression of the 2 changes 210 zs0c0 (:,:) = MAX( rzero, zs0c0 (:,:) * area(:,:) ) !! the last digit of the results 222 211 zs0c1 (:,:) = MAX( rzero, zs0c1 (:,:) * area(:,:) ) 223 212 zs0c2 (:,:) = MAX( rzero, zs0c2 (:,:) * area(:,:) ) … … 225 214 226 215 227 ! -------------------------------------------------------------------! 228 ! Up-dating and limitation of sea ice properties after transport ! 229 ! -------------------------------------------------------------------! 230 231 ! Up-dating and limitation of sea ice properties after transport. 216 !-------------------------------------------------------------------! 217 ! Updating and limitation of sea ice properties after transport ! 218 !-------------------------------------------------------------------! 232 219 DO jj = 1, jpj 233 !!!iii zindhe = REAL( MAX( 0, isign(1, jj - njeq ) ) ) !ibug mpp !!bugmpp njeq!234 220 zindhe = MAX( 0.e0, SIGN( 1.e0, fcor(1,jj) ) ) ! = 0 for SH, =1 for NH 235 221 DO ji = 1, jpi 236 222 ! 237 223 ! Recover mean values over the grid squares. 238 224 zs0sn (ji,jj) = MAX( rzero, zs0sn (ji,jj)/area(ji,jj) ) … … 272 258 END DO 273 259 END DO 274 260 ! 275 261 ENDIF 276 262 ! 277 263 END SUBROUTINE lim_trp_2 278 264 … … 284 270 !! ** Purpose : initialization of ice advection parameters 285 271 !! 286 !! ** Method : Read the namicetrp namelist and check the parameter287 !! values called at the first timestep (nit000)272 !! ** Method : Read the namicetrp namelist and check the parameter 273 !! values called at the first timestep (nit000) 288 274 !! 289 275 !! ** input : Namelist namicetrp 290 !!291 !! history :292 !! 2.0 ! 03-08 (C. Ethe) Original code293 276 !!------------------------------------------------------------------- 294 277 NAMELIST/namicetrp/ bound 295 278 !!------------------------------------------------------------------- 296 297 ! Read Namelist namicetrp 298 REWIND ( numnam_ice ) 279 ! 280 REWIND ( numnam_ice ) ! Read Namelist namicetrp 299 281 READ ( numnam_ice , namicetrp ) 300 282 IF(lwp) THEN … … 304 286 WRITE(numout,*) ' boundary conditions (0. no-slip, 1. free-slip) bound = ', bound 305 287 ENDIF 306 288 ! 307 289 END SUBROUTINE lim_trp_init_2 308 290
Note: See TracChangeset
for help on using the changeset viewer.