00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042 module ice_dyn_evp
00043
00044
00045
00046 use ice_kinds_mod
00047 use ice_fileunits
00048 use ice_communicate, only: my_task, master_task, MPI_COMM_ICE
00049 use ice_domain_size
00050 use ice_constants
00051 use perf_mod, only: t_startf, t_stopf, t_barrierf
00052
00053
00054
00055 implicit none
00056 save
00057
00058
00059
00060 integer (kind=int_kind) ::
00061 kdyn ,
00062 ndte
00063
00064 logical (kind=log_kind) ::
00065 evp_damping
00066
00067
00068
00069 character (len=char_len) ::
00070 yield_curve
00071
00072 real (kind=dbl_kind), parameter ::
00073 dragw = dragio * rhow,
00074
00075 eyc = 0.36_dbl_kind, &
00076
00077 cosw = c1 , &
00078 sinw = c0 , &
00079 a_min = p001, &
00080 m_min = p01
00081
00082 real (kind=dbl_kind) ::
00083 ecci ,
00084 dtei ,
00085 dte2T ,
00086 denom1 ,
00087 denom2 ,
00088 rcon
00089
00090 real (kind=dbl_kind), allocatable ::
00091 fcor_blk(:,:,:)
00092
00093
00094
00095 contains
00096
00097
00098
00099
00100
00101
00102
00103
00104 subroutine evp (dt)
00105
00106
00107
00108
00109
00110 #ifdef CICE_IN_NEMO
00111
00112
00113
00114
00115 #endif
00116
00117
00118
00119
00120
00121
00122
00123 use ice_boundary
00124 use ice_blocks
00125 use ice_domain
00126 use ice_state
00127 use ice_flux
00128 use ice_grid
00129 use ice_timers
00130 use ice_mechred, only: ice_strength
00131
00132
00133
00134 real (kind=dbl_kind), intent(in) ::
00135 dt
00136
00137
00138
00139 integer (kind=int_kind) ::
00140 ksub ,
00141 iblk ,
00142 ilo,ihi,jlo,jhi,
00143 i, j
00144
00145 integer (kind=int_kind), dimension(max_blocks) ::
00146 icellt ,
00147 icellu
00148
00149 integer (kind=int_kind), dimension (nx_block*ny_block, max_blocks) ::
00150 indxti ,
00151 indxtj ,
00152 indxui ,
00153 indxuj
00154
00155 real (kind=dbl_kind), dimension (nx_block,ny_block,max_blocks) ::
00156 tmass ,
00157 waterx ,
00158 watery ,
00159 forcex ,
00160 forcey ,
00161 aiu ,
00162 umass ,
00163 umassdtei
00164
00165 real (kind=dbl_kind), allocatable :: fld2(:,:,:,:)
00166
00167 real (kind=dbl_kind), dimension(nx_block,ny_block,8)::
00168 str
00169
00170 integer (kind=int_kind), dimension (nx_block,ny_block,max_blocks) ::
00171 icetmask
00172
00173 type (block) ::
00174 this_block
00175
00176 call ice_timer_start(timer_dynamics)
00177
00178
00179
00180
00181
00182
00183
00184
00185
00186
00187
00188
00189
00190
00191
00192
00193
00194
00195
00196
00197
00198
00199
00200
00201
00202
00203
00204
00205 do iblk = 1, nblocks
00206
00207 do j = 1, ny_block
00208 do i = 1, nx_block
00209 rdg_conv (i,j,iblk) = c0
00210 rdg_shear(i,j,iblk) = c0
00211 divu (i,j,iblk) = c0
00212 shear(i,j,iblk) = c0
00213 prs_sig(i,j,iblk) = c0
00214 enddo
00215 enddo
00216
00217
00218
00219
00220
00221 this_block = get_block(blocks_ice(iblk),iblk)
00222 ilo = this_block%ilo
00223 ihi = this_block%ihi
00224 jlo = this_block%jlo
00225 jhi = this_block%jhi
00226
00227 call evp_prep1 (nx_block, ny_block, &
00228 ilo, ihi, jlo, jhi, &
00229 aice (:,:,iblk), vice (:,:,iblk), &
00230 vsno (:,:,iblk), tmask (:,:,iblk), &
00231 strairxT_accum(:,:,iblk), strairyT_accum(:,:,iblk), &
00232 strairx (:,:,iblk), strairy (:,:,iblk), &
00233 tmass (:,:,iblk), icetmask(:,:,iblk))
00234
00235 enddo
00236
00237
00238
00239
00240
00241
00242 call ice_timer_start(timer_bound)
00243 call ice_HaloUpdate (icetmask, halo_info, &
00244 field_loc_center, field_type_scalar)
00245 call ice_timer_stop(timer_bound)
00246
00247
00248
00249
00250
00251
00252
00253
00254
00255 call to_ugrid(tmass,umass)
00256 call to_ugrid(aice, aiu)
00257
00258 #ifdef CICE_IN_NEMO
00259
00260
00261
00262
00263
00264 strairx(:,:,:) = strax(:,:,:)
00265 strairy(:,:,:) = stray(:,:,:)
00266 #else
00267 call t2ugrid_vector(strairx)
00268 call t2ugrid_vector(strairy)
00269 #endif
00270
00271
00272
00273
00274
00275
00276 do iblk = 1, nblocks
00277
00278
00279
00280
00281
00282 this_block = get_block(blocks_ice(iblk),iblk)
00283 ilo = this_block%ilo
00284 ihi = this_block%ihi
00285 jlo = this_block%jlo
00286 jhi = this_block%jhi
00287
00288 call evp_prep2 (nx_block, ny_block, &
00289 ilo, ihi, jlo, jhi, &
00290 icellt(iblk), icellu(iblk), &
00291 indxti (:,iblk), indxtj (:,iblk), &
00292 indxui (:,iblk), indxuj (:,iblk), &
00293 aiu (:,:,iblk), umass (:,:,iblk), &
00294 umassdtei (:,:,iblk), fcor_blk (:,:,iblk), &
00295 umask (:,:,iblk), &
00296 uocn (:,:,iblk), vocn (:,:,iblk), &
00297 strairx (:,:,iblk), strairy (:,:,iblk), &
00298 ss_tltx (:,:,iblk), ss_tlty (:,:,iblk), &
00299 icetmask (:,:,iblk), iceumask (:,:,iblk), &
00300 fm (:,:,iblk), &
00301 strtltx (:,:,iblk), strtlty (:,:,iblk), &
00302 strocnx (:,:,iblk), strocny (:,:,iblk), &
00303 strintx (:,:,iblk), strinty (:,:,iblk), &
00304 waterx (:,:,iblk), watery (:,:,iblk), &
00305 forcex (:,:,iblk), forcey (:,:,iblk), &
00306 stressp_1 (:,:,iblk), stressp_2 (:,:,iblk), &
00307 stressp_3 (:,:,iblk), stressp_4 (:,:,iblk), &
00308 stressm_1 (:,:,iblk), stressm_2 (:,:,iblk), &
00309 stressm_3 (:,:,iblk), stressm_4 (:,:,iblk), &
00310 stress12_1(:,:,iblk), stress12_2(:,:,iblk), &
00311 stress12_3(:,:,iblk), stress12_4(:,:,iblk), &
00312 uvel (:,:,iblk), vvel (:,:,iblk))
00313
00314
00315
00316
00317
00318 call ice_strength (nx_block, ny_block, &
00319 ilo, ihi, jlo, jhi, &
00320 icellt(iblk), &
00321 indxti (:,iblk), &
00322 indxtj (:,iblk), &
00323 aice (:,:, iblk), &
00324 vice (:,:, iblk), &
00325 aice0 (:,:, iblk), &
00326 aicen (:,:,:,iblk), &
00327 vicen (:,:,:,iblk), &
00328 strength(:,:, iblk) )
00329
00330 enddo
00331
00332
00333 allocate(fld2(nx_block,ny_block,2,max_blocks))
00334
00335
00336 do iblk = 1,nblocks
00337 fld2(1:nx_block,1:ny_block,1,iblk) = uvel(1:nx_block,1:ny_block,iblk)
00338 fld2(1:nx_block,1:ny_block,2,iblk) = vvel(1:nx_block,1:ny_block,iblk)
00339 enddo
00340
00341
00342
00343
00344
00345
00346 call ice_timer_start(timer_bound)
00347 call ice_HaloUpdate (strength, halo_info, &
00348 field_loc_center, field_type_scalar)
00349
00350 call ice_HaloUpdate (fld2, halo_info, &
00351 field_loc_NEcorner, field_type_vector)
00352
00353
00354
00355
00356 call ice_timer_stop(timer_bound)
00357
00358
00359 do iblk = 1,nblocks
00360 uvel(1:nx_block,1:ny_block,iblk) = fld2(1:nx_block,1:ny_block,1,iblk)
00361 vvel(1:nx_block,1:ny_block,iblk) = fld2(1:nx_block,1:ny_block,2,iblk)
00362 enddo
00363
00364
00365
00366
00367
00368 do ksub = 1,ndte
00369
00370
00371
00372
00373
00374
00375
00376
00377 do iblk = 1, nblocks
00378
00379
00380 call stress (nx_block, ny_block, &
00381 ksub, icellt(iblk), &
00382 indxti (:,iblk), indxtj (:,iblk), &
00383 uvel (:,:,iblk), vvel (:,:,iblk), &
00384 dxt (:,:,iblk), dyt (:,:,iblk), &
00385 dxhy (:,:,iblk), dyhx (:,:,iblk), &
00386 cxp (:,:,iblk), cyp (:,:,iblk), &
00387 cxm (:,:,iblk), cym (:,:,iblk), &
00388 tarear (:,:,iblk), tinyarea (:,:,iblk), &
00389 strength (:,:,iblk), &
00390 stressp_1 (:,:,iblk), stressp_2 (:,:,iblk), &
00391 stressp_3 (:,:,iblk), stressp_4 (:,:,iblk), &
00392 stressm_1 (:,:,iblk), stressm_2 (:,:,iblk), &
00393 stressm_3 (:,:,iblk), stressm_4 (:,:,iblk), &
00394 stress12_1(:,:,iblk), stress12_2(:,:,iblk), &
00395 stress12_3(:,:,iblk), stress12_4(:,:,iblk), &
00396 shear (:,:,iblk), divu (:,:,iblk), &
00397 prs_sig (:,:,iblk), &
00398 rdg_conv (:,:,iblk), rdg_shear (:,:,iblk), &
00399 str (:,:,:) )
00400
00401
00402
00403
00404
00405
00406 call stepu (nx_block, ny_block, &
00407 icellu (iblk), &
00408 indxui (:,iblk), indxuj (:,iblk), &
00409 aiu (:,:,iblk), str (:,:,:), &
00410 uocn (:,:,iblk), vocn (:,:,iblk), &
00411 waterx (:,:,iblk), watery (:,:,iblk), &
00412 forcex (:,:,iblk), forcey (:,:,iblk), &
00413 umassdtei(:,:,iblk), fm (:,:,iblk), &
00414 uarear (:,:,iblk), &
00415 strocnx (:,:,iblk), strocny (:,:,iblk), &
00416 strintx (:,:,iblk), strinty (:,:,iblk), &
00417 uvel (:,:,iblk), vvel (:,:,iblk))
00418
00419 enddo
00420
00421
00422
00423
00424
00425
00426
00427 do iblk = 1,nblocks
00428 fld2(1:nx_block,1:ny_block,1,iblk) = uvel(1:nx_block,1:ny_block,iblk)
00429 fld2(1:nx_block,1:ny_block,2,iblk) = vvel(1:nx_block,1:ny_block,iblk)
00430 enddo
00431
00432
00433 call ice_timer_start(timer_bound)
00434 call ice_HaloUpdate (fld2, halo_info, &
00435 field_loc_NEcorner, field_type_vector)
00436
00437
00438
00439
00440 call ice_timer_stop(timer_bound)
00441
00442
00443 do iblk = 1,nblocks
00444 uvel(1:nx_block,1:ny_block,iblk) = fld2(1:nx_block,1:ny_block,1,iblk)
00445 vvel(1:nx_block,1:ny_block,iblk) = fld2(1:nx_block,1:ny_block,2,iblk)
00446 enddo
00447
00448
00449
00450
00451 enddo
00452
00453 deallocate(fld2)
00454
00455
00456 if (trim(grid_type) == 'tripole') then
00457
00458 call ice_HaloUpdate_stress(stressp_1, stressp_3, halo_info, &
00459 field_loc_center, field_type_scalar)
00460 call ice_HaloUpdate_stress(stressp_3, stressp_1, halo_info, &
00461 field_loc_center, field_type_scalar)
00462 call ice_HaloUpdate_stress(stressp_2, stressp_4, halo_info, &
00463 field_loc_center, field_type_scalar)
00464 call ice_HaloUpdate_stress(stressp_4, stressp_2, halo_info, &
00465 field_loc_center, field_type_scalar)
00466
00467 call ice_HaloUpdate_stress(stressm_1, stressm_3, halo_info, &
00468 field_loc_center, field_type_scalar)
00469 call ice_HaloUpdate_stress(stressm_3, stressm_1, halo_info, &
00470 field_loc_center, field_type_scalar)
00471 call ice_HaloUpdate_stress(stressm_2, stressm_4, halo_info, &
00472 field_loc_center, field_type_scalar)
00473 call ice_HaloUpdate_stress(stressm_4, stressm_2, halo_info, &
00474 field_loc_center, field_type_scalar)
00475
00476 call ice_HaloUpdate_stress(stress12_1, stress12_3, halo_info, &
00477 field_loc_center, field_type_scalar)
00478 call ice_HaloUpdate_stress(stress12_3, stress12_1, halo_info, &
00479 field_loc_center, field_type_scalar)
00480 call ice_HaloUpdate_stress(stress12_2, stress12_4, halo_info, &
00481 field_loc_center, field_type_scalar)
00482 call ice_HaloUpdate_stress(stress12_4, stress12_2, halo_info, &
00483 field_loc_center, field_type_scalar)
00484
00485 endif
00486
00487
00488
00489
00490
00491
00492
00493
00494
00495 do iblk = 1, nblocks
00496
00497 call evp_finish &
00498 (nx_block, ny_block, &
00499 icellu (iblk), &
00500 indxui (:,iblk), indxuj (:,iblk), &
00501 uvel (:,:,iblk), vvel (:,:,iblk), &
00502 uocn (:,:,iblk), vocn (:,:,iblk), &
00503 aiu (:,:,iblk), &
00504 strocnx (:,:,iblk), strocny (:,:,iblk), &
00505 strocnxT(:,:,iblk), strocnyT(:,:,iblk))
00506
00507 enddo
00508
00509
00510 call u2tgrid_vector(strocnxT)
00511 call u2tgrid_vector(strocnyT)
00512
00513 call ice_timer_stop(timer_dynamics)
00514
00515
00516 end subroutine evp
00517
00518
00519
00520
00521
00522
00523
00524
00525 subroutine init_evp (dt)
00526
00527
00528
00529
00530
00531
00532
00533
00534
00535
00536
00537 use ice_boundary
00538 use ice_blocks
00539 use ice_domain
00540 use ice_state
00541 use ice_flux
00542 use ice_grid
00543 use ice_fileunits
00544
00545
00546
00547 real (kind=dbl_kind), intent(in) ::
00548 dt
00549
00550
00551
00552 integer (kind=int_kind) ::
00553 i, j, k,
00554 iblk
00555
00556 real (kind=dbl_kind) ::
00557 dte ,
00558 ecc ,
00559 tdamp2
00560
00561 call set_evp_parameters (dt)
00562
00563 if (my_task == master_task) then
00564 write(nu_diag,*) 'dt_dyn = ',dt
00565 write(nu_diag,*) 'dte = ',dt/real(ndte,kind=dbl_kind)
00566 write(nu_diag,*) 'tdamp =', eyc*dt
00567 endif
00568
00569 allocate(fcor_blk(nx_block,ny_block,max_blocks))
00570
00571
00572 do iblk = 1, nblocks
00573 do j = 1, ny_block
00574 do i = 1, nx_block
00575
00576
00577 uvel(i,j,iblk) = c0
00578 vvel(i,j,iblk) = c0
00579
00580
00581 divu (i,j,iblk) = c0
00582 shear(i,j,iblk) = c0
00583 rdg_conv (i,j,iblk) = c0
00584 rdg_shear(i,j,iblk) = c0
00585
00586
00587
00588 fcor_blk(i,j,iblk) = c2*omega*sin(ULAT(i,j,iblk))
00589
00590
00591 stressp_1 (i,j,iblk) = c0
00592 stressp_2 (i,j,iblk) = c0
00593 stressp_3 (i,j,iblk) = c0
00594 stressp_4 (i,j,iblk) = c0
00595 stressm_1 (i,j,iblk) = c0
00596 stressm_2 (i,j,iblk) = c0
00597 stressm_3 (i,j,iblk) = c0
00598 stressm_4 (i,j,iblk) = c0
00599 stress12_1(i,j,iblk) = c0
00600 stress12_2(i,j,iblk) = c0
00601 stress12_3(i,j,iblk) = c0
00602 stress12_4(i,j,iblk) = c0
00603
00604
00605 iceumask(i,j,iblk) = .false.
00606
00607 enddo
00608 enddo
00609 enddo
00610
00611
00612 end subroutine init_evp
00613
00614
00615
00616
00617
00618
00619
00620
00621 subroutine set_evp_parameters (dt)
00622
00623
00624
00625
00626
00627
00628
00629
00630
00631
00632
00633
00634
00635
00636
00637
00638 real (kind=dbl_kind), intent(in) ::
00639 dt
00640
00641
00642
00643 real (kind=dbl_kind) ::
00644 dte ,
00645 ecc ,
00646 tdamp2
00647
00648
00649 dte = dt/real(ndte,kind=dbl_kind)
00650 dtei = c1/dte
00651
00652
00653 ecc = c4
00654 ecci = p25
00655
00656
00657 tdamp2 = c2*eyc*dt
00658 dte2T = dte/tdamp2
00659 denom1 = c1/(c1+dte2T)
00660 denom2 = c1/(c1+dte2T*ecc)
00661 rcon = 1230._dbl_kind*eyc*dt*dtei**2
00662
00663 end subroutine set_evp_parameters
00664
00665
00666
00667
00668
00669
00670
00671
00672 subroutine evp_prep1 (nx_block, ny_block, &
00673 ilo, ihi, jlo, jhi, &
00674 aice, vice, &
00675 vsno, tmask, &
00676 strairxT, strairyT, &
00677 strairx, strairy, &
00678 tmass, icetmask)
00679
00680
00681
00682
00683
00684
00685
00686
00687
00688
00689
00690
00691
00692
00693
00694
00695 integer (kind=int_kind), intent(in) ::
00696 nx_block, ny_block,
00697 ilo,ihi,jlo,jhi
00698
00699 real (kind=dbl_kind), dimension (nx_block,ny_block),
00700 intent(in) ::
00701 aice ,
00702 vice ,
00703 vsno ,
00704 strairxT,
00705 strairyT
00706
00707 logical (kind=log_kind), dimension (nx_block,ny_block),
00708 intent(in) ::
00709 tmask
00710
00711 real (kind=dbl_kind), dimension (nx_block,ny_block),
00712 intent(out) ::
00713 strairx ,
00714 strairy ,
00715 tmass
00716
00717 integer (kind=int_kind), dimension (nx_block,ny_block),
00718 intent(out) ::
00719 icetmask
00720
00721
00722
00723 integer (kind=int_kind) ::
00724 i, j
00725
00726 logical (kind=log_kind), dimension(nx_block,ny_block) ::
00727 tmphm
00728
00729 do j = 1, ny_block
00730 do i = 1, nx_block
00731
00732
00733
00734
00735
00736
00737 if (tmask(i,j)) then
00738 tmass(i,j) = (rhoi*vice(i,j) + rhos*vsno(i,j))
00739 else
00740 tmass(i,j) = c0
00741 endif
00742
00743
00744
00745
00746 tmphm(i,j) = tmask(i,j) .and. (aice (i,j) > a_min) &
00747 .and. (tmass(i,j) > m_min)
00748
00749
00750
00751
00752
00753
00754 strairx(i,j) = strairxT(i,j)
00755 strairy(i,j) = strairyT(i,j)
00756
00757
00758
00759
00760 icetmask (i,j) = 0
00761
00762 enddo
00763 enddo
00764
00765 do j = jlo, jhi
00766 do i = ilo, ihi
00767
00768
00769 if (tmphm(i-1,j+1) .or. tmphm(i,j+1) .or. tmphm(i+1,j+1) .or. &
00770 tmphm(i-1,j) .or. tmphm(i,j) .or. tmphm(i+1,j) .or. &
00771 tmphm(i-1,j-1) .or. tmphm(i,j-1) .or. tmphm(i+1,j-1) ) then
00772 icetmask(i,j) = 1
00773 endif
00774
00775 if (.not.tmask(i,j)) icetmask(i,j) = 0
00776
00777 enddo
00778 enddo
00779
00780 end subroutine evp_prep1
00781
00782
00783
00784
00785
00786
00787
00788
00789 subroutine evp_prep2 (nx_block, ny_block, &
00790 ilo, ihi, jlo, jhi, &
00791 icellt, icellu, &
00792 indxti, indxtj, &
00793 indxui, indxuj, &
00794 aiu, umass, &
00795 umassdtei, fcor, &
00796 umask, &
00797 uocn, vocn, &
00798 strairx, strairy, &
00799 ss_tltx, ss_tlty, &
00800 icetmask, iceumask, &
00801 fm, &
00802 strtltx, strtlty, &
00803 strocnx, strocny, &
00804 strintx, strinty, &
00805 waterx, watery, &
00806 forcex, forcey, &
00807 stressp_1, stressp_2, &
00808 stressp_3, stressp_4, &
00809 stressm_1, stressm_2, &
00810 stressm_3, stressm_4, &
00811 stress12_1, stress12_2, &
00812 stress12_3, stress12_4, &
00813 uvel, vvel)
00814
00815
00816
00817
00818
00819
00820
00821
00822
00823
00824
00825
00826
00827
00828
00829
00830
00831
00832 integer (kind=int_kind), intent(in) ::
00833 nx_block, ny_block,
00834 ilo,ihi,jlo,jhi
00835
00836 integer (kind=int_kind), intent(out) ::
00837 icellt ,
00838 icellu
00839
00840 integer (kind=int_kind), dimension (nx_block*ny_block),
00841 intent(out) ::
00842 indxti ,
00843 indxtj ,
00844 indxui ,
00845 indxuj
00846
00847 logical (kind=log_kind), dimension (nx_block,ny_block),
00848 intent(in) ::
00849 umask
00850
00851 integer (kind=int_kind), dimension (nx_block,ny_block),
00852 intent(in) ::
00853 icetmask
00854
00855 logical (kind=log_kind), dimension (nx_block,ny_block),
00856 intent(inout) ::
00857 iceumask
00858
00859 real (kind=dbl_kind), dimension (nx_block,ny_block), intent(in) ::
00860 aiu ,
00861 umass ,
00862 fcor ,
00863 strairx ,
00864 strairy ,
00865 uocn ,
00866 vocn ,
00867 ss_tltx ,
00868 ss_tlty
00869
00870 real (kind=dbl_kind), dimension (nx_block,ny_block),
00871 intent(out) ::
00872 umassdtei,
00873 waterx ,
00874 watery ,
00875 forcex ,
00876 forcey
00877
00878 real (kind=dbl_kind), dimension (nx_block,ny_block),
00879 intent(inout) ::
00880 fm ,
00881 stressp_1, stressp_2, stressp_3, stressp_4 ,
00882 stressm_1, stressm_2, stressm_3, stressm_4 ,
00883 stress12_1,stress12_2,stress12_3,stress12_4,
00884 uvel ,
00885 vvel ,
00886 strtltx ,
00887 strtlty ,
00888 strocnx ,
00889 strocny ,
00890 strintx ,
00891 strinty
00892
00893
00894
00895 integer (kind=int_kind) ::
00896 i, j, ij
00897
00898 logical (kind=log_kind), dimension(nx_block,ny_block) ::
00899 iceumask_old
00900
00901
00902
00903
00904
00905 do j = 1, ny_block
00906 do i = 1, nx_block
00907 waterx (i,j) = c0
00908 watery (i,j) = c0
00909 forcex (i,j) = c0
00910 forcey (i,j) = c0
00911 umassdtei(i,j) = c0
00912
00913 if (icetmask(i,j)==0) then
00914 stressp_1 (i,j) = c0
00915 stressp_2 (i,j) = c0
00916 stressp_3 (i,j) = c0
00917 stressp_4 (i,j) = c0
00918 stressm_1 (i,j) = c0
00919 stressm_2 (i,j) = c0
00920 stressm_3 (i,j) = c0
00921 stressm_4 (i,j) = c0
00922 stress12_1(i,j) = c0
00923 stress12_2(i,j) = c0
00924 stress12_3(i,j) = c0
00925 stress12_4(i,j) = c0
00926 endif
00927 enddo
00928 enddo
00929
00930
00931
00932
00933
00934
00935
00936 icellt = 0
00937 do j = jlo, jhi+1
00938 do i = ilo, ihi+1
00939 if (icetmask(i,j) == 1) then
00940 icellt = icellt + 1
00941 indxti(icellt) = i
00942 indxtj(icellt) = j
00943 endif
00944 enddo
00945 enddo
00946
00947
00948
00949
00950
00951
00952
00953 icellu = 0
00954 do j = jlo, jhi
00955 do i = ilo, ihi
00956
00957
00958 iceumask_old(i,j) = iceumask(i,j)
00959 iceumask(i,j) = (umask(i,j)) .and. (aiu (i,j) > a_min) &
00960 .and. (umass(i,j) > m_min)
00961
00962 if (iceumask(i,j)) then
00963 icellu = icellu + 1
00964 indxui(icellu) = i
00965 indxuj(icellu) = j
00966
00967
00968 if (.not. iceumask_old(i,j)) then
00969 uvel(i,j) = uocn(i,j)
00970 vvel(i,j) = vocn(i,j)
00971 endif
00972 else
00973
00974 uvel(i,j) = c0
00975 vvel(i,j) = c0
00976 strintx(i,j) = c0
00977 strinty(i,j) = c0
00978 strocnx(i,j) = c0
00979 strocny(i,j) = c0
00980 endif
00981 enddo
00982 enddo
00983
00984
00985
00986
00987
00988 do ij = 1, icellu
00989 i = indxui(ij)
00990 j = indxuj(ij)
00991
00992 umassdtei(i,j) = umass(i,j)*dtei
00993 fm(i,j) = fcor(i,j)*umass(i,j)
00994
00995
00996 waterx(i,j) = uocn(i,j)*cosw - vocn(i,j)*sinw
00997 watery(i,j) = vocn(i,j)*cosw + uocn(i,j)*sinw
00998
00999
01000 #ifndef coupled
01001
01002 strtltx(i,j) = -fm(i,j)*vocn(i,j)
01003 strtlty(i,j) = fm(i,j)*uocn(i,j)
01004 #else
01005 strtltx(i,j) = -gravit*umass(i,j)*ss_tltx(i,j)
01006 strtlty(i,j) = -gravit*umass(i,j)*ss_tlty(i,j)
01007 #endif
01008 forcex(i,j) = strairx(i,j) + strtltx(i,j)
01009 forcey(i,j) = strairy(i,j) + strtlty(i,j)
01010 enddo
01011
01012 end subroutine evp_prep2
01013
01014
01015
01016
01017
01018
01019
01020
01021 subroutine stress (nx_block, ny_block, &
01022 ksub, icellt, &
01023 indxti, indxtj, &
01024 uvel, vvel, &
01025 dxt, dyt, &
01026 dxhy, dyhx, &
01027 cxp, cyp, &
01028 cxm, cym, &
01029 tarear, tinyarea, &
01030 strength, &
01031 stressp_1, stressp_2, &
01032 stressp_3, stressp_4, &
01033 stressm_1, stressm_2, &
01034 stressm_3, stressm_4, &
01035 stress12_1, stress12_2, &
01036 stress12_3, stress12_4, &
01037 shear, divu, &
01038 prs_sig, &
01039 rdg_conv, rdg_shear, &
01040 str )
01041
01042
01043
01044
01045
01046
01047
01048
01049
01050
01051
01052
01053
01054
01055
01056 integer (kind=int_kind), intent(in) ::
01057 nx_block, ny_block,
01058 ksub ,
01059 icellt
01060
01061 integer (kind=int_kind), dimension (nx_block*ny_block),
01062 intent(in) ::
01063 indxti ,
01064 indxtj
01065
01066 real (kind=dbl_kind), dimension (nx_block,ny_block), intent(in) ::
01067 strength ,
01068 uvel ,
01069 vvel ,
01070 dxt ,
01071 dyt ,
01072 dxhy ,
01073 dyhx ,
01074 cyp ,
01075 cxp ,
01076 cym ,
01077 cxm ,
01078 tarear ,
01079 tinyarea
01080
01081 real (kind=dbl_kind), dimension (nx_block,ny_block),
01082 intent(inout) ::
01083 stressp_1, stressp_2, stressp_3, stressp_4 ,
01084 stressm_1, stressm_2, stressm_3, stressm_4 ,
01085 stress12_1,stress12_2,stress12_3,stress12_4
01086
01087 real (kind=dbl_kind), dimension (nx_block,ny_block),
01088 intent(inout) ::
01089 prs_sig ,
01090 shear ,
01091 divu ,
01092 rdg_conv ,
01093 rdg_shear
01094
01095 real (kind=dbl_kind), dimension(nx_block,ny_block,8),
01096 intent(out) ::
01097 str
01098
01099
01100
01101 integer (kind=int_kind) ::
01102 i, j, ij
01103
01104 real (kind=dbl_kind) ::
01105 divune, divunw, divuse, divusw ,
01106 tensionne, tensionnw, tensionse, tensionsw,
01107 shearne, shearnw, shearse, shearsw ,
01108 Deltane, Deltanw, Deltase, Deltasw ,
01109 c0ne, c0nw, c0se, c0sw ,
01110 c1ne, c1nw, c1se, c1sw ,
01111 ssigpn, ssigps, ssigpe, ssigpw ,
01112 ssigmn, ssigms, ssigme, ssigmw ,
01113 ssig12n, ssig12s, ssig12e, ssig12w ,
01114 ssigp1, ssigp2, ssigm1, ssigm2, ssig121, ssig122,
01115 csigpne, csigpnw, csigpse, csigpsw ,
01116 csigmne, csigmnw, csigmse, csigmsw ,
01117 csig12ne, csig12nw, csig12se, csig12sw ,
01118 str12ew, str12we, str12ns, str12sn ,
01119 strp_tmp, strm_tmp, tmp
01120
01121
01122
01123
01124
01125 str(:,:,:) = c0
01126
01127
01128
01129
01130 do ij = 1, icellt
01131 i = indxti(ij)
01132 j = indxtj(ij)
01133
01134
01135
01136
01137
01138
01139 divune = cyp(i,j)*uvel(i ,j ) - dyt(i,j)*uvel(i-1,j ) &
01140 + cxp(i,j)*vvel(i ,j ) - dxt(i,j)*vvel(i ,j-1)
01141 divunw = cym(i,j)*uvel(i-1,j ) + dyt(i,j)*uvel(i ,j ) &
01142 + cxp(i,j)*vvel(i-1,j ) - dxt(i,j)*vvel(i-1,j-1)
01143 divusw = cym(i,j)*uvel(i-1,j-1) + dyt(i,j)*uvel(i ,j-1) &
01144 + cxm(i,j)*vvel(i-1,j-1) + dxt(i,j)*vvel(i-1,j )
01145 divuse = cyp(i,j)*uvel(i ,j-1) - dyt(i,j)*uvel(i-1,j-1) &
01146 + cxm(i,j)*vvel(i ,j-1) + dxt(i,j)*vvel(i ,j )
01147
01148
01149 tensionne = -cym(i,j)*uvel(i ,j ) - dyt(i,j)*uvel(i-1,j ) &
01150 + cxm(i,j)*vvel(i ,j ) + dxt(i,j)*vvel(i ,j-1)
01151 tensionnw = -cyp(i,j)*uvel(i-1,j ) + dyt(i,j)*uvel(i ,j ) &
01152 + cxm(i,j)*vvel(i-1,j ) + dxt(i,j)*vvel(i-1,j-1)
01153 tensionsw = -cyp(i,j)*uvel(i-1,j-1) + dyt(i,j)*uvel(i ,j-1) &
01154 + cxp(i,j)*vvel(i-1,j-1) - dxt(i,j)*vvel(i-1,j )
01155 tensionse = -cym(i,j)*uvel(i ,j-1) - dyt(i,j)*uvel(i-1,j-1) &
01156 + cxp(i,j)*vvel(i ,j-1) - dxt(i,j)*vvel(i ,j )
01157
01158
01159 shearne = -cym(i,j)*vvel(i ,j ) - dyt(i,j)*vvel(i-1,j ) &
01160 - cxm(i,j)*uvel(i ,j ) - dxt(i,j)*uvel(i ,j-1)
01161 shearnw = -cyp(i,j)*vvel(i-1,j ) + dyt(i,j)*vvel(i ,j ) &
01162 - cxm(i,j)*uvel(i-1,j ) - dxt(i,j)*uvel(i-1,j-1)
01163 shearsw = -cyp(i,j)*vvel(i-1,j-1) + dyt(i,j)*vvel(i ,j-1) &
01164 - cxp(i,j)*uvel(i-1,j-1) + dxt(i,j)*uvel(i-1,j )
01165 shearse = -cym(i,j)*vvel(i ,j-1) - dyt(i,j)*vvel(i-1,j-1) &
01166 - cxp(i,j)*uvel(i ,j-1) + dxt(i,j)*uvel(i ,j )
01167
01168
01169 Deltane = sqrt(divune**2 + ecci*(tensionne**2 + shearne**2))
01170 Deltanw = sqrt(divunw**2 + ecci*(tensionnw**2 + shearnw**2))
01171 Deltase = sqrt(divuse**2 + ecci*(tensionse**2 + shearse**2))
01172 Deltasw = sqrt(divusw**2 + ecci*(tensionsw**2 + shearsw**2))
01173
01174
01175
01176
01177 if (ksub == ndte) then
01178 divu(i,j) = p25*(divune + divunw + divuse + divusw) * tarear(i,j)
01179 tmp = p25*(Deltane + Deltanw + Deltase + Deltasw) * tarear(i,j)
01180 rdg_conv(i,j) = -min(divu(i,j),c0)
01181 rdg_shear(i,j) = p5*(tmp-abs(divu(i,j)))
01182
01183
01184
01185 shear(i,j) = p25*tarear(i,j)*sqrt( &
01186 (tensionne + tensionnw + tensionse + tensionsw)**2 &
01187 + (shearne + shearnw + shearse + shearsw)**2)
01188
01189 endif
01190
01191
01192
01193
01194
01195 if (evp_damping) then
01196
01197 c0ne = min(strength(i,j)/max(Deltane,c4*tinyarea(i,j)),rcon)
01198 c0nw = min(strength(i,j)/max(Deltanw,c4*tinyarea(i,j)),rcon)
01199 c0sw = min(strength(i,j)/max(Deltasw,c4*tinyarea(i,j)),rcon)
01200 c0se = min(strength(i,j)/max(Deltase,c4*tinyarea(i,j)),rcon)
01201 prs_sig(i,j) = strength(i,j)* &
01202 Deltane/max(Deltane,c4*tinyarea(i,j))
01203 else
01204
01205 c0ne = strength(i,j)/max(Deltane,tinyarea(i,j))
01206 c0nw = strength(i,j)/max(Deltanw,tinyarea(i,j))
01207 c0sw = strength(i,j)/max(Deltasw,tinyarea(i,j))
01208 c0se = strength(i,j)/max(Deltase,tinyarea(i,j))
01209 prs_sig(i,j) = c0ne*Deltane
01210 endif
01211
01212 c1ne = c0ne*dte2T
01213 c1nw = c0nw*dte2T
01214 c1sw = c0sw*dte2T
01215 c1se = c0se*dte2T
01216
01217
01218
01219
01220
01221
01222 stressp_1(i,j) = (stressp_1(i,j) + c1ne*(divune - Deltane)) &
01223 * denom1
01224 stressp_2(i,j) = (stressp_2(i,j) + c1nw*(divunw - Deltanw)) &
01225 * denom1
01226 stressp_3(i,j) = (stressp_3(i,j) + c1sw*(divusw - Deltasw)) &
01227 * denom1
01228 stressp_4(i,j) = (stressp_4(i,j) + c1se*(divuse - Deltase)) &
01229 * denom1
01230
01231 stressm_1(i,j) = (stressm_1(i,j) + c1ne*tensionne) * denom2
01232 stressm_2(i,j) = (stressm_2(i,j) + c1nw*tensionnw) * denom2
01233 stressm_3(i,j) = (stressm_3(i,j) + c1sw*tensionsw) * denom2
01234 stressm_4(i,j) = (stressm_4(i,j) + c1se*tensionse) * denom2
01235
01236 stress12_1(i,j) = (stress12_1(i,j) + c1ne*shearne*p5) * denom2
01237 stress12_2(i,j) = (stress12_2(i,j) + c1nw*shearnw*p5) * denom2
01238 stress12_3(i,j) = (stress12_3(i,j) + c1sw*shearsw*p5) * denom2
01239 stress12_4(i,j) = (stress12_4(i,j) + c1se*shearse*p5) * denom2
01240
01241
01242
01243
01244
01245
01246
01247
01248
01249
01250
01251
01252
01253
01254
01255
01256
01257
01258
01259
01260
01261
01262
01263
01264
01265
01266
01267
01268
01269
01270 ssigpn = stressp_1(i,j) + stressp_2(i,j)
01271 ssigps = stressp_3(i,j) + stressp_4(i,j)
01272 ssigpe = stressp_1(i,j) + stressp_4(i,j)
01273 ssigpw = stressp_2(i,j) + stressp_3(i,j)
01274 ssigp1 =(stressp_1(i,j) + stressp_3(i,j))*p055
01275 ssigp2 =(stressp_2(i,j) + stressp_4(i,j))*p055
01276
01277 ssigmn = stressm_1(i,j) + stressm_2(i,j)
01278 ssigms = stressm_3(i,j) + stressm_4(i,j)
01279 ssigme = stressm_1(i,j) + stressm_4(i,j)
01280 ssigmw = stressm_2(i,j) + stressm_3(i,j)
01281 ssigm1 =(stressm_1(i,j) + stressm_3(i,j))*p055
01282 ssigm2 =(stressm_2(i,j) + stressm_4(i,j))*p055
01283
01284 ssig12n = stress12_1(i,j) + stress12_2(i,j)
01285 ssig12s = stress12_3(i,j) + stress12_4(i,j)
01286 ssig12e = stress12_1(i,j) + stress12_4(i,j)
01287 ssig12w = stress12_2(i,j) + stress12_3(i,j)
01288 ssig121 =(stress12_1(i,j) + stress12_3(i,j))*p111
01289 ssig122 =(stress12_2(i,j) + stress12_4(i,j))*p111
01290
01291 csigpne = p111*stressp_1(i,j) + ssigp2 + p027*stressp_3(i,j)
01292 csigpnw = p111*stressp_2(i,j) + ssigp1 + p027*stressp_4(i,j)
01293 csigpsw = p111*stressp_3(i,j) + ssigp2 + p027*stressp_1(i,j)
01294 csigpse = p111*stressp_4(i,j) + ssigp1 + p027*stressp_2(i,j)
01295
01296 csigmne = p111*stressm_1(i,j) + ssigm2 + p027*stressm_3(i,j)
01297 csigmnw = p111*stressm_2(i,j) + ssigm1 + p027*stressm_4(i,j)
01298 csigmsw = p111*stressm_3(i,j) + ssigm2 + p027*stressm_1(i,j)
01299 csigmse = p111*stressm_4(i,j) + ssigm1 + p027*stressm_2(i,j)
01300
01301 csig12ne = p222*stress12_1(i,j) + ssig122 &
01302 + p055*stress12_3(i,j)
01303 csig12nw = p222*stress12_2(i,j) + ssig121 &
01304 + p055*stress12_4(i,j)
01305 csig12sw = p222*stress12_3(i,j) + ssig122 &
01306 + p055*stress12_1(i,j)
01307 csig12se = p222*stress12_4(i,j) + ssig121 &
01308 + p055*stress12_2(i,j)
01309
01310 str12ew = p5*dxt(i,j)*(p333*ssig12e + p166*ssig12w)
01311 str12we = p5*dxt(i,j)*(p333*ssig12w + p166*ssig12e)
01312 str12ns = p5*dyt(i,j)*(p333*ssig12n + p166*ssig12s)
01313 str12sn = p5*dyt(i,j)*(p333*ssig12s + p166*ssig12n)
01314
01315
01316
01317
01318 strp_tmp = p25*dyt(i,j)*(p333*ssigpn + p166*ssigps)
01319 strm_tmp = p25*dyt(i,j)*(p333*ssigmn + p166*ssigms)
01320
01321
01322 str(i,j,1) = -strp_tmp - strm_tmp - str12ew &
01323 + dxhy(i,j)*(-csigpne + csigmne) + dyhx(i,j)*csig12ne
01324
01325
01326 str(i,j,2) = strp_tmp + strm_tmp - str12we &
01327 + dxhy(i,j)*(-csigpnw + csigmnw) + dyhx(i,j)*csig12nw
01328
01329 strp_tmp = p25*dyt(i,j)*(p333*ssigps + p166*ssigpn)
01330 strm_tmp = p25*dyt(i,j)*(p333*ssigms + p166*ssigmn)
01331
01332
01333 str(i,j,3) = -strp_tmp - strm_tmp + str12ew &
01334 + dxhy(i,j)*(-csigpse + csigmse) + dyhx(i,j)*csig12se
01335
01336
01337 str(i,j,4) = strp_tmp + strm_tmp + str12we &
01338 + dxhy(i,j)*(-csigpsw + csigmsw) + dyhx(i,j)*csig12sw
01339
01340
01341
01342
01343 strp_tmp = p25*dxt(i,j)*(p333*ssigpe + p166*ssigpw)
01344 strm_tmp = p25*dxt(i,j)*(p333*ssigme + p166*ssigmw)
01345
01346
01347 str(i,j,5) = -strp_tmp + strm_tmp - str12ns &
01348 - dyhx(i,j)*(csigpne + csigmne) + dxhy(i,j)*csig12ne
01349
01350
01351 str(i,j,6) = strp_tmp - strm_tmp - str12sn &
01352 - dyhx(i,j)*(csigpse + csigmse) + dxhy(i,j)*csig12se
01353
01354 strp_tmp = p25*dxt(i,j)*(p333*ssigpw + p166*ssigpe)
01355 strm_tmp = p25*dxt(i,j)*(p333*ssigmw + p166*ssigme)
01356
01357
01358 str(i,j,7) = -strp_tmp + strm_tmp + str12ns &
01359 - dyhx(i,j)*(csigpnw + csigmnw) + dxhy(i,j)*csig12nw
01360
01361
01362 str(i,j,8) = strp_tmp - strm_tmp + str12sn &
01363 - dyhx(i,j)*(csigpsw + csigmsw) + dxhy(i,j)*csig12sw
01364
01365 enddo
01366
01367 end subroutine stress
01368
01369
01370
01371
01372
01373
01374
01375
01376 subroutine stepu (nx_block, ny_block, &
01377 icellu, &
01378 indxui, indxuj, &
01379 aiu, str, &
01380 uocn, vocn, &
01381 waterx, watery, &
01382 forcex, forcey, &
01383 umassdtei, fm, &
01384 uarear, &
01385 strocnx, strocny, &
01386 strintx, strinty, &
01387 uvel, vvel)
01388
01389
01390
01391
01392
01393
01394
01395
01396
01397
01398
01399
01400
01401
01402 integer (kind=int_kind), intent(in) ::
01403 nx_block, ny_block,
01404 icellu
01405
01406 integer (kind=int_kind), dimension (nx_block*ny_block),
01407 intent(in) ::
01408 indxui ,
01409 indxuj
01410
01411 real (kind=dbl_kind), dimension (nx_block,ny_block), intent(in) ::
01412 aiu ,
01413 waterx ,
01414 watery ,
01415 forcex ,
01416 forcey ,
01417 umassdtei,
01418 uocn ,
01419 vocn ,
01420 fm ,
01421 uarear
01422
01423 real (kind=dbl_kind), dimension(nx_block,ny_block,8),
01424 intent(in) ::
01425 str
01426
01427 real (kind=dbl_kind), dimension (nx_block,ny_block),
01428 intent(inout) ::
01429 uvel ,
01430 vvel
01431
01432 real (kind=dbl_kind), dimension (nx_block,ny_block),
01433 intent(inout) ::
01434 strocnx ,
01435 strocny ,
01436 strintx ,
01437 strinty
01438
01439
01440
01441 integer (kind=int_kind) ::
01442 i, j, ij
01443
01444 real (kind=dbl_kind) ::
01445 uold, vold ,
01446 vrel ,
01447 cca,ccb,ab2,cc1,cc2,
01448 taux, tauy
01449
01450
01451
01452
01453
01454 do ij =1, icellu
01455 i = indxui(ij)
01456 j = indxuj(ij)
01457
01458 uold = uvel(i,j)
01459 vold = vvel(i,j)
01460
01461
01462 vrel = aiu(i,j)*dragw*sqrt((uocn(i,j) - uold)**2 + &
01463 (vocn(i,j) - vold)**2)
01464
01465 taux = vrel*waterx(i,j)
01466 tauy = vrel*watery(i,j)
01467
01468
01469 cca = umassdtei(i,j) + vrel * cosw
01470 ccb = fm(i,j) + vrel * sinw
01471 ab2 = cca**2 + ccb**2
01472
01473
01474 strintx(i,j) = uarear(i,j)* &
01475 (str(i,j,1) + str(i+1,j,2) + str(i,j+1,3) + str(i+1,j+1,4))
01476 strinty(i,j) = uarear(i,j)* &
01477 (str(i,j,5) + str(i,j+1,6) + str(i+1,j,7) + str(i+1,j+1,8))
01478
01479
01480 cc1 = strintx(i,j) + forcex(i,j) + taux &
01481 + umassdtei(i,j)*uold
01482 cc2 = strinty(i,j) + forcey(i,j) + tauy &
01483 + umassdtei(i,j)*vold
01484
01485 uvel(i,j) = (cca*cc1 + ccb*cc2) / ab2
01486 vvel(i,j) = (cca*cc2 - ccb*cc1) / ab2
01487
01488
01489
01490
01491
01492 strocnx(i,j) = taux
01493 strocny(i,j) = tauy
01494
01495 enddo
01496
01497 end subroutine stepu
01498
01499
01500
01501
01502
01503
01504
01505
01506 subroutine evp_finish (nx_block, ny_block, &
01507 icellu, &
01508 indxui, indxuj, &
01509 uvel, vvel, &
01510 uocn, vocn, &
01511 aiu, &
01512 strocnx, strocny, &
01513 strocnxT, strocnyT)
01514
01515
01516
01517
01518
01519
01520
01521
01522
01523
01524
01525
01526
01527
01528 integer (kind=int_kind), intent(in) ::
01529 nx_block, ny_block,
01530 icellu
01531
01532 integer (kind=int_kind), dimension (nx_block*ny_block),
01533 intent(in) ::
01534 indxui ,
01535 indxuj
01536
01537 real (kind=dbl_kind), dimension (nx_block,ny_block), intent(in) ::
01538 uvel ,
01539 vvel ,
01540 uocn ,
01541 vocn ,
01542 aiu
01543
01544 real (kind=dbl_kind), dimension (nx_block,ny_block),
01545 intent(inout) ::
01546 strocnx ,
01547 strocny ,
01548 strocnxT,
01549 strocnyT
01550
01551
01552
01553 integer (kind=int_kind) ::
01554 i, j, ij
01555
01556 real (kind=dbl_kind) :: vrel
01557
01558 do j = 1, ny_block
01559 do i = 1, nx_block
01560 strocnxT(i,j) = c0
01561 strocnyT(i,j) = c0
01562 enddo
01563 enddo
01564
01565
01566 do ij =1, icellu
01567 i = indxui(ij)
01568 j = indxuj(ij)
01569
01570 vrel = dragw*sqrt((uocn(i,j) - uvel(i,j))**2 + &
01571 (vocn(i,j) - vvel(i,j))**2)
01572 strocnx(i,j) = strocnx(i,j) &
01573 - vrel*(uvel(i,j)*cosw - vvel(i,j)*sinw) * aiu(i,j)
01574 strocny(i,j) = strocny(i,j) &
01575 - vrel*(vvel(i,j)*cosw + uvel(i,j)*sinw) * aiu(i,j)
01576
01577
01578 strocnxT(i,j) = strocnx(i,j) / aiu(i,j)
01579 strocnyT(i,j) = strocny(i,j) / aiu(i,j)
01580 enddo
01581
01582 end subroutine evp_finish
01583
01584
01585
01586
01587
01588
01589
01590
01591 subroutine principal_stress(nx_block, ny_block, &
01592 stressp_1, stressm_1, &
01593 stress12_1, prs_sig, &
01594 sig1, sig2)
01595
01596
01597
01598
01599
01600
01601
01602
01603
01604
01605
01606
01607
01608
01609 integer (kind=int_kind), intent(in) ::
01610 nx_block, ny_block
01611
01612 real (kind=dbl_kind), dimension (nx_block,ny_block), intent(in) ::
01613 stressp_1 ,
01614 stressm_1 ,
01615 stress12_1,
01616 prs_sig
01617
01618 real (kind=dbl_kind), dimension (nx_block,ny_block), intent(out)::
01619 sig1 ,
01620 sig2
01621
01622
01623
01624 integer (kind=int_kind) :: i, j
01625
01626 do j = 1, ny_block
01627 do i = 1, nx_block
01628 if (prs_sig(i,j) > puny) then
01629 sig1(i,j) = (p5*(stressp_1(i,j) &
01630 + sqrt(stressm_1(i,j)**2+c4*stress12_1(i,j)**2))) &
01631 / prs_sig(i,j)
01632 sig2(i,j) = (p5*(stressp_1(i,j) &
01633 - sqrt(stressm_1(i,j)**2+c4*stress12_1(i,j)**2))) &
01634 / prs_sig(i,j)
01635 else
01636 sig1(i,j) = spval_dbl
01637 sig2(i,j) = spval_dbl
01638 endif
01639 enddo
01640 enddo
01641
01642 end subroutine principal_stress
01643
01644
01645
01646 end module ice_dyn_evp
01647
01648