- Timestamp:
- 2017-09-13T18:46:56+02:00 (7 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/iceupdate.F90
r8514 r8518 36 36 !!gm end 37 37 USE sbccpl ! Surface boundary condition: coupled interface 38 USE icealb ! albedo parameters38 USE icealb ! sea-ice: albedo parameters 39 39 USE traqsr ! add penetration of solar flux in the calculation of heat budget 40 USE domvvl ! Variable volume 41 USE icectl ! ??? 40 USE icectl ! sea-ice: control prints 42 41 USE bdy_oce , ONLY : ln_bdy 43 42 ! … … 58 57 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: utau_oce, vtau_oce ! air-ocean surface i- & j-stress [N/m2] 59 58 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: tmod_io ! modulus of the ice-ocean velocity [m/s] 60 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: soce_0 , sice_0 ! cst SSS and ice salinity (levitating sea-ice)61 59 62 60 !! * Substitutions … … 73 71 !! *** ROUTINE ice_update_alloc *** 74 72 !!------------------------------------------------------------------- 75 ALLOCATE( soce_0(jpi,jpj) , utau_oce(jpi,jpj) , & 76 & sice_0(jpi,jpj) , vtau_oce(jpi,jpj) , tmod_io(jpi,jpj), STAT=ice_update_alloc) 77 ! 73 ALLOCATE( utau_oce(jpi,jpj), vtau_oce(jpi,jpj), tmod_io(jpi,jpj), STAT=ice_update_alloc ) 74 ! 78 75 IF( lk_mpp ) CALL mpp_sum( ice_update_alloc ) 79 76 IF( ice_update_alloc /= 0 ) CALL ctl_warn('ice_update_alloc: failed to allocate arrays') 77 ! 80 78 END FUNCTION ice_update_alloc 81 79 … … 206 204 ! 207 205 alb_ice(:,:,:) = ( 1._wp - cldf_ice ) * zalb_cs(:,:,:) + cldf_ice * zalb_os(:,:,:) 208 209 ! ! conservation test 210 IF( ln_icediachk .AND. .NOT. ln_bdy) CALL ice_cons_final( 'iceupdate' ) 211 ! ! control prints 212 IF( ln_icectl ) CALL ice_prt( kt, iiceprt, jiceprt, 3, ' - Final state ice_update - ' ) 213 IF( ln_ctl ) CALL ice_prt3D( 'iceupdate' ) 214 ! 215 IF( nn_timing == 1 ) CALL timing_stop('ice_update_flx') 206 ! 207 IF( lrst_ice ) THEN !* write snwice_mass fields in the restart file 208 CALL update_rst( 'WRITE', kt ) 209 ENDIF 210 ! 211 ! controls 212 IF( ln_icediachk .AND. .NOT. ln_bdy) CALL ice_cons_final('iceupdate') ! conservation 213 IF( ln_icectl ) CALL ice_prt (kt, iiceprt, jiceprt, 3, 'Final state ice_update') ! prints 214 IF( ln_ctl ) CALL ice_prt3D ('iceupdate') ! prints 215 IF( nn_timing == 1 ) CALL timing_stop ('ice_update_flx') ! timing 216 216 ! 217 217 END SUBROUTINE ice_update_flx … … 320 320 IF( ice_update_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'ice_update_init : unable to allocate standard arrays' ) 321 321 ! 322 soce_0(:,:) = soce ! constant SSS and ice salinity used in levitating case 0 (i.e. virtual salt flux) 323 sice_0(:,:) = sice 324 WHERE( 14._wp <= glamt(:,:) .AND. glamt(:,:) <= 32._wp .AND. & ! reduced values in the Baltic Sea area 325 & 54._wp <= gphit(:,:) .AND. gphit(:,:) <= 66._wp ) 326 soce_0(:,:) = 4._wp 327 sice_0(:,:) = 2._wp 328 END WHERE 329 ! 330 IF( .NOT.ln_rstart ) THEN ! set 331 ! 332 snwice_mass (:,:) = tmask(:,:,1) * ( rhosn * vt_s(:,:) + rhoic * vt_i(:,:) ) ! snow+ice mass 333 snwice_mass_b(:,:) = snwice_mass(:,:) 334 ! 335 IF( ln_ice_embd ) THEN ! embedded sea-ice: deplete the initial ssh below sea-ice area 336 sshn(:,:) = sshn(:,:) - snwice_mass(:,:) * r1_rau0 337 sshb(:,:) = sshb(:,:) - snwice_mass(:,:) * r1_rau0 338 339 !!gm I really don't like this stuff here... Find a way to put that elsewhere or differently 340 !!gm 341 IF( .NOT.ln_linssh ) THEN 342 DO jk = 1,jpkm1 ! adjust initial vertical scale factors 343 e3t_n(:,:,jk) = e3t_0(:,:,jk) * ( 1._wp + sshn(:,:)*tmask(:,:,1) / (ht_0(:,:) + 1._wp - tmask(:,:,1) ) ) 344 e3t_b(:,:,jk) = e3t_0(:,:,jk) * ( 1._wp + sshb(:,:)*tmask(:,:,1) / (ht_0(:,:) + 1._wp - tmask(:,:,1) ) ) 345 END DO 346 e3t_a(:,:,:) = e3t_b(:,:,:) 347 !!gm we are in no-restart case, so sshn=sshb ==>> faster calculation: 348 !! REAL(wp) :: ze3t ! local scalar 349 !! REAL(wp), DIMENSION(jpi,jpj) :: z2d ! workspace 350 !! 351 !! WHERE( ht_0(:,:) > 0 ) ; z2d(:,:) = 1._wp + sshn(:,:)*tmask(:,:,1) / ht_0(:,:) 352 !! ELSEWHERE ; z2d(:,:) = 1._wp 353 !! END WHERE 354 !! DO jk = 1,jpkm1 ! adjust initial vertical scale factors 355 !! ze3t = e3t_0(:,:,jk) * z2d(:,:) 356 !! e3t_n(:,:,jk) = ze3t 357 !! e3t_b(:,:,jk) = ze3t 358 !! e3t_a(:,:,jk) = ze3t 359 !! END DO 360 !!gm but since it is only done at the initialisation.... just the following can be acceptable: 361 ! DO jk = 1,jpkm1 ! adjust initial vertical scale factors 362 ! e3t_n(:,:,jk) = e3t_0(:,:,jk) * ( 1._wp + sshn(:,:)*tmask(:,:,1) / (ht_0(:,:) + 1._wp - tmask(:,:,1)) ) 363 ! END DO 364 ! e3t_b(:,:,:) = e3t_n(:,:,:) 365 ! e3t_a(:,:,:) = e3t_n(:,:,:) 366 !!gm end 367 ! Reconstruction of all vertical scale factors at now and before time-steps 368 ! ========================================================================= 369 ! Horizontal scale factor interpolations 370 ! -------------------------------------- 371 CALL dom_vvl_interpol( e3t_b(:,:,:), e3u_b(:,:,:), 'U' ) 372 CALL dom_vvl_interpol( e3t_b(:,:,:), e3v_b(:,:,:), 'V' ) 373 CALL dom_vvl_interpol( e3t_n(:,:,:), e3u_n(:,:,:), 'U' ) 374 CALL dom_vvl_interpol( e3t_n(:,:,:), e3v_n(:,:,:), 'V' ) 375 CALL dom_vvl_interpol( e3u_n(:,:,:), e3f_n(:,:,:), 'F' ) 376 ! Vertical scale factor interpolations 377 ! ------------------------------------ 378 CALL dom_vvl_interpol( e3t_n(:,:,:), e3w_n (:,:,:), 'W' ) 379 CALL dom_vvl_interpol( e3u_n(:,:,:), e3uw_n(:,:,:), 'UW' ) 380 CALL dom_vvl_interpol( e3v_n(:,:,:), e3vw_n(:,:,:), 'VW' ) 381 CALL dom_vvl_interpol( e3u_b(:,:,:), e3uw_b(:,:,:), 'UW' ) 382 CALL dom_vvl_interpol( e3v_b(:,:,:), e3vw_b(:,:,:), 'VW' ) 383 ! t- and w- points depth 384 ! ---------------------- 385 !!gm not sure of that.... 386 gdept_n(:,:,1) = 0.5_wp * e3w_n(:,:,1) 387 gdepw_n(:,:,1) = 0.0_wp 388 gde3w_n(:,:,1) = gdept_n(:,:,1) - sshn(:,:) 389 DO jk = 2, jpk 390 gdept_n(:,:,jk) = gdept_n(:,:,jk-1) + e3w_n(:,:,jk ) 391 gdepw_n(:,:,jk) = gdepw_n(:,:,jk-1) + e3t_n(:,:,jk-1) 392 gde3w_n(:,:,jk) = gdept_n(:,:,jk ) - sshn (:,:) 393 END DO 322 CALL update_rst( 'READ' ) !* read or initialize all required files 323 ! 324 END SUBROUTINE ice_update_init 325 326 SUBROUTINE update_rst( cdrw, kt ) 327 !!--------------------------------------------------------------------- 328 !! *** ROUTINE rhg_evp_rst *** 329 !! 330 !! ** Purpose : Read or write RHG file in restart file 331 !! 332 !! ** Method : use of IOM library 333 !!---------------------------------------------------------------------- 334 CHARACTER(len=*) , INTENT(in) :: cdrw ! "READ"/"WRITE" flag 335 INTEGER, OPTIONAL, INTENT(in) :: kt ! ice time-step 336 ! 337 INTEGER :: iter ! local integer 338 INTEGER :: id1 ! local integer 339 !!---------------------------------------------------------------------- 340 ! 341 IF( TRIM(cdrw) == 'READ' ) THEN ! Read/initialize 342 ! ! --------------- 343 IF( ln_rstart ) THEN !* Read the restart file 344 ! 345 id1 = iom_varid( numrir, 'snwice_mass' , ldstop = .FALSE. ) 346 ! 347 IF( id1 > 0 ) THEN ! fields exist 348 CALL iom_get( numrir, jpdom_autoglo, 'snwice_mass' , snwice_mass ) 349 CALL iom_get( numrir, jpdom_autoglo, 'snwice_mass_b', snwice_mass_b ) 350 ELSE ! start from rest 351 IF(lwp) WRITE(numout,*) ' ==>> previous run without snow-ice mass output then set it' 352 snwice_mass (:,:) = tmask(:,:,1) * ( rhosn * vt_s(:,:) + rhoic * vt_i(:,:) ) 353 snwice_mass_b(:,:) = snwice_mass(:,:) 394 354 ENDIF 355 ELSE !* Start from rest 356 IF(lwp) WRITE(numout,*) ' ==>> start from rest: set the snow-ice mass' 357 snwice_mass (:,:) = tmask(:,:,1) * ( rhosn * vt_s(:,:) + rhoic * vt_i(:,:) ) 358 snwice_mass_b(:,:) = snwice_mass(:,:) 395 359 ENDIF 396 ENDIF ! .NOT. ln_rstart 397 ! 398 END SUBROUTINE ice_update_init 360 ! 361 ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN ! Create restart file 362 ! ! ------------------- 363 IF(lwp) WRITE(numout,*) '---- update-rst ----' 364 iter = kt + nn_fsbc - 1 ! ice restarts are written at kt == nitrst - nn_fsbc + 1 365 ! 366 CALL iom_rstput( iter, nitrst, numriw, 'snwice_mass' , snwice_mass ) 367 CALL iom_rstput( iter, nitrst, numriw, 'snwice_mass_b', snwice_mass_b ) 368 ! 369 ENDIF 370 ! 371 END SUBROUTINE update_rst 399 372 400 373 #else
Note: See TracChangeset
for help on using the changeset viewer.
