Changeset 77 for trunk/NEMO/LIM_SRC
- Timestamp:
- 2004-04-22T14:49:55+02:00 (20 years ago)
- Location:
- trunk/NEMO/LIM_SRC
- Files:
-
- 4 edited
-
dom_ice.F90 (modified) (1 diff)
-
limdia.F90 (modified) (12 diffs)
-
limistate.F90 (modified) (10 diffs)
-
limrhg.F90 (modified) (31 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/NEMO/LIM_SRC/dom_ice.F90
r12 r77 13 13 14 14 IMPLICIT NONE 15 P UBLIC15 PRIVATE 16 16 17 REAL(wp) :: & !: 18 dz !: depht of the 1st ocean level 17 !! * Share module variables 18 LOGICAL, PUBLIC :: & !: 19 l_jeq = .TRUE. !: Equator inside the domain flag 19 20 20 REAL(wp),DIMENSION(jpi,jpj) :: & !: 21 area !: Surface of a grid square 21 INTEGER, PUBLIC :: & !: 22 njeq , njeqm1 !: j-index of the equator if it is inside the domain 23 ! ! (otherwise = jpj+10 (SH) or -10 (SH) ) 22 24 23 REAL(wp),DIMENSION(jpi,jpj,2,2) :: & !: 24 wght , & !: weight of the 4 neighbours to compute averages 25 akappa , & !: first group of metric coefficients 26 bkappa !: third group of metric coefficients 25 REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: & !: 26 fs2cor , & !: coriolis factor 27 fcor , & !: coriolis coefficient 28 covrai , & !: sine of geographic latitude 29 area , & !: surface of grid cell 30 tms , tmu !: temperature and velocity points masks 27 31 28 REAL(wp),DIMENSION(jpi,jpj) :: & !: 29 tms , tmu !: masks arrays 32 REAL(wp), PUBLIC, DIMENSION(jpi,jpj,2,2) :: & !: 33 wght , & !: weight of the 4 neighbours to compute averages 34 akappa , & !: first group of metric coefficients 35 bkappa !: third group of metric coefficients 30 36 31 REAL(wp),DIMENSION(jpi,jpj,2,2,2,2) :: & !: 32 alambd !: second group of metric coefficients 33 34 REAL(wp),DIMENSION(jpi,jpj) :: & !: 35 fs2cor , & !: coriolis factor 36 fcor , & !: coriolis coefficient 37 covrai , & !: sine of geographic latitude 38 aire !: total area 39 40 INTEGER :: & !: 41 jeq , jeqm1 !: 37 REAL(wp), PUBLIC, DIMENSION(jpi,jpj,2,2,2,2) :: & !: 38 alambd !: second group of metric coefficients 42 39 43 40 !!====================================================================== -
trunk/NEMO/LIM_SRC/limdia.F90
r12 r77 12 12 !!---------------------------------------------------------------------- 13 13 !! * Modules used 14 USE phycst 15 USE in_out_manager14 USE phycst ! 15 USE par_ice ! ice parameters 16 16 USE ice_oce ! ice variables 17 USE daymod 18 USE dom_ice 19 USE ice 20 USE iceini 21 USE limistate 17 USE daymod ! 18 USE dom_ice ! 19 USE ice ! 20 USE iceini ! 21 USE limistate ! 22 USE in_out_manager ! I/O manager 22 23 23 24 IMPLICIT NONE … … 28 29 29 30 !! * Shared module variables 30 INTEGER, PUBLIC :: & 31 INTEGER, PUBLIC :: & !: 31 32 ntmoy = 1 , & !: instantaneous values of ice evolution or averaging ntmoy 32 33 ninfo = 1 !: frequency of ouputs on file ice_evolu in case of averaging 33 34 34 35 !! * Module variables 35 INTEGER :: & 36 nfrinf = 4 ! number of variables written in one line 37 38 INTEGER :: & 36 INTEGER, PARAMETER :: & ! Parameters for outputs to files "evolu" 37 jpinfmx = 100 , & ! maximum number of key variables 38 jpchinf = 5 , & ! ??? 39 jpchsep = jpchinf + 2 ! ??? 40 41 INTEGER :: & 42 nfrinf = 4 , & ! number of variables written in one line 39 43 nferme , & ! last time step at which the var. are written on file 40 44 nvinfo , & ! number of total variables 41 45 nbvt , & ! number of time variables 42 46 naveg ! number of step for accumulation before averaging 43 44 REAL(wp), DIMENSION(ninfmx) :: & 45 vinfom ! temporary working space 46 47 CHARACTER(len=8) :: & 47 48 CHARACTER(len=8) :: & 48 49 fmtinf = '1PE13.5 ' ! format of the output values 49 50 CHARACTER(len=30) :: & 50 CHARACTER(len=30) :: & 51 51 fmtw , & ! formats 52 52 fmtr , & ! ??? 53 53 fmtitr ! ??? 54 55 CHARACTER(len=nchsep), DIMENSION(ninfmx) :: & 54 CHARACTER(len=jpchsep), DIMENSION(jpinfmx) :: & 56 55 titvar ! title of key variables 57 58 REAL(wp) :: & 59 epsi06 = 1.e-06 56 57 REAL(wp) :: & 58 epsi06 = 1.e-06 ! ??? 59 REAL(wp), DIMENSION(jpinfmx) :: & 60 vinfom ! temporary working space 61 REAL(wp), DIMENSION(jpi,jpj) :: & 62 aire ! masked grid cell area 60 63 61 64 !! * Substitutions … … 81 84 INTEGER :: jv,ji, jj ! dummy loop indices 82 85 INTEGER :: nv ! indice of variable 83 REAL(wp), DIMENSION( ninfmx) :: &86 REAL(wp), DIMENSION(jpinfmx) :: & 84 87 vinfor ! temporary working space 85 88 REAL(wp) :: & … … 98 101 99 102 nv = 1 100 vinfor(nv) = REAL( numit)103 vinfor(nv) = REAL( numit ) 101 104 nv = nv + 1 102 105 vinfor(nv) = nyear 103 106 104 107 DO jv = nbvt + 1, nvinfo 105 vinfor(jv) = 0. 0108 vinfor(jv) = 0.e0 106 109 END DO 107 110 … … 109 112 zextent85 = 0.e0 110 113 ! variables in northern Hemis 111 DO jj = jeq, jpjm1114 DO jj = njeq, jpjm1 112 115 DO ji = fs_2, fs_jpim1 ! vector opt. 113 116 IF( tms(ji,jj) == 1 ) THEN … … 135 138 ! variables in southern Hemis 136 139 nv = nv + 1 137 DO jj = 2, jeqm1140 DO jj = 2, njeqm1 138 141 DO ji = fs_2, fs_jpim1 ! vector opt. 139 142 IF( tms(ji,jj) == 1 ) THEN … … 169 172 naveg = 0 170 173 DO jv = 1, nvinfo 171 vinfom(jv) =0.0174 vinfom(jv) = 0.e0 172 175 END DO 173 176 ENDIF … … 199 202 REAL(wp) :: zxx0, zxx1 ! temporary scalars 200 203 201 CHARACTER(len= nchinf) :: titinf204 CHARACTER(len=jpchinf) :: titinf 202 205 !!------------------------------------------------------------------- 203 204 206 205 207 ! Read Namelist namicedia … … 216 218 ENDIF 217 219 220 ! masked grid cell area 221 aire(:,:) = area(:,:) * tms(:,:) 222 218 223 ! Titles of ice key variables : 219 224 nv = 1 … … 265 270 266 271 ! definition of formats 267 WRITE( fmtw , '(A,I3,A2,I1,A)' ) '(', nfrinf, '(A', nchsep, ','//fmtinf//'))'268 WRITE( fmtr , '(A,I3,A,I1,A)' ) '(', nfrinf, '(', nchsep, 'X,'//fmtinf//'))'269 WRITE( fmtitr, '(A,I3,A,I1,A)' ) '(', nvinfo, 'A', nchinf, ')'272 WRITE( fmtw , '(A,I3,A2,I1,A)' ) '(', nfrinf, '(A', jpchsep, ','//fmtinf//'))' 273 WRITE( fmtr , '(A,I3,A,I1,A)' ) '(', nfrinf, '(', jpchsep, 'X,'//fmtinf//'))' 274 WRITE( fmtitr, '(A,I3,A,I1,A)' ) '(', nvinfo, 'A', jpchinf, ')' 270 275 271 276 ! opening "ice_evolu" file 272 irecl = ( nchinf + 1 ) * nvinfo277 irecl = ( jpchinf + 1 ) * nvinfo 273 278 OPEN( numevo_ice, file='ice.evolu', status='unknown', RECL = irecl) 274 279 OPEN( numevo_ice, file='ice.evolu', status='unknown') … … 276 281 !- ecriture de 2 lignes d''entete : 277 282 WRITE(numevo_ice,1000) fmtr, fmtw, fmtitr, nvinfo, ntot, 0, nfrinf 278 zxx0 = 0.001 * REAL( ninfo)279 zxx1 = 0.001 * REAL( ndeb)280 WRITE(numevo_ice,1111) REAL( nchinf), 0., zxx1, zxx0, 0., 0., 0283 zxx0 = 0.001 * REAL( ninfo ) 284 zxx1 = 0.001 * REAL( ndeb ) 285 WRITE(numevo_ice,1111) REAL(jpchinf), 0., zxx1, zxx0, 0., 0., 0 281 286 282 287 !- ecriture de 2 lignes de titre : … … 289 294 !--preparation de "titvar" pour l''ecriture parmi les valeurs numeriques : 290 295 DO jv = 2 , nvinfo 291 titinf = titvar(jv)(: nchinf)296 titinf = titvar(jv)(:jpchinf) 292 297 titvar(jv) = ' '//titinf 293 298 END DO -
trunk/NEMO/LIM_SRC/limistate.F90
r3 r77 16 16 USE oce ! dynamics and tracers variables 17 17 USE dom_oce 18 USE par_ice ! ice parameters 18 19 USE ice_oce ! ice variables 19 20 USE in_out_manager … … 38 39 alins = 0.1 ! initial leads area in the south 39 40 40 REAL(wp) :: 41 zzero = 0. 0, &42 zone = 1. 041 REAL(wp) :: & ! constant values 42 zzero = 0.e0 , & 43 zone = 1.e0 43 44 !!---------------------------------------------------------------------- 44 45 !! LIM 2.0 , UCL-LODYC-IPSL (2003) … … 50 51 !!------------------------------------------------------------------- 51 52 !! *** ROUTINE lim_istate *** 52 !! Purpose : 53 !! restart from a state defined in a binary file and arbitrary sea 54 !! ice conditions 53 !! 54 !! ** Purpose : defined the sea-ice initial state 55 !! 56 !! ** Method : restart from a state defined in a binary file 57 !! or from arbitrary sea-ice conditions 55 58 !! 56 59 !! History : … … 59 62 !! * Local variables 60 63 INTEGER :: ji, jj, jk ! dummy loop indices 61 62 64 REAL(wp) :: zidto, & ! temporary scalar 63 65 zs0, ztf, zbin … … 82 84 83 85 84 u_io (:,:) = 0. 85 v_io (:,:) = 0. 86 u_io (:,:) = 0.e0 87 v_io (:,:) = 0.e0 86 88 sst_io(:,:) = ( nfice - 1 ) * ( tn(:,:,1) + rt0 ) ! use the ocean initial values 87 89 sss_io(:,:) = ( nfice - 1 ) * sn(:,:,1) ! tricky trick *(nfice-1) ! … … 96 98 tfu(:,:) = ztf 97 99 98 !-- Northern hemisphere. 99 DO jj = jeq , jpj 100 DO ji = 1 , jpi 100 DO jj = 1, jpj 101 DO ji = 1, jpi 101 102 !--- Criterion for presence (zidto=1) or absence (zidto=0) of ice 102 103 zidto = tms(ji,jj) * ( 1.0 - MAX(zzero, SIGN( zone, ztn(ji,jj) - tfu(ji,jj) - ttest) ) ) 103 104 104 hicif(ji,jj) = zidto * hginn 105 frld(ji,jj) = zidto * alinn + ( 1.0 - zidto ) * 1.0 106 hsnif(ji,jj) = zidto * hninn 107 ENDDO 108 ENDDO 109 110 !--- Southern hemisphere. 111 DO jj = 1, jeqm1 112 DO ji = 1, jpi 113 !--- Criterion for presence (zidto=1) or absence (zidto=0) of ice. 114 zidto = tms(ji,jj) * ( 1.0 - MAX(zzero, SIGN( zone, ztn(ji,jj) - tfu(ji,jj) - ttest) ) ) 115 116 hicif(ji,jj) = zidto * hgins 117 frld(ji,jj) = zidto * alins + ( 1.0 - zidto ) * 1.0 118 hsnif(ji,jj) = zidto * hnins 119 ENDDO 120 ENDDO 105 IF( fcor(ji,jj) >= 0.e0 ) THEN !-- Northern hemisphere. 106 hicif(ji,jj) = zidto * hginn 107 frld(ji,jj) = zidto * alinn + ( 1.0 - zidto ) * 1.0 108 hsnif(ji,jj) = zidto * hninn 109 ELSE !--- Southern hemisphere. 110 hicif(ji,jj) = zidto * hgins 111 frld(ji,jj) = zidto * alins + ( 1.0 - zidto ) * 1.0 112 hsnif(ji,jj) = zidto * hnins 113 ENDIF 114 END DO 115 END DO 121 116 122 117 sist (:,:) = tfu(:,:) … … 132 127 # endif 133 128 134 135 !--- Moments for advection. 136 137 sxice (:,:) = 0.e0 ; sxsn (:,:) = 0.e0 ; sxa (:,:) = 0.e0 138 syice (:,:) = 0.e0 ; sysn (:,:) = 0.e0 ; sya (:,:) = 0.e0 139 sxxice(:,:) = 0.e0 ; sxxsn(:,:) = 0.e0 ; sxxa (:,:) = 0.e0 140 syyice(:,:) = 0.e0 ; syysn(:,:) = 0.e0 ; syya (:,:) = 0.e0 141 sxyice(:,:) = 0.e0 ; sxysn(:,:) = 0.e0 ; sxya (:,:) = 0.e0 142 143 sxc0 (:,:) = 0.e0 ; sxc1 (:,:) = 0.e0 ; sxc2 (:,:) = 0.e0 144 syc0 (:,:) = 0.e0 ; syc1 (:,:) = 0.e0 ; syc2 (:,:) = 0.e0 145 sxxc0 (:,:) = 0.e0 ; sxxc1(:,:) = 0.e0 ; sxxc2(:,:) = 0.e0 146 syyc0 (:,:) = 0.e0 ; syyc1(:,:) = 0.e0 ; syyc2(:,:) = 0.e0 147 sxyc0 (:,:) = 0.e0 ; sxyc1(:,:) = 0.e0 ; sxyc2(:,:) = 0.e0 148 149 sxst (:,:) = 0.e0 150 syst (:,:) = 0.e0 151 sxxst (:,:) = 0.e0 152 syyst (:,:) = 0.e0 153 sxyst (:,:) = 0.e0 154 155 !-- lateral boundary conditions 156 CALL lbc_lnk( hicif, 'T', 1. ) 157 CALL lbc_lnk( frld , 'T', 1. ) 129 !--- Moments for advection. 130 131 sxice (:,:) = 0.e0 ; sxsn (:,:) = 0.e0 ; sxa (:,:) = 0.e0 132 syice (:,:) = 0.e0 ; sysn (:,:) = 0.e0 ; sya (:,:) = 0.e0 133 sxxice(:,:) = 0.e0 ; sxxsn(:,:) = 0.e0 ; sxxa (:,:) = 0.e0 134 syyice(:,:) = 0.e0 ; syysn(:,:) = 0.e0 ; syya (:,:) = 0.e0 135 sxyice(:,:) = 0.e0 ; sxysn(:,:) = 0.e0 ; sxya (:,:) = 0.e0 136 137 sxc0 (:,:) = 0.e0 ; sxc1 (:,:) = 0.e0 ; sxc2 (:,:) = 0.e0 138 syc0 (:,:) = 0.e0 ; syc1 (:,:) = 0.e0 ; syc2 (:,:) = 0.e0 139 sxxc0 (:,:) = 0.e0 ; sxxc1(:,:) = 0.e0 ; sxxc2(:,:) = 0.e0 140 syyc0 (:,:) = 0.e0 ; syyc1(:,:) = 0.e0 ; syyc2(:,:) = 0.e0 141 sxyc0 (:,:) = 0.e0 ; sxyc1(:,:) = 0.e0 ; sxyc2(:,:) = 0.e0 142 143 sxst (:,:) = 0.e0 144 syst (:,:) = 0.e0 145 sxxst (:,:) = 0.e0 146 syyst (:,:) = 0.e0 147 sxyst (:,:) = 0.e0 148 149 !-- lateral boundary conditions 150 CALL lbc_lnk( hicif, 'T', 1. ) 151 CALL lbc_lnk( frld , 'T', 1. ) 158 152 !i bug frld = 1 over land 159 153 frld(:,:) = tms(:,:) * frld(:,:) + ( 1. - tms(:,:) ) ! put 1 over land 160 154 !i 161 162 163 DO jk = 1, nlayersp1164 165 166 167 168 155 CALL lbc_lnk( hsnif, 'T', 1. ) 156 CALL lbc_lnk( sist , 'T', 1. ) 157 DO jk = 1, jplayersp1 158 CALL lbc_lnk(tbif(:,:,jk), 'T', 1. ) 159 END DO 160 CALL lbc_lnk( fsbbq , 'T', 1. ) 161 CALL lbc_lnk( qstoif , 'T', 1. ) 162 CALL lbc_lnk( sss_io , 'T', 1. ) 169 163 170 164 END SUBROUTINE lim_istate … … 175 169 !! *** ROUTINE lim_istate_init *** 176 170 !! 177 !! ** Purpose : Definition of initial state of the ice178 !! 179 !! ** Method :Read the namiceini namelist and check the parameter171 !! ** Purpose : Definition of initial state of the ice 172 !! 173 !! ** Method : Read the namiceini namelist and check the parameter 180 174 !! values called at the first timestep (nit000) 181 175 !! 182 !! ** input : 183 !! Namelist namiceini 184 !! 185 !! history : 176 !! ** input : Namelist namiceini 177 !! 178 !! history 186 179 !! 8.5 ! 03-08 (C. Ethe) original code 187 180 !!------------------------------------------------------------------- … … 189 182 !!------------------------------------------------------------------- 190 183 191 ! Define the initial parameters192 ! -------------------------193 194 184 ! Read Namelist namiceini 185 195 186 REWIND ( numnam_ice ) 196 187 READ ( numnam_ice , namiceini ) … … 199 190 WRITE(numout,*) 'lim_istate_init : ice parameters inititialisation ' 200 191 WRITE(numout,*) '~~~~~~~~~~~~~~~' 201 WRITE(numout,*) ' threshold water temp. for initial sea-ice ttest = ', ttest202 WRITE(numout,*) ' initial snow thickness in the north hninn = ', hninn203 WRITE(numout,*) ' initial ice thickness in the north hginn = ', hginn204 WRITE(numout,*) ' initial leads area in the north alinn = ', alinn205 WRITE(numout,*) ' initial snow thickness in the south hnins = ', hnins206 WRITE(numout,*) ' initial ice thickness in the south hgins = ', hgins207 WRITE(numout,*) ' initial leads area in the south alins = ', alins192 WRITE(numout,*) ' threshold water temp. for initial sea-ice ttest = ', ttest 193 WRITE(numout,*) ' initial snow thickness in the north hninn = ', hninn 194 WRITE(numout,*) ' initial ice thickness in the north hginn = ', hginn 195 WRITE(numout,*) ' initial leads area in the north alinn = ', alinn 196 WRITE(numout,*) ' initial snow thickness in the south hnins = ', hnins 197 WRITE(numout,*) ' initial ice thickness in the south hgins = ', hgins 198 WRITE(numout,*) ' initial leads area in the south alins = ', alins 208 199 ENDIF 209 200 -
trunk/NEMO/LIM_SRC/limrhg.F90
r12 r77 1 1 MODULE limrhg 2 #if defined key_ice_lim3 2 !!====================================================================== 4 3 !! *** MODULE limrhg *** 5 4 !! Ice rheology : performs sea ice rheology 6 5 !!====================================================================== 7 6 #if defined key_ice_lim 7 !!---------------------------------------------------------------------- 8 !! 'key_ice_lim' LIM sea-ice model 8 9 !!---------------------------------------------------------------------- 9 10 !! lim_rhg : computes ice velocities … … 11 12 !! * Modules used 12 13 USE phycst 14 USE par_oce 13 15 USE ice_oce ! ice variables 14 16 USE dom_ice 15 17 USE ice 16 18 USE lbclnk 19 USE lib_mpp 17 20 USE in_out_manager ! I/O manager 18 21 … … 24 27 25 28 !! * Module variables 26 27 28 29 REAL(wp) :: & ! constant values 30 rzero = 0.e0 , & 31 rone = 1.e0 29 32 !!---------------------------------------------------------------------- 30 33 !! LIM 2.0 , UCL-LODYC-IPSL (2003) … … 33 36 CONTAINS 34 37 35 SUBROUTINE lim_rhg( k hemi)38 SUBROUTINE lim_rhg( k_j1, k_jpj ) 36 39 !!------------------------------------------------------------------- 37 40 !! *** SUBROUTINR lim_rhg *** 41 !! 38 42 !! ** purpose : determines the velocity field of sea ice by using 39 43 !! atmospheric (wind stress) and oceanic (water stress and surface … … 49 53 !!------------------------------------------------------------------- 50 54 ! * Arguments 51 INTEGER, INTENT(in) :: khemi ! -1/1 = South/North hemisphere flag 55 INTEGER, INTENT(in) :: k_j1 , & ! southern j-index for ice computation 56 & k_jpj ! northern j-index for ice computation 52 57 53 58 ! * Local variables … … 55 60 56 61 INTEGER :: & 57 i_j1, i_j2, i_jpj, i_jpjm1, & ! ????58 62 iim1, ijm1, iip1 , ijp1 , & ! temporary integers 59 63 iter, jter ! " " … … 69 73 zmpzas, zstms, zindu, zindu1, zusdtp, zmassdt, zcorlal, & 70 74 ztrace2, zdeter, zdelta, zsang, zmask, zdgp, zdgi, zdiag 71 REAL(wp),DIMENSION(jpj) :: &72 zind ! i-averaged indicator of sea-ice73 75 REAL(wp),DIMENSION(jpi,jpj) :: & 74 76 zpresh, zfrld, zmass, zcorl, & … … 77 79 zc1, zc2, zd1, zd2, zden, zu_ice, zv_ice, zresr 78 80 REAL(wp),DIMENSION(jpi,jpj,2,2) :: & 79 zs igm11, zsigm12, zsigm22, zsigm2181 zs11, zs12, zs22, zs21 80 82 !!------------------------------------------------------------------- 81 83 82 ! Store initial velocities. 84 ! Store initial velocities 85 ! ------------------------ 83 86 zu0(:,:) = u_ice(:,:) 84 87 zv0(:,:) = v_ice(:,:) 85 88 86 ! Numerical characteristics. 87 ! -------------------------- 88 89 ! Define the j-limits where ice dynamics is computed 90 ! --------------------------------------------------- 91 92 DO jj = 1, jpj 93 zind(jj) = SUM( frld(:,jj) ) ! = FLOAT(jpj) if ocean everwhere on a j-line 94 END DO 95 96 IF( khemi == 1 ) THEN ! Northern hemisphere 97 i_j1 = jeq 98 i_jpj = jpj 99 DO jj = jpj, jeq, -1 100 IF( zind(jj) < FLOAT(jpi) ) i_j1 = jj 101 END DO 102 i_j1 = MAX( 1, i_j1-1) 103 IF( l_ctl .AND. lwp ) WRITE(numout,*) 'lim_rhg : NH i_j1 = ', i_j1, ' ij_pj = ', i_jpj 104 ELSE ! Southern hemisphere 105 i_j1 = 2 106 i_jpj = jpj 107 DO jj = 1, jeq 108 IF( zind(jj) < FLOAT(jpi) ) i_jpj = jj 109 END DO 110 i_jpj = MIN( jpj, i_jpj+2) 111 IF( l_ctl .AND. lwp ) WRITE(numout,*) 'lim_rhg : SH i_j1 = ', i_j1, ' ij_pj = ', i_jpj 112 ENDIF 113 i_j2 = i_j1 + 1 ! inner domain indices 114 i_jpjm1 = i_jpj - 1 115 116 117 ! 2) Sign of turning angle for oceanic drag. | 118 !----------------------------------------------------------------------- 119 120 zsang = REAL( khemi ) * sangvg ! only the sinus changes its sign with the hemisphere 121 122 123 ! 3) Ice mass, ice strength, and wind stress at the center | 124 ! of the grid squares. | 125 !----------------------------------------------------------------------- 126 127 DO jj = i_j1 , i_jpjm1 89 ! Ice mass, ice strength, and wind stress at the center | 90 ! of the grid squares. | 91 !------------------------------------------------------------------- 92 93 DO jj = k_j1 , k_jpj-1 128 94 DO ji = 1 , jpi 129 95 za1(ji,jj) = tms(ji,jj) * ( rhosn * hsnm(ji,jj) + rhoic * hicm(ji,jj) ) … … 145 111 !--------------------------------------------------------------------------- 146 112 147 DO jj = i_j2, i_jpjm1113 DO jj = k_j1+1, k_jpj-1 148 114 DO ji = 2, jpi 149 115 zstms = tms(ji,jj ) * wght(ji,jj,2,2) + tms(ji-1,jj ) * wght(ji,jj,1,2) & … … 211 177 ! Computation of free drift field for free slip boundary conditions. 212 178 213 DO jj = i_j1, i_jpjm1179 DO jj = k_j1, k_jpj-1 214 180 DO ji = 1, jpim1 215 181 !- Rate of strain tensor. … … 242 208 243 209 !- Determination of zc1u, zc2u, zc1v and zc2v. 244 DO jj = i_j2, i_jpjm1210 DO jj = k_j1+1, k_jpj-1 245 211 DO ji = 2, jpim1 246 212 ze11 = akappa(ji-1,jj-1,1,1) … … 254 220 255 221 zdiag = zvis22 * ( ze11 + ze22 ) 256 zs igm11(ji,jj,1,1) = zvis11 * ze11 + zdiag257 zs igm12(ji,jj,1,1) = zvis12 * ze12 + zvis21 * ze21258 zs igm22(ji,jj,1,1) = zvis11 * ze22 + zdiag259 zs igm21(ji,jj,1,1) = zvis12 * ze21 + zvis21 * ze12222 zs11(ji,jj,1,1) = zvis11 * ze11 + zdiag 223 zs12(ji,jj,1,1) = zvis12 * ze12 + zvis21 * ze21 224 zs22(ji,jj,1,1) = zvis11 * ze22 + zdiag 225 zs21(ji,jj,1,1) = zvis12 * ze21 + zvis21 * ze12 260 226 261 227 ze11 = -akappa(ji,jj-1,1,1) … … 269 235 270 236 zdiag = zvis22 * ( ze11 + ze22 ) 271 zs igm11(ji,jj,2,1) = zvis11 * ze11 + zdiag272 zs igm12(ji,jj,2,1) = zvis12 * ze12 + zvis21 * ze21273 zs igm22(ji,jj,2,1) = zvis11 * ze22 + zdiag274 zs igm21(ji,jj,2,1) = zvis12 * ze21 + zvis21 * ze12237 zs11(ji,jj,2,1) = zvis11 * ze11 + zdiag 238 zs12(ji,jj,2,1) = zvis12 * ze12 + zvis21 * ze21 239 zs22(ji,jj,2,1) = zvis11 * ze22 + zdiag 240 zs21(ji,jj,2,1) = zvis12 * ze21 + zvis21 * ze12 275 241 276 242 ze11 = akappa(ji-1,jj,1,1) … … 284 250 285 251 zdiag = zvis22 * ( ze11 + ze22 ) 286 zs igm11(ji,jj,1,2) = zvis11 * ze11 + zdiag287 zs igm12(ji,jj,1,2) = zvis12 * ze12 + zvis21 * ze21288 zs igm22(ji,jj,1,2) = zvis11 * ze22 + zdiag289 zs igm21(ji,jj,1,2) = zvis12 * ze21 + zvis21 * ze12252 zs11(ji,jj,1,2) = zvis11 * ze11 + zdiag 253 zs12(ji,jj,1,2) = zvis12 * ze12 + zvis21 * ze21 254 zs22(ji,jj,1,2) = zvis11 * ze22 + zdiag 255 zs21(ji,jj,1,2) = zvis12 * ze21 + zvis21 * ze12 290 256 291 257 ze11 = -akappa(ji,jj,1,1) … … 299 265 300 266 zdiag = zvis22 * ( ze11 + ze22 ) 301 zs igm11(ji,jj,2,2) = zvis11 * ze11 + zdiag302 zs igm12(ji,jj,2,2) = zvis12 * ze12 + zvis21 * ze21303 zs igm22(ji,jj,2,2) = zvis11 * ze22 + zdiag304 zs igm21(ji,jj,2,2) = zvis12 * ze21 + zvis21 * ze12267 zs11(ji,jj,2,2) = zvis11 * ze11 + zdiag 268 zs12(ji,jj,2,2) = zvis12 * ze12 + zvis21 * ze21 269 zs22(ji,jj,2,2) = zvis11 * ze22 + zdiag 270 zs21(ji,jj,2,2) = zvis12 * ze21 + zvis21 * ze12 305 271 END DO 306 272 END DO 307 273 308 DO jj = i_j2, i_jpjm1274 DO jj = k_j1+1, k_jpj-1 309 275 DO ji = 2, jpim1 310 zc1u(ji,jj) = alambd(ji,jj,2,2,2,1) * zsigm11(ji,jj,2,1) + alambd(ji,jj,2,2,2,2) * zsigm11(ji,jj,2,2) & 311 & - alambd(ji,jj,2,2,1,1) * zsigm11(ji,jj,1,1) - alambd(ji,jj,2,2,1,2) * zsigm11(ji,jj,1,2) & 312 & - alambd(ji,jj,1,1,2,1) * zsigm12(ji,jj,2,1) - alambd(ji,jj,1,1,1,1) * zsigm12(ji,jj,1,1) & 313 & + alambd(ji,jj,1,1,2,2) * zsigm12(ji,jj,2,2) + alambd(ji,jj,1,1,1,2) * zsigm12(ji,jj,1,2) & 314 & + alambd(ji,jj,1,2,1,1) * zsigm21(ji,jj,1,1) + alambd(ji,jj,1,2,2,1) * zsigm21(ji,jj,2,1) & 315 & + alambd(ji,jj,1,2,1,2) * zsigm21(ji,jj,1,2) + alambd(ji,jj,1,2,2,2) * zsigm21(ji,jj,2,2) & 316 & - alambd(ji,jj,2,1,1,1) * zsigm22(ji,jj,1,1) - alambd(ji,jj,2,1,2,1) * zsigm22(ji,jj,2,1) & 317 & - alambd(ji,jj,2,1,1,2) * zsigm22(ji,jj,1,2) - alambd(ji,jj,2,1,2,2) * zsigm22(ji,jj,2,2) 276 zc1u(ji,jj) = & 277 + alambd(ji,jj,2,2,2,1) * zs11(ji,jj,2,1) + alambd(ji,jj,2,2,2,2) * zs11(ji,jj,2,2) & 278 - alambd(ji,jj,2,2,1,1) * zs11(ji,jj,1,1) - alambd(ji,jj,2,2,1,2) * zs11(ji,jj,1,2) & 279 - alambd(ji,jj,1,1,2,1) * zs12(ji,jj,2,1) - alambd(ji,jj,1,1,1,1) * zs12(ji,jj,1,1) & 280 + alambd(ji,jj,1,1,2,2) * zs12(ji,jj,2,2) + alambd(ji,jj,1,1,1,2) * zs12(ji,jj,1,2) & 281 + alambd(ji,jj,1,2,1,1) * zs21(ji,jj,1,1) + alambd(ji,jj,1,2,2,1) * zs21(ji,jj,2,1) & 282 + alambd(ji,jj,1,2,1,2) * zs21(ji,jj,1,2) + alambd(ji,jj,1,2,2,2) * zs21(ji,jj,2,2) & 283 - alambd(ji,jj,2,1,1,1) * zs22(ji,jj,1,1) - alambd(ji,jj,2,1,2,1) * zs22(ji,jj,2,1) & 284 - alambd(ji,jj,2,1,1,2) * zs22(ji,jj,1,2) - alambd(ji,jj,2,1,2,2) * zs22(ji,jj,2,2) 318 285 319 zc2u(ji,jj) = alambd(ji,jj,2,2,2,1) * zsigm21(ji,jj,2,1) + alambd(ji,jj,2,2,2,2) * zsigm21(ji,jj,2,2) & 320 & - alambd(ji,jj,2,2,1,1) * zsigm21(ji,jj,1,1) - alambd(ji,jj,2,2,1,2) * zsigm21(ji,jj,1,2) & 321 & - alambd(ji,jj,1,1,2,1) * zsigm22(ji,jj,2,1) - alambd(ji,jj,1,1,1,1) * zsigm22(ji,jj,1,1) & 322 & + alambd(ji,jj,1,1,2,2) * zsigm22(ji,jj,2,2) + alambd(ji,jj,1,1,1,2) * zsigm22(ji,jj,1,2) & 323 & - alambd(ji,jj,1,2,1,1) * zsigm11(ji,jj,1,1) - alambd(ji,jj,1,2,2,1) * zsigm11(ji,jj,2,1) & 324 & - alambd(ji,jj,1,2,1,2) * zsigm11(ji,jj,1,2) - alambd(ji,jj,1,2,2,2) * zsigm11(ji,jj,2,2) & 325 & + alambd(ji,jj,2,1,1,1) * zsigm12(ji,jj,1,1) + alambd(ji,jj,2,1,2,1) * zsigm12(ji,jj,2,1) & 326 & + alambd(ji,jj,2,1,1,2) * zsigm12(ji,jj,1,2) + alambd(ji,jj,2,1,2,2) * zsigm12(ji,jj,2,2) 286 zc2u(ji,jj) = & 287 + alambd(ji,jj,2,2,2,1) * zs21(ji,jj,2,1) + alambd(ji,jj,2,2,2,2) * zs21(ji,jj,2,2) & 288 - alambd(ji,jj,2,2,1,1) * zs21(ji,jj,1,1) - alambd(ji,jj,2,2,1,2) * zs21(ji,jj,1,2) & 289 - alambd(ji,jj,1,1,2,1) * zs22(ji,jj,2,1) - alambd(ji,jj,1,1,1,1) * zs22(ji,jj,1,1) & 290 + alambd(ji,jj,1,1,2,2) * zs22(ji,jj,2,2) + alambd(ji,jj,1,1,1,2) * zs22(ji,jj,1,2) & 291 - alambd(ji,jj,1,2,1,1) * zs11(ji,jj,1,1) - alambd(ji,jj,1,2,2,1) * zs11(ji,jj,2,1) & 292 - alambd(ji,jj,1,2,1,2) * zs11(ji,jj,1,2) - alambd(ji,jj,1,2,2,2) * zs11(ji,jj,2,2) & 293 + alambd(ji,jj,2,1,1,1) * zs12(ji,jj,1,1) + alambd(ji,jj,2,1,2,1) * zs12(ji,jj,2,1) & 294 + alambd(ji,jj,2,1,1,2) * zs12(ji,jj,1,2) + alambd(ji,jj,2,1,2,2) * zs12(ji,jj,2,2) 327 295 END DO 328 296 END DO 329 297 330 DO jj = i_j2, i_jpjm1298 DO jj = k_j1+1, k_jpj-1 331 299 DO ji = 2, jpim1 332 300 ! zc1v , zc2v. … … 341 309 342 310 zdiag = zvis22 * ( ze11 + ze22 ) 343 zs igm11(ji,jj,1,1) = zvis11 * ze11 + zdiag344 zs igm12(ji,jj,1,1) = zvis12 * ze12 + zvis21 * ze21345 zs igm22(ji,jj,1,1) = zvis11 * ze22 + zdiag346 zs igm21(ji,jj,1,1) = zvis12 * ze21 + zvis21 * ze12311 zs11(ji,jj,1,1) = zvis11 * ze11 + zdiag 312 zs12(ji,jj,1,1) = zvis12 * ze12 + zvis21 * ze21 313 zs22(ji,jj,1,1) = zvis11 * ze22 + zdiag 314 zs21(ji,jj,1,1) = zvis12 * ze21 + zvis21 * ze12 347 315 348 316 ze11 = akappa(ji,jj-1,1,2) … … 356 324 357 325 zdiag = zvis22 * ( ze11 + ze22 ) 358 zs igm11(ji,jj,2,1) = zvis11 * ze11 + zdiag359 zs igm12(ji,jj,2,1) = zvis12 * ze12 + zvis21 * ze21360 zs igm22(ji,jj,2,1) = zvis11 * ze22 + zdiag361 zs igm21(ji,jj,2,1) = zvis12 * ze21 + zvis21 * ze12326 zs11(ji,jj,2,1) = zvis11 * ze11 + zdiag 327 zs12(ji,jj,2,1) = zvis12 * ze12 + zvis21 * ze21 328 zs22(ji,jj,2,1) = zvis11 * ze22 + zdiag 329 zs21(ji,jj,2,1) = zvis12 * ze21 + zvis21 * ze12 362 330 363 331 ze11 = akappa(ji-1,jj,1,2) … … 371 339 372 340 zdiag = zvis22 * ( ze11 + ze22 ) 373 zs igm11(ji,jj,1,2) = zvis11 * ze11 + zdiag374 zs igm12(ji,jj,1,2) = zvis12 * ze12 + zvis21 * ze21375 zs igm22(ji,jj,1,2) = zvis11 * ze22 + zdiag376 zs igm21(ji,jj,1,2) = zvis12 * ze21 + zvis21 * ze12341 zs11(ji,jj,1,2) = zvis11 * ze11 + zdiag 342 zs12(ji,jj,1,2) = zvis12 * ze12 + zvis21 * ze21 343 zs22(ji,jj,1,2) = zvis11 * ze22 + zdiag 344 zs21(ji,jj,1,2) = zvis12 * ze21 + zvis21 * ze12 377 345 378 346 ze11 = akappa(ji,jj,1,2) … … 386 354 387 355 zdiag = zvis22 * ( ze11 + ze22 ) 388 zs igm11(ji,jj,2,2) = zvis11 * ze11 + zdiag389 zs igm12(ji,jj,2,2) = zvis12 * ze12 + zvis21 * ze21390 zs igm22(ji,jj,2,2) = zvis11 * ze22 + zdiag391 zs igm21(ji,jj,2,2) = zvis12 * ze21 + zvis21 * ze12356 zs11(ji,jj,2,2) = zvis11 * ze11 + zdiag 357 zs12(ji,jj,2,2) = zvis12 * ze12 + zvis21 * ze21 358 zs22(ji,jj,2,2) = zvis11 * ze22 + zdiag 359 zs21(ji,jj,2,2) = zvis12 * ze21 + zvis21 * ze12 392 360 393 361 END DO 394 362 END DO 395 363 396 DO jj = i_j2, i_jpjm1364 DO jj = k_j1+1, k_jpj-1 397 365 DO ji = 2, jpim1 398 zc1v(ji,jj) = alambd(ji,jj,2,2,2,1) * zsigm11(ji,jj,2,1) + alambd(ji,jj,2,2,2,2) * zsigm11(ji,jj,2,2) & 399 & - alambd(ji,jj,2,2,1,1) * zsigm11(ji,jj,1,1) - alambd(ji,jj,2,2,1,2) * zsigm11(ji,jj,1,2) & 400 & - alambd(ji,jj,1,1,2,1) * zsigm12(ji,jj,2,1) - alambd(ji,jj,1,1,1,1) * zsigm12(ji,jj,1,1) & 401 & + alambd(ji,jj,1,1,2,2) * zsigm12(ji,jj,2,2) + alambd(ji,jj,1,1,1,2) * zsigm12(ji,jj,1,2) & 402 & + alambd(ji,jj,1,2,1,1) * zsigm21(ji,jj,1,1) + alambd(ji,jj,1,2,2,1) * zsigm21(ji,jj,2,1) & 403 & + alambd(ji,jj,1,2,1,2) * zsigm21(ji,jj,1,2) + alambd(ji,jj,1,2,2,2) * zsigm21(ji,jj,2,2) & 404 & - alambd(ji,jj,2,1,1,1) * zsigm22(ji,jj,1,1) - alambd(ji,jj,2,1,2,1) * zsigm22(ji,jj,2,1) & 405 & - alambd(ji,jj,2,1,1,2) * zsigm22(ji,jj,1,2) - alambd(ji,jj,2,1,2,2) * zsigm22(ji,jj,2,2) 406 zc2v(ji,jj) = alambd(ji,jj,2,2,2,1) * zsigm21(ji,jj,2,1) + alambd(ji,jj,2,2,2,2) * zsigm21(ji,jj,2,2) & 407 & - alambd(ji,jj,2,2,1,1) * zsigm21(ji,jj,1,1) - alambd(ji,jj,2,2,1,2) * zsigm21(ji,jj,1,2) & 408 & - alambd(ji,jj,1,1,2,1) * zsigm22(ji,jj,2,1) - alambd(ji,jj,1,1,1,1) * zsigm22(ji,jj,1,1) & 409 & + alambd(ji,jj,1,1,2,2) * zsigm22(ji,jj,2,2) + alambd(ji,jj,1,1,1,2) * zsigm22(ji,jj,1,2) & 410 & - alambd(ji,jj,1,2,1,1) * zsigm11(ji,jj,1,1) - alambd(ji,jj,1,2,2,1) * zsigm11(ji,jj,2,1) & 411 & - alambd(ji,jj,1,2,1,2) * zsigm11(ji,jj,1,2) - alambd(ji,jj,1,2,2,2) * zsigm11(ji,jj,2,2) & 412 & + alambd(ji,jj,2,1,1,1) * zsigm12(ji,jj,1,1) + alambd(ji,jj,2,1,2,1) * zsigm12(ji,jj,2,1) & 413 & + alambd(ji,jj,2,1,1,2) * zsigm12(ji,jj,1,2) + alambd(ji,jj,2,1,2,2) * zsigm12(ji,jj,2,2) 366 zc1v(ji,jj) = & 367 + alambd(ji,jj,2,2,2,1) * zs11(ji,jj,2,1) + alambd(ji,jj,2,2,2,2) * zs11(ji,jj,2,2) & 368 - alambd(ji,jj,2,2,1,1) * zs11(ji,jj,1,1) - alambd(ji,jj,2,2,1,2) * zs11(ji,jj,1,2) & 369 - alambd(ji,jj,1,1,2,1) * zs12(ji,jj,2,1) - alambd(ji,jj,1,1,1,1) * zs12(ji,jj,1,1) & 370 + alambd(ji,jj,1,1,2,2) * zs12(ji,jj,2,2) + alambd(ji,jj,1,1,1,2) * zs12(ji,jj,1,2) & 371 + alambd(ji,jj,1,2,1,1) * zs21(ji,jj,1,1) + alambd(ji,jj,1,2,2,1) * zs21(ji,jj,2,1) & 372 + alambd(ji,jj,1,2,1,2) * zs21(ji,jj,1,2) + alambd(ji,jj,1,2,2,2) * zs21(ji,jj,2,2) & 373 - alambd(ji,jj,2,1,1,1) * zs22(ji,jj,1,1) - alambd(ji,jj,2,1,2,1) * zs22(ji,jj,2,1) & 374 - alambd(ji,jj,2,1,1,2) * zs22(ji,jj,1,2) - alambd(ji,jj,2,1,2,2) * zs22(ji,jj,2,2) 375 zc2v(ji,jj) = & 376 + alambd(ji,jj,2,2,2,1) * zs21(ji,jj,2,1) + alambd(ji,jj,2,2,2,2) * zs21(ji,jj,2,2) & 377 - alambd(ji,jj,2,2,1,1) * zs21(ji,jj,1,1) - alambd(ji,jj,2,2,1,2) * zs21(ji,jj,1,2) & 378 - alambd(ji,jj,1,1,2,1) * zs22(ji,jj,2,1) - alambd(ji,jj,1,1,1,1) * zs22(ji,jj,1,1) & 379 + alambd(ji,jj,1,1,2,2) * zs22(ji,jj,2,2) + alambd(ji,jj,1,1,1,2) * zs22(ji,jj,1,2) & 380 - alambd(ji,jj,1,2,1,1) * zs11(ji,jj,1,1) - alambd(ji,jj,1,2,2,1) * zs11(ji,jj,2,1) & 381 - alambd(ji,jj,1,2,1,2) * zs11(ji,jj,1,2) - alambd(ji,jj,1,2,2,2) * zs11(ji,jj,2,2) & 382 + alambd(ji,jj,2,1,1,1) * zs12(ji,jj,1,1) + alambd(ji,jj,2,1,2,1) * zs12(ji,jj,2,1) & 383 + alambd(ji,jj,2,1,1,2) * zs12(ji,jj,1,2) + alambd(ji,jj,2,1,2,2) * zs12(ji,jj,2,2) 414 384 END DO 415 385 END DO … … 420 390 421 391 ! Store previous drift field. 422 DO jj = i_j1, i_jpjm1392 DO jj = k_j1, k_jpj-1 423 393 zu_ice(:,jj) = u_ice(:,jj) 424 394 zv_ice(:,jj) = v_ice(:,jj) 425 395 END DO 426 396 427 DO jj = i_j2, i_jpjm1 428 DO ji = 2, jpim1 429 zur = u_ice(ji,jj) - u_oce(ji,jj) 430 zvr = v_ice(ji,jj) - v_oce(ji,jj) 431 zmod = SQRT( zur * zur + zvr * zvr) * ( 1.0 - zfrld(ji,jj) ) 432 za = rhoco * zmod 433 zac = za * cangvg 397 DO jj = k_j1+1, k_jpj-1 398 zsang = SIGN( 1.e0, fcor(1,jj) ) * sangvg ! only the sinus changes its sign with the hemisphere 399 DO ji = 2, jpim1 400 zur = u_ice(ji,jj) - u_oce(ji,jj) 401 zvr = v_ice(ji,jj) - v_oce(ji,jj) 402 zmod = SQRT( zur * zur + zvr * zvr) * ( 1.0 - zfrld(ji,jj) ) 403 za = rhoco * zmod 404 zac = za * cangvg 434 405 zmpzas = alpha * zcorl(ji,jj) + za * zsang 435 406 zmassdt = zusdtp * zmass(ji,jj) … … 456 427 ! Terms that do not involve already up-dated velocities. 457 428 458 DO jj = i_j2, i_jpjm1429 DO jj = k_j1+1, k_jpj-1 459 430 DO ji = 2, jpim1 460 431 iim1 = ji … … 471 442 zvis21 = zviseta (iim1,ijm1) 472 443 zdiag = zvis22 * ( ze11 + ze22 ) 473 zs igm11(ji,jj,2,1) = zvis11 * ze11 + zdiag474 zs igm12(ji,jj,2,1) = zvis12 * ze12 + zvis21 * ze21475 zs igm22(ji,jj,2,1) = zvis11 * ze22 + zdiag476 zs igm21(ji,jj,2,1) = zvis12 * ze21 + zvis21 * ze12444 zs11(ji,jj,2,1) = zvis11 * ze11 + zdiag 445 zs12(ji,jj,2,1) = zvis12 * ze12 + zvis21 * ze21 446 zs22(ji,jj,2,1) = zvis11 * ze22 + zdiag 447 zs21(ji,jj,2,1) = zvis12 * ze21 + zvis21 * ze12 477 448 478 449 … … 494 465 zvis21 = zviseta (iim1,ijm1) 495 466 zdiag = zvis22 * ( ze11 + ze22 ) 496 zs igm11(ji,jj,1,2) = zvis11 * ze11 + zdiag497 zs igm12(ji,jj,1,2) = zvis12 * ze12 + zvis21 * ze21498 zs igm22(ji,jj,1,2) = zvis11 * ze22 + zdiag499 zs igm21(ji,jj,1,2) = zvis12 * ze21 + zvis21 * ze12467 zs11(ji,jj,1,2) = zvis11 * ze11 + zdiag 468 zs12(ji,jj,1,2) = zvis12 * ze12 + zvis21 * ze21 469 zs22(ji,jj,1,2) = zvis11 * ze22 + zdiag 470 zs21(ji,jj,1,2) = zvis12 * ze21 + zvis21 * ze12 500 471 501 472 iim1 = ji … … 517 488 518 489 zdiag = zvis22 * ( ze11 + ze22 ) 519 zs igm11(ji,jj,2,2) = zvis11 * ze11 + zdiag520 zs igm12(ji,jj,2,2) = zvis12 * ze12 + zvis21 * ze21521 zs igm22(ji,jj,2,2) = zvis11 * ze22 + zdiag522 zs igm21(ji,jj,2,2) = zvis12 * ze21 + zvis21 * ze12490 zs11(ji,jj,2,2) = zvis11 * ze11 + zdiag 491 zs12(ji,jj,2,2) = zvis12 * ze12 + zvis21 * ze21 492 zs22(ji,jj,2,2) = zvis11 * ze22 + zdiag 493 zs21(ji,jj,2,2) = zvis12 * ze21 + zvis21 * ze12 523 494 524 495 END DO … … 529 500 ! Using arrays u and v in the computation of the terms ze leads to GAUSS-SEIDEL method. 530 501 531 DO jj = i_j2, i_jpjm1502 DO jj = k_j1+1, k_jpj-1 532 503 DO ji = 2, jpim1 533 504 iim1 = ji - 1 … … 549 520 550 521 zdiag = zvis22 * ( ze11 + ze22 ) 551 zs igm11(ji,jj,1,1) = zvis11 * ze11 + zdiag552 zs igm12(ji,jj,1,1) = zvis12 * ze12 + zvis21 * ze21553 zs igm22(ji,jj,1,1) = zvis11 * ze22 + zdiag554 zs igm21(ji,jj,1,1) = zvis12 * ze21 + zvis21 * ze12522 zs11(ji,jj,1,1) = zvis11 * ze11 + zdiag 523 zs12(ji,jj,1,1) = zvis12 * ze12 + zvis21 * ze21 524 zs22(ji,jj,1,1) = zvis11 * ze22 + zdiag 525 zs21(ji,jj,1,1) = zvis12 * ze21 + zvis21 * ze12 555 526 556 527 … … 572 543 573 544 zdiag = zvis22 * ( ze11 + ze22 ) 574 zs igm11(ji,jj,2,1) = zsigm11(ji,jj,2,1) + zvis11 * ze11 + zdiag575 zs igm12(ji,jj,2,1) = zsigm12(ji,jj,2,1) + zvis12 * ze12 + zvis21 * ze21576 zs igm22(ji,jj,2,1) = zsigm22(ji,jj,2,1) + zvis11 * ze22 + zdiag577 zs igm21(ji,jj,2,1) = zsigm21(ji,jj,2,1) + zvis12 * ze21 + zvis21 * ze12545 zs11(ji,jj,2,1) = zs11(ji,jj,2,1) + zvis11 * ze11 + zdiag 546 zs12(ji,jj,2,1) = zs12(ji,jj,2,1) + zvis12 * ze12 + zvis21 * ze21 547 zs22(ji,jj,2,1) = zs22(ji,jj,2,1) + zvis11 * ze22 + zdiag 548 zs21(ji,jj,2,1) = zs21(ji,jj,2,1) + zvis12 * ze21 + zvis21 * ze12 578 549 579 550 … … 590 561 591 562 zdiag = zvis22 * ( ze11 + ze22 ) 592 zs igm11(ji,jj,1,2) = zsigm11(ji,jj,1,2) + zvis11 * ze11 + zdiag593 zs igm12(ji,jj,1,2) = zsigm12(ji,jj,1,2) + zvis12 * ze12 + zvis21 * ze21594 zs igm22(ji,jj,1,2) = zsigm22(ji,jj,1,2) + zvis11 * ze22 + zdiag595 zs igm21(ji,jj,1,2) = zsigm21(ji,jj,1,2) + zvis12 * ze21 + zvis21 * ze12563 zs11(ji,jj,1,2) = zs11(ji,jj,1,2) + zvis11 * ze11 + zdiag 564 zs12(ji,jj,1,2) = zs12(ji,jj,1,2) + zvis12 * ze12 + zvis21 * ze21 565 zs22(ji,jj,1,2) = zs22(ji,jj,1,2) + zvis11 * ze22 + zdiag 566 zs21(ji,jj,1,2) = zs21(ji,jj,1,2) + zvis12 * ze21 + zvis21 * ze12 596 567 597 568 !i END DO 598 569 !i END DO 599 570 600 !i DO jj = i_j2, i_jpjm1571 !i DO jj = k_j1+1, k_jpj-1 601 572 !i DO ji = 2, jpim1 602 zd1(ji,jj) = alambd(ji,jj,2,2,2,1) * zsigm11(ji,jj,2,1) + alambd(ji,jj,2,2,2,2) * zsigm11(ji,jj,2,2) & 603 & - alambd(ji,jj,2,2,1,1) * zsigm11(ji,jj,1,1) - alambd(ji,jj,2,2,1,2) * zsigm11(ji,jj,1,2) & 604 & - alambd(ji,jj,1,1,2,1) * zsigm12(ji,jj,2,1) - alambd(ji,jj,1,1,1,1) * zsigm12(ji,jj,1,1) & 605 & + alambd(ji,jj,1,1,2,2) * zsigm12(ji,jj,2,2) + alambd(ji,jj,1,1,1,2) * zsigm12(ji,jj,1,2) & 606 & + alambd(ji,jj,1,2,1,1) * zsigm21(ji,jj,1,1) + alambd(ji,jj,1,2,2,1) * zsigm21(ji,jj,2,1) & 607 & + alambd(ji,jj,1,2,1,2) * zsigm21(ji,jj,1,2) + alambd(ji,jj,1,2,2,2) * zsigm21(ji,jj,2,2) & 608 & - alambd(ji,jj,2,1,1,1) * zsigm22(ji,jj,1,1) - alambd(ji,jj,2,1,2,1) * zsigm22(ji,jj,2,1) & 609 & - alambd(ji,jj,2,1,1,2) * zsigm22(ji,jj,1,2) - alambd(ji,jj,2,1,2,2) * zsigm22(ji,jj,2,2) 610 611 zd2(ji,jj) = alambd(ji,jj,2,2,2,1) * zsigm21(ji,jj,2,1) + alambd(ji,jj,2,2,2,2) * zsigm21(ji,jj,2,2) & 612 & - alambd(ji,jj,2,2,1,1) * zsigm21(ji,jj,1,1) - alambd(ji,jj,2,2,1,2) * zsigm21(ji,jj,1,2) & 613 & - alambd(ji,jj,1,1,2,1) * zsigm22(ji,jj,2,1) - alambd(ji,jj,1,1,1,1) * zsigm22(ji,jj,1,1) & 614 & + alambd(ji,jj,1,1,2,2) * zsigm22(ji,jj,2,2) + alambd(ji,jj,1,1,1,2) * zsigm22(ji,jj,1,2) & 615 & - alambd(ji,jj,1,2,1,1) * zsigm11(ji,jj,1,1) - alambd(ji,jj,1,2,2,1) * zsigm11(ji,jj,2,1) & 616 & - alambd(ji,jj,1,2,1,2) * zsigm11(ji,jj,1,2) - alambd(ji,jj,1,2,2,2) * zsigm11(ji,jj,2,2) & 617 & + alambd(ji,jj,2,1,1,1) * zsigm12(ji,jj,1,1) + alambd(ji,jj,2,1,2,1) * zsigm12(ji,jj,2,1) & 618 & + alambd(ji,jj,2,1,1,2) * zsigm12(ji,jj,1,2) + alambd(ji,jj,2,1,2,2) * zsigm12(ji,jj,2,2) 573 zd1(ji,jj) = & 574 + alambd(ji,jj,2,2,2,1) * zs11(ji,jj,2,1) + alambd(ji,jj,2,2,2,2) * zs11(ji,jj,2,2) & 575 - alambd(ji,jj,2,2,1,1) * zs11(ji,jj,1,1) - alambd(ji,jj,2,2,1,2) * zs11(ji,jj,1,2) & 576 - alambd(ji,jj,1,1,2,1) * zs12(ji,jj,2,1) - alambd(ji,jj,1,1,1,1) * zs12(ji,jj,1,1) & 577 + alambd(ji,jj,1,1,2,2) * zs12(ji,jj,2,2) + alambd(ji,jj,1,1,1,2) * zs12(ji,jj,1,2) & 578 + alambd(ji,jj,1,2,1,1) * zs21(ji,jj,1,1) + alambd(ji,jj,1,2,2,1) * zs21(ji,jj,2,1) & 579 + alambd(ji,jj,1,2,1,2) * zs21(ji,jj,1,2) + alambd(ji,jj,1,2,2,2) * zs21(ji,jj,2,2) & 580 - alambd(ji,jj,2,1,1,1) * zs22(ji,jj,1,1) - alambd(ji,jj,2,1,2,1) * zs22(ji,jj,2,1) & 581 - alambd(ji,jj,2,1,1,2) * zs22(ji,jj,1,2) - alambd(ji,jj,2,1,2,2) * zs22(ji,jj,2,2) 582 zd2(ji,jj) = & 583 + alambd(ji,jj,2,2,2,1) * zs21(ji,jj,2,1) + alambd(ji,jj,2,2,2,2) * zs21(ji,jj,2,2) & 584 - alambd(ji,jj,2,2,1,1) * zs21(ji,jj,1,1) - alambd(ji,jj,2,2,1,2) * zs21(ji,jj,1,2) & 585 - alambd(ji,jj,1,1,2,1) * zs22(ji,jj,2,1) - alambd(ji,jj,1,1,1,1) * zs22(ji,jj,1,1) & 586 + alambd(ji,jj,1,1,2,2) * zs22(ji,jj,2,2) + alambd(ji,jj,1,1,1,2) * zs22(ji,jj,1,2) & 587 - alambd(ji,jj,1,2,1,1) * zs11(ji,jj,1,1) - alambd(ji,jj,1,2,2,1) * zs11(ji,jj,2,1) & 588 - alambd(ji,jj,1,2,1,2) * zs11(ji,jj,1,2) - alambd(ji,jj,1,2,2,2) * zs11(ji,jj,2,2) & 589 + alambd(ji,jj,2,1,1,1) * zs12(ji,jj,1,1) + alambd(ji,jj,2,1,2,1) * zs12(ji,jj,2,1) & 590 + alambd(ji,jj,2,1,1,2) * zs12(ji,jj,1,2) + alambd(ji,jj,2,1,2,2) * zs12(ji,jj,2,2) 619 591 END DO 620 592 END DO 621 593 622 DO jj = i_j2, i_jpjm1594 DO jj = k_j1+1, k_jpj-1 623 595 DO ji = 2, jpim1 624 596 zunw = ( ( za1(ji,jj) + zd1(ji,jj) ) * zc2(ji,jj) & … … 639 611 640 612 !--- 5.2.5.4. Convergence test. 641 DO jj = i_j2 , i_jpjm1613 DO jj = k_j1+1 , k_jpj-1 642 614 zresr(:,jj) = MAX( ABS( u_ice(:,jj) - zu_ice(:,jj) ) , ABS( v_ice(:,jj) - zv_ice(:,jj) ) ) 643 615 END DO 644 zresm = MAXVAL( zresr( 1:jpi , i_j2:i_jpjm1 ) ) 616 zresm = MAXVAL( zresr( 1:jpi , k_j1+1:k_jpj-1 ) ) 617 IF( lk_mpp ) CALL mpp_max( zresm ) ! max over the global domain 645 618 646 619 IF ( zresm <= resl) EXIT iflag … … 649 622 650 623 zindu1 = 1.0 - zindu 651 DO jj = i_j1 , i_jpjm1624 DO jj = k_j1 , k_jpj-1 652 625 zu0(:,jj) = zindu * zu0(:,jj) + zindu1 * u_ice(:,jj) 653 626 zv0(:,jj) = zindu * zv0(:,jj) + zindu1 * v_ice(:,jj) … … 657 630 ! ! ==================== ! 658 631 659 IF( l_ctl .AND. lwp) THEN632 IF(l_ctl) THEN 660 633 WRITE(numout,*) ' lim_rhg : res= ', zresm, 'iter= ', jter,' u_ice ', SUM( u_ice ) , ' v_ice ', SUM( v_ice ) 661 634 ENDIF 662 635 663 636 END SUBROUTINE lim_rhg 637 664 638 #else 639 !!---------------------------------------------------------------------- 640 !! Default option Dummy module NO LIM sea-ice model 641 !!---------------------------------------------------------------------- 642 CONTAINS 643 SUBROUTINE lim_rhg( k1 , k2 ) ! Dummy routine 644 WRITE(*,*) 'lim_rhg: You should not have seen this print! error?', k1, k2 645 END SUBROUTINE lim_rhg 646 #endif 647 665 648 !!============================================================================== 666 !! *** MODULE limrhg ***667 !! No sea ice668 !!==============================================================================669 CONTAINS670 SUBROUTINE lim_rhg ! Empty routine671 END SUBROUTINE lim_rhg672 673 #endif674 675 649 END MODULE limrhg
Note: See TracChangeset
for help on using the changeset viewer.
