Changeset 8518
- Timestamp:
- 2017-09-13T18:46:56+02:00 (7 years ago)
- Location:
- branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO
- Files:
-
- 19 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/ice.F90
r8517 r8518 170 170 171 171 ! !!** ice-dynamics namelist (namice_dyn) ** 172 LOGICAL, PUBLIC :: ln_dynFULL !: dyn.: full ice dynamics (rheology + advection + ridging/rafting + correction)173 LOGICAL, PUBLIC :: ln_dynRHGADV !: dyn.: no ridge/raft & no corrections (rheology + advection)174 LOGICAL, PUBLIC :: ln_dynADV !: dyn.: only advection w prescribed vel.(rn_uvice + advection)175 172 REAL(wp), PUBLIC :: rn_ishlat !: lateral boundary condition for sea-ice 176 173 REAL(wp), PUBLIC :: rn_cio !: drag coefficient for oceanic stress … … 180 177 REAL(wp), PUBLIC :: rn_lfrelax !: relaxation time scale (s-1) to reach static friction (landfast ice) 181 178 ! 182 ! !!** ice-rdige/raft namelist (namice_rdgrft) **183 LOGICAL , PUBLIC :: ln_str_H79 !: ice strength parameterization (Hibler79)184 REAL(wp), PUBLIC :: rn_pstar !: determines ice strength, Hibler JPO79185 REAL(wp), PUBLIC :: rn_crhg !: determines changes in ice strength186 LOGICAL , PUBLIC :: ln_str_R75 !: ice strength parameterization (Rothrock75)187 REAL(wp), PUBLIC :: rn_perdg !: ridging work divided by pot. energy change in ridging188 !189 179 ! !!** ice-rheology namelist (namice_rhg) ** 190 LOGICAL , PUBLIC :: ln_rhg_EVP !: EVP rheology191 180 REAL(wp), PUBLIC :: rn_creepl !: creep limit : has to be under 1.0e-9 192 181 REAL(wp), PUBLIC :: rn_ecc !: eccentricity of the elliptical yield curve 193 182 INTEGER , PUBLIC :: nn_nevp !: number of iterations for subcycling 194 183 REAL(wp), PUBLIC :: rn_relast !: ratio => telast/rdt_ice (1/3 or 1/9 depending on nb of subcycling nevp) 195 !196 ! !!** ice-advection namelist (namice_adv) **197 LOGICAL , PUBLIC :: ln_adv_Pra !: Prather advection scheme198 LOGICAL , PUBLIC :: ln_adv_UMx !: Ultimate-Macho advection scheme199 INTEGER , PUBLIC :: nn_UMx !: order of the UMx advection scheme200 184 ! 201 185 ! !!** ice-thermodynamics namelist (namice_thd) ** … … 230 214 231 215 ! !!** ice-salinity namelist (namice_sal) ** 232 LOGICAL , PUBLIC :: ln_icedS !: activate gravity drainage and flushing (T) or not (F)233 216 INTEGER , PUBLIC :: nn_icesal !: salinity configuration used in the model 234 217 ! ! 1 - constant salinity in both space and time … … 236 219 ! ! 3 - salinity profile, constant in time 237 220 REAL(wp), PUBLIC :: rn_icesal !: bulk salinity (ppt) in case of constant salinity 238 REAL(wp), PUBLIC :: rn_sal_gd !: restoring salinity for gravity drainage [PSU]239 REAL(wp), PUBLIC :: rn_time_gd !: restoring time constant for gravity drainage (= 20 days) [s]240 REAL(wp), PUBLIC :: rn_sal_fl !: restoring salinity for flushing [PSU]241 REAL(wp), PUBLIC :: rn_time_fl !: restoring time constant for gravity drainage (= 10 days) [s]242 221 REAL(wp), PUBLIC :: rn_simax !: maximum ice salinity [PSU] 243 222 REAL(wp), PUBLIC :: rn_simin !: minimum ice salinity [PSU] … … 329 308 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hfx_dif !: total heat flux causing Temp change in the ice [W.m-2] 330 309 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hfx_snw !: heat flux for snow melt [W.m-2] 331 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hfx_err !: heat flux error after heat diffusion [W.m-2]332 310 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hfx_err_dif !: heat flux remaining due to change in non-solar flux [W.m-2] 333 311 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hfx_err_rem !: heat flux error after heat remapping [W.m-2] … … 494 472 & sfx_res(jpi,jpj) , sfx_bri(jpi,jpj) , sfx_dyn(jpi,jpj) , sfx_sub(jpi,jpj) , sfx_lam(jpi,jpj) , & 495 473 & sfx_bog(jpi,jpj) , sfx_bom(jpi,jpj) , sfx_sum(jpi,jpj) , sfx_sni(jpi,jpj) , sfx_opw(jpi,jpj) , & 496 & hfx_res(jpi,jpj) , hfx_snw(jpi,jpj) , hfx_sub(jpi,jpj) , hfx_err(jpi,jpj) ,&474 & hfx_res(jpi,jpj) , hfx_snw(jpi,jpj) , hfx_sub(jpi,jpj) , & 497 475 & hfx_in (jpi,jpj) , hfx_out(jpi,jpj) , fhld (jpi,jpj) , & 498 476 & hfx_sum(jpi,jpj) , hfx_bom(jpi,jpj) , hfx_bog(jpi,jpj) , hfx_dif(jpi,jpj) , & -
branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/ice1D.F90
r8486 r8518 47 47 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: hfx_opw_1d 48 48 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: hfx_snw_1d 49 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: hfx_err_1d50 49 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: hfx_err_rem_1d 51 50 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: hfx_err_dif_1d … … 178 177 & rn_amax_1d(jpij) , & 179 178 & hfx_thd_1d(jpij) , hfx_spr_1d(jpij) , & 180 & hfx_snw_1d(jpij) , hfx_sub_1d(jpij) , hfx_err_1d(jpij) ,&179 & hfx_snw_1d(jpij) , hfx_sub_1d(jpij) , & 181 180 & hfx_res_1d(jpij) , hfx_err_rem_1d(jpij) , hfx_err_dif_1d(jpij) , hfx_out_1d(jpij), STAT=ierr(ii) ) 182 181 ! -
branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/iceadv.F90
r8517 r8518 41 41 INTEGER, PARAMETER :: np_advPRA = 1 ! Prather scheme 42 42 INTEGER, PARAMETER :: np_advUMx = 2 ! Ultimate-Macho scheme 43 43 ! 44 ! ** namelist (namice_adv) ** 45 LOGICAL :: ln_adv_Pra ! Prather advection scheme 46 LOGICAL :: ln_adv_UMx ! Ultimate-Macho advection scheme 47 INTEGER :: nn_UMx ! order of the UMx advection scheme 48 ! 44 49 !! * Substitution 45 50 # include "vectopt_loop_substitute.h90" … … 85 90 CASE( np_advUMx ) ! ULTIMATE-MACHO scheme ! 86 91 ! !-----------------------! 87 CALL ice_adv_umx( kt, u_ice, v_ice, &92 CALL ice_adv_umx( nn_UMx, kt, u_ice, v_ice, & 88 93 & ato_i, v_i, v_s, smv_i, oa_i, a_i, a_ip, v_ip, e_s, e_i ) 89 94 ! !-----------------------! … … 107 112 IF( iom_use('deitrp') ) CALL iom_put( "deitrp" , diag_trp_ei ) ! advected ice enthalpy (W/m2) 108 113 IF( iom_use('destrp') ) CALL iom_put( "destrp" , diag_trp_es ) ! advected snw enthalpy (W/m2) 114 115 IF( lrst_ice ) THEN !* write Prather fields in the restart file 116 IF( ln_adv_Pra ) CALL adv_pra_rst( 'WRITE', kt ) 117 ENDIF 109 118 110 119 ! controls … … 158 167 IF( ioptio /= 1 ) CALL ctl_stop( 'ice_adv_init: choose one and only one ice advection scheme (ln_adv_Pra or ln_adv_UMx)' ) 159 168 ! 169 IF( ln_adv_Pra ) CALL adv_pra_rst( 'READ' ) !* read or initialize all required files 170 ! 160 171 END SUBROUTINE ice_adv_init 161 172 -
branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/iceadv_prather.F90
r8512 r8518 20 20 USE lbclnk ! lateral boundary condition - MPP exchanges 21 21 USE in_out_manager ! I/O manager 22 USE iom ! I/O library 22 23 USE prtctl ! Print control 23 24 USE lib_mpp ! MPP library … … 28 29 29 30 PUBLIC ice_adv_prather ! called by iceadv 31 PUBLIC adv_pra_rst ! called by iceadv 30 32 31 33 !! * Substitutions … … 592 594 END SUBROUTINE adv_y 593 595 596 SUBROUTINE adv_pra_rst( cdrw, kt ) 597 !!--------------------------------------------------------------------- 598 !! *** ROUTINE adv_pra_rst *** 599 !! 600 !! ** Purpose : Read or write RHG file in restart file 601 !! 602 !! ** Method : use of IOM library 603 !!---------------------------------------------------------------------- 604 CHARACTER(len=*) , INTENT(in) :: cdrw ! "READ"/"WRITE" flag 605 INTEGER, OPTIONAL, INTENT(in) :: kt ! ice time-step 606 ! 607 INTEGER :: jk, jl ! dummy loop indices 608 INTEGER :: iter ! local integer 609 INTEGER :: id1 ! local integer 610 CHARACTER(len=25) :: znam 611 CHARACTER(len=2) :: zchar, zchar1 612 REAL(wp), DIMENSION(jpi,jpj) :: z2d 613 !!---------------------------------------------------------------------- 614 ! 615 IF( TRIM(cdrw) == 'READ' ) THEN ! Read/initialize 616 ! ! --------------- 617 IF( ln_rstart ) THEN !* Read the restart file 618 ! 619 id1 = iom_varid( numrir, 'sxopw' , ldstop = .FALSE. ) 620 ! 621 IF( id1 > 0 ) THEN ! fields exist 622 DO jl = 1, jpl 623 WRITE(zchar,'(I2.2)') jl 624 znam = 'sxice'//'_htc'//zchar 625 CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 626 sxice(:,:,jl) = z2d(:,:) 627 znam = 'syice'//'_htc'//zchar 628 CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 629 syice(:,:,jl) = z2d(:,:) 630 znam = 'sxxice'//'_htc'//zchar 631 CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 632 sxxice(:,:,jl) = z2d(:,:) 633 znam = 'syyice'//'_htc'//zchar 634 CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 635 syyice(:,:,jl) = z2d(:,:) 636 znam = 'sxyice'//'_htc'//zchar 637 CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 638 sxyice(:,:,jl) = z2d(:,:) 639 znam = 'sxsn'//'_htc'//zchar 640 CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 641 sxsn(:,:,jl) = z2d(:,:) 642 znam = 'sysn'//'_htc'//zchar 643 CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 644 sysn(:,:,jl) = z2d(:,:) 645 znam = 'sxxsn'//'_htc'//zchar 646 CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 647 sxxsn(:,:,jl) = z2d(:,:) 648 znam = 'syysn'//'_htc'//zchar 649 CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 650 syysn(:,:,jl) = z2d(:,:) 651 znam = 'sxysn'//'_htc'//zchar 652 CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 653 sxysn(:,:,jl) = z2d(:,:) 654 znam = 'sxa'//'_htc'//zchar 655 CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 656 sxa(:,:,jl) = z2d(:,:) 657 znam = 'sya'//'_htc'//zchar 658 CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 659 sya(:,:,jl) = z2d(:,:) 660 znam = 'sxxa'//'_htc'//zchar 661 CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 662 sxxa(:,:,jl) = z2d(:,:) 663 znam = 'syya'//'_htc'//zchar 664 CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 665 syya(:,:,jl) = z2d(:,:) 666 znam = 'sxya'//'_htc'//zchar 667 CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 668 sxya(:,:,jl) = z2d(:,:) 669 znam = 'sxc0'//'_htc'//zchar 670 CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 671 sxc0(:,:,jl) = z2d(:,:) 672 znam = 'syc0'//'_htc'//zchar 673 CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 674 syc0(:,:,jl) = z2d(:,:) 675 znam = 'sxxc0'//'_htc'//zchar 676 CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 677 sxxc0(:,:,jl) = z2d(:,:) 678 znam = 'syyc0'//'_htc'//zchar 679 CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 680 syyc0(:,:,jl) = z2d(:,:) 681 znam = 'sxyc0'//'_htc'//zchar 682 CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 683 sxyc0(:,:,jl) = z2d(:,:) 684 znam = 'sxsal'//'_htc'//zchar 685 CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 686 sxsal(:,:,jl) = z2d(:,:) 687 znam = 'sysal'//'_htc'//zchar 688 CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 689 sysal(:,:,jl) = z2d(:,:) 690 znam = 'sxxsal'//'_htc'//zchar 691 CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 692 sxxsal(:,:,jl) = z2d(:,:) 693 znam = 'syysal'//'_htc'//zchar 694 CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 695 syysal(:,:,jl) = z2d(:,:) 696 znam = 'sxysal'//'_htc'//zchar 697 CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 698 sxysal(:,:,jl) = z2d(:,:) 699 znam = 'sxage'//'_htc'//zchar 700 CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 701 sxage(:,:,jl) = z2d(:,:) 702 znam = 'syage'//'_htc'//zchar 703 CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 704 syage(:,:,jl) = z2d(:,:) 705 znam = 'sxxage'//'_htc'//zchar 706 CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 707 sxxage(:,:,jl) = z2d(:,:) 708 znam = 'syyage'//'_htc'//zchar 709 CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 710 syyage(:,:,jl) = z2d(:,:) 711 znam = 'sxyage'//'_htc'//zchar 712 CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 713 sxyage(:,:,jl)= z2d(:,:) 714 END DO 715 ! MV MP 2016 716 IF ( nn_pnd_scheme > 0 ) THEN 717 DO jl = 1, jpl 718 WRITE(zchar,'(I2.2)') jl 719 znam = 'sxap'//'_htc'//zchar 720 CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 721 sxap(:,:,jl) = z2d(:,:) 722 znam = 'syap'//'_htc'//zchar 723 CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 724 syap(:,:,jl) = z2d(:,:) 725 znam = 'sxxap'//'_htc'//zchar 726 CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 727 sxxap(:,:,jl) = z2d(:,:) 728 znam = 'syyap'//'_htc'//zchar 729 CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 730 syyap(:,:,jl) = z2d(:,:) 731 znam = 'sxyap'//'_htc'//zchar 732 CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 733 sxyap(:,:,jl) = z2d(:,:) 734 735 znam = 'sxvp'//'_htc'//zchar 736 CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 737 sxvp(:,:,jl) = z2d(:,:) 738 znam = 'syvp'//'_htc'//zchar 739 CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 740 syvp(:,:,jl) = z2d(:,:) 741 znam = 'sxxvp'//'_htc'//zchar 742 CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 743 sxxvp(:,:,jl) = z2d(:,:) 744 znam = 'syyvp'//'_htc'//zchar 745 CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 746 syyvp(:,:,jl) = z2d(:,:) 747 znam = 'sxyvp'//'_htc'//zchar 748 CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 749 sxyvp(:,:,jl) = z2d(:,:) 750 END DO 751 ENDIF 752 ! END MV MP 2016 753 754 CALL iom_get( numrir, jpdom_autoglo, 'sxopw ' , sxopw ) 755 CALL iom_get( numrir, jpdom_autoglo, 'syopw ' , syopw ) 756 CALL iom_get( numrir, jpdom_autoglo, 'sxxopw' , sxxopw ) 757 CALL iom_get( numrir, jpdom_autoglo, 'syyopw' , syyopw ) 758 CALL iom_get( numrir, jpdom_autoglo, 'sxyopw' , sxyopw ) 759 760 DO jl = 1, jpl 761 WRITE(zchar,'(I2.2)') jl 762 DO jk = 1, nlay_i 763 WRITE(zchar1,'(I2.2)') jk 764 znam = 'sxe'//'_il'//zchar1//'_htc'//zchar 765 CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 766 sxe(:,:,jk,jl) = z2d(:,:) 767 znam = 'sye'//'_il'//zchar1//'_htc'//zchar 768 CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 769 sye(:,:,jk,jl) = z2d(:,:) 770 znam = 'sxxe'//'_il'//zchar1//'_htc'//zchar 771 CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 772 sxxe(:,:,jk,jl) = z2d(:,:) 773 znam = 'syye'//'_il'//zchar1//'_htc'//zchar 774 CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 775 syye(:,:,jk,jl) = z2d(:,:) 776 znam = 'sxye'//'_il'//zchar1//'_htc'//zchar 777 CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 778 sxye(:,:,jk,jl) = z2d(:,:) 779 END DO 780 END DO 781 ! 782 ELSE ! start rheology from rest 783 IF(lwp) WRITE(numout,*) ' ==>> previous run without Prather, set moments to 0' 784 sxopw (:,:) = 0._wp ; sxice (:,:,:) = 0._wp ; sxsn (:,:,:) = 0._wp ; sxa (:,:,:) = 0._wp 785 syopw (:,:) = 0._wp ; syice (:,:,:) = 0._wp ; sysn (:,:,:) = 0._wp ; sya (:,:,:) = 0._wp 786 sxxopw(:,:) = 0._wp ; sxxice(:,:,:) = 0._wp ; sxxsn(:,:,:) = 0._wp ; sxxa (:,:,:) = 0._wp 787 syyopw(:,:) = 0._wp ; syyice(:,:,:) = 0._wp ; syysn(:,:,:) = 0._wp ; syya (:,:,:) = 0._wp 788 sxyopw(:,:) = 0._wp ; sxyice(:,:,:) = 0._wp ; sxysn(:,:,:) = 0._wp ; sxya (:,:,:) = 0._wp 789 ! 790 sxc0 (:,:,:) = 0._wp ; sxe (:,:,:,:) = 0._wp ; sxsal (:,:,:) = 0._wp ; sxage (:,:,:) = 0._wp 791 syc0 (:,:,:) = 0._wp ; sye (:,:,:,:) = 0._wp ; sysal (:,:,:) = 0._wp ; syage (:,:,:) = 0._wp 792 sxxc0 (:,:,:) = 0._wp ; sxxe (:,:,:,:) = 0._wp ; sxxsal (:,:,:) = 0._wp ; sxxage (:,:,:) = 0._wp 793 syyc0 (:,:,:) = 0._wp ; syye (:,:,:,:) = 0._wp ; syysal (:,:,:) = 0._wp ; syyage (:,:,:) = 0._wp 794 sxyc0 (:,:,:) = 0._wp ; sxye (:,:,:,:) = 0._wp ; sxysal (:,:,:) = 0._wp ; sxyage (:,:,:) = 0._wp 795 ! MV MP 2016 796 IF ( nn_pnd_scheme > 0 ) THEN 797 sxap (:,:,:) = 0._wp ; sxvp (:,:,:) = 0._wp 798 syap (:,:,:) = 0._wp ; syvp (:,:,:) = 0._wp 799 sxxap (:,:,:) = 0._wp ; sxxvp (:,:,:) = 0._wp 800 syyap (:,:,:) = 0._wp ; syyvp (:,:,:) = 0._wp 801 sxyap (:,:,:) = 0._wp ; sxyvp (:,:,:) = 0._wp 802 ENDIF 803 ! END MV MP 2016 804 ENDIF 805 ELSE !* Start from rest 806 IF(lwp) WRITE(numout,*) ' ==>> start from rest: set moments to 0' 807 sxopw (:,:) = 0._wp ; sxice (:,:,:) = 0._wp ; sxsn (:,:,:) = 0._wp ; sxa (:,:,:) = 0._wp 808 syopw (:,:) = 0._wp ; syice (:,:,:) = 0._wp ; sysn (:,:,:) = 0._wp ; sya (:,:,:) = 0._wp 809 sxxopw(:,:) = 0._wp ; sxxice(:,:,:) = 0._wp ; sxxsn(:,:,:) = 0._wp ; sxxa (:,:,:) = 0._wp 810 syyopw(:,:) = 0._wp ; syyice(:,:,:) = 0._wp ; syysn(:,:,:) = 0._wp ; syya (:,:,:) = 0._wp 811 sxyopw(:,:) = 0._wp ; sxyice(:,:,:) = 0._wp ; sxysn(:,:,:) = 0._wp ; sxya (:,:,:) = 0._wp 812 ! 813 sxc0 (:,:,:) = 0._wp ; sxe (:,:,:,:) = 0._wp ; sxsal (:,:,:) = 0._wp ; sxage (:,:,:) = 0._wp 814 syc0 (:,:,:) = 0._wp ; sye (:,:,:,:) = 0._wp ; sysal (:,:,:) = 0._wp ; syage (:,:,:) = 0._wp 815 sxxc0 (:,:,:) = 0._wp ; sxxe (:,:,:,:) = 0._wp ; sxxsal (:,:,:) = 0._wp ; sxxage (:,:,:) = 0._wp 816 syyc0 (:,:,:) = 0._wp ; syye (:,:,:,:) = 0._wp ; syysal (:,:,:) = 0._wp ; syyage (:,:,:) = 0._wp 817 sxyc0 (:,:,:) = 0._wp ; sxye (:,:,:,:) = 0._wp ; sxysal (:,:,:) = 0._wp ; sxyage (:,:,:) = 0._wp 818 ! MV MP 2016 819 IF ( nn_pnd_scheme > 0 ) THEN 820 sxap (:,:,:) = 0._wp ; sxvp (:,:,:) = 0._wp 821 syap (:,:,:) = 0._wp ; syvp (:,:,:) = 0._wp 822 sxxap (:,:,:) = 0._wp ; sxxvp (:,:,:) = 0._wp 823 syyap (:,:,:) = 0._wp ; syyvp (:,:,:) = 0._wp 824 sxyap (:,:,:) = 0._wp ; sxyvp (:,:,:) = 0._wp 825 ENDIF 826 ! END MV MP 2016 827 ENDIF 828 ! 829 ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN ! Create restart file 830 ! ! ------------------- 831 IF(lwp) WRITE(numout,*) '---- adv-rst ----' 832 iter = kt + nn_fsbc - 1 ! ice restarts are written at kt == nitrst - nn_fsbc + 1 833 ! 834 DO jl = 1, jpl 835 WRITE(zchar,'(I2.2)') jl 836 znam = 'sxice'//'_htc'//zchar 837 z2d(:,:) = sxice(:,:,jl) 838 CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 839 znam = 'syice'//'_htc'//zchar 840 z2d(:,:) = syice(:,:,jl) 841 CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 842 znam = 'sxxice'//'_htc'//zchar 843 z2d(:,:) = sxxice(:,:,jl) 844 CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 845 znam = 'syyice'//'_htc'//zchar 846 z2d(:,:) = syyice(:,:,jl) 847 CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 848 znam = 'sxyice'//'_htc'//zchar 849 z2d(:,:) = sxyice(:,:,jl) 850 CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 851 znam = 'sxsn'//'_htc'//zchar 852 z2d(:,:) = sxsn(:,:,jl) 853 CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 854 znam = 'sysn'//'_htc'//zchar 855 z2d(:,:) = sysn(:,:,jl) 856 CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 857 znam = 'sxxsn'//'_htc'//zchar 858 z2d(:,:) = sxxsn(:,:,jl) 859 CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 860 znam = 'syysn'//'_htc'//zchar 861 z2d(:,:) = syysn(:,:,jl) 862 CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 863 znam = 'sxysn'//'_htc'//zchar 864 z2d(:,:) = sxysn(:,:,jl) 865 CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 866 znam = 'sxa'//'_htc'//zchar 867 z2d(:,:) = sxa(:,:,jl) 868 CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 869 znam = 'sya'//'_htc'//zchar 870 z2d(:,:) = sya(:,:,jl) 871 CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 872 znam = 'sxxa'//'_htc'//zchar 873 z2d(:,:) = sxxa(:,:,jl) 874 CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 875 znam = 'syya'//'_htc'//zchar 876 z2d(:,:) = syya(:,:,jl) 877 CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 878 znam = 'sxya'//'_htc'//zchar 879 z2d(:,:) = sxya(:,:,jl) 880 CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 881 znam = 'sxc0'//'_htc'//zchar 882 z2d(:,:) = sxc0(:,:,jl) 883 CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 884 znam = 'syc0'//'_htc'//zchar 885 z2d(:,:) = syc0(:,:,jl) 886 CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 887 znam = 'sxxc0'//'_htc'//zchar 888 z2d(:,:) = sxxc0(:,:,jl) 889 CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 890 znam = 'syyc0'//'_htc'//zchar 891 z2d(:,:) = syyc0(:,:,jl) 892 CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 893 znam = 'sxyc0'//'_htc'//zchar 894 z2d(:,:) = sxyc0(:,:,jl) 895 CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 896 znam = 'sxsal'//'_htc'//zchar 897 z2d(:,:) = sxsal(:,:,jl) 898 CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 899 znam = 'sysal'//'_htc'//zchar 900 z2d(:,:) = sysal(:,:,jl) 901 CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 902 znam = 'sxxsal'//'_htc'//zchar 903 z2d(:,:) = sxxsal(:,:,jl) 904 CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 905 znam = 'syysal'//'_htc'//zchar 906 z2d(:,:) = syysal(:,:,jl) 907 CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 908 znam = 'sxysal'//'_htc'//zchar 909 z2d(:,:) = sxysal(:,:,jl) 910 CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 911 znam = 'sxage'//'_htc'//zchar 912 z2d(:,:) = sxage(:,:,jl) 913 CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 914 znam = 'syage'//'_htc'//zchar 915 z2d(:,:) = syage(:,:,jl) 916 CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 917 znam = 'sxxage'//'_htc'//zchar 918 z2d(:,:) = sxxage(:,:,jl) 919 CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 920 znam = 'syyage'//'_htc'//zchar 921 z2d(:,:) = syyage(:,:,jl) 922 CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 923 znam = 'sxyage'//'_htc'//zchar 924 z2d(:,:) = sxyage(:,:,jl) 925 CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 926 END DO 927 928 CALL iom_rstput( iter, nitrst, numriw, 'sxopw ' , sxopw ) 929 CALL iom_rstput( iter, nitrst, numriw, 'syopw ' , syopw ) 930 CALL iom_rstput( iter, nitrst, numriw, 'sxxopw' , sxxopw ) 931 CALL iom_rstput( iter, nitrst, numriw, 'syyopw' , syyopw ) 932 CALL iom_rstput( iter, nitrst, numriw, 'sxyopw' , sxyopw ) 933 934 DO jl = 1, jpl 935 WRITE(zchar,'(I2.2)') jl 936 DO jk = 1, nlay_i 937 WRITE(zchar1,'(I2.2)') jk 938 znam = 'sxe'//'_il'//zchar1//'_htc'//zchar 939 z2d(:,:) = sxe(:,:,jk,jl) 940 CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 941 znam = 'sye'//'_il'//zchar1//'_htc'//zchar 942 z2d(:,:) = sye(:,:,jk,jl) 943 CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 944 znam = 'sxxe'//'_il'//zchar1//'_htc'//zchar 945 z2d(:,:) = sxxe(:,:,jk,jl) 946 CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 947 znam = 'syye'//'_il'//zchar1//'_htc'//zchar 948 z2d(:,:) = syye(:,:,jk,jl) 949 CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 950 znam = 'sxye'//'_il'//zchar1//'_htc'//zchar 951 z2d(:,:) = sxye(:,:,jk,jl) 952 CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 953 END DO 954 END DO 955 ! MV MP 2016 956 IF ( nn_pnd_scheme > 0 ) THEN 957 DO jl = 1, jpl 958 WRITE(zchar,'(I2.2)') jl 959 znam = 'sxap'//'_htc'//zchar 960 z2d(:,:) = sxap(:,:,jl) 961 CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 962 znam = 'syap'//'_htc'//zchar 963 z2d(:,:) = syap(:,:,jl) 964 CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 965 znam = 'sxxap'//'_htc'//zchar 966 z2d(:,:) = sxxap(:,:,jl) 967 CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 968 znam = 'syyap'//'_htc'//zchar 969 z2d(:,:) = syyap(:,:,jl) 970 CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 971 znam = 'sxyap'//'_htc'//zchar 972 z2d(:,:) = sxyap(:,:,jl) 973 CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 974 975 znam = 'sxvp'//'_htc'//zchar 976 z2d(:,:) = sxvp(:,:,jl) 977 CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 978 znam = 'syvp'//'_htc'//zchar 979 z2d(:,:) = syvp(:,:,jl) 980 CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 981 znam = 'sxxvp'//'_htc'//zchar 982 z2d(:,:) = sxxvp(:,:,jl) 983 CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 984 znam = 'syyvp'//'_htc'//zchar 985 z2d(:,:) = syyvp(:,:,jl) 986 CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 987 znam = 'sxyvp'//'_htc'//zchar 988 z2d(:,:) = sxyvp(:,:,jl) 989 CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 990 END DO 991 ENDIF 992 ! 993 ENDIF 994 ! 995 END SUBROUTINE adv_pra_rst 996 594 997 #else 595 998 !!---------------------------------------------------------------------- -
branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/iceadv_umx.F90
r8512 r8518 29 29 PRIVATE 30 30 31 PUBLIC ice_adv_umx ! routinecalled by iceadv.F9031 PUBLIC ice_adv_umx ! called by iceadv.F90 32 32 33 33 REAL(wp) :: z1_6 = 1._wp / 6._wp ! =1/6 … … 43 43 CONTAINS 44 44 45 SUBROUTINE ice_adv_umx( k t, pu_ice, pv_ice, &45 SUBROUTINE ice_adv_umx( k_order, kt, pu_ice, pv_ice, & 46 46 & pato_i, pv_i, pv_s, psmv_i, poa_i, pa_i, pa_ip, pv_ip, pe_s, pe_i ) 47 47 !!---------------------------------------------------------------------- … … 54 54 !! Reference : Leonard, B.P., 1991, Comput. Methods Appl. Mech. Eng., 88, 17-74. 55 55 !!---------------------------------------------------------------------- 56 INTEGER , INTENT(in ) :: k_order ! order of the scheme (1-5 or 20) 56 57 INTEGER , INTENT(in ) :: kt ! time step 57 58 REAL(wp), DIMENSION(:,:) , INTENT(in ) :: pu_ice ! ice i-velocity … … 112 113 !---------------! 113 114 DO jt = 1, initad 114 CALL adv_umx( k t, zdt, zudy, zvdx, zcu_box, zcv_box, pato_i(:,:) ) ! Open water area115 CALL adv_umx( k_order, kt, zdt, zudy, zvdx, zcu_box, zcv_box, pato_i(:,:) ) ! Open water area 115 116 DO jl = 1, jpl 116 CALL adv_umx( k t, zdt, zudy, zvdx, zcu_box, zcv_box, pa_i(:,:,jl) ) ! Ice area117 CALL adv_umx( k t, zdt, zudy, zvdx, zcu_box, zcv_box, pv_i(:,:,jl) ) ! Ice volume118 CALL adv_umx( k t, zdt, zudy, zvdx, zcu_box, zcv_box, psmv_i(:,:,jl) ) ! Salt content119 CALL adv_umx( k t, zdt, zudy, zvdx, zcu_box, zcv_box, poa_i (:,:,jl) ) ! Age content117 CALL adv_umx( k_order, kt, zdt, zudy, zvdx, zcu_box, zcv_box, pa_i(:,:,jl) ) ! Ice area 118 CALL adv_umx( k_order, kt, zdt, zudy, zvdx, zcu_box, zcv_box, pv_i(:,:,jl) ) ! Ice volume 119 CALL adv_umx( k_order, kt, zdt, zudy, zvdx, zcu_box, zcv_box, psmv_i(:,:,jl) ) ! Salt content 120 CALL adv_umx( k_order, kt, zdt, zudy, zvdx, zcu_box, zcv_box, poa_i (:,:,jl) ) ! Age content 120 121 DO jk = 1, nlay_i 121 CALL adv_umx( k t, zdt, zudy, zvdx, zcu_box, zcv_box, pe_i(:,:,jk,jl) ) ! Ice heat content122 END DO 123 CALL adv_umx( k t, zdt, zudy, zvdx, zcu_box, zcv_box, pv_s(:,:,jl) ) ! Snow volume124 CALL adv_umx( k t, zdt, zudy, zvdx, zcu_box, zcv_box, pe_s(:,:,1,jl) ) ! Snow heat content122 CALL adv_umx( k_order, kt, zdt, zudy, zvdx, zcu_box, zcv_box, pe_i(:,:,jk,jl) ) ! Ice heat content 123 END DO 124 CALL adv_umx( k_order, kt, zdt, zudy, zvdx, zcu_box, zcv_box, pv_s(:,:,jl) ) ! Snow volume 125 CALL adv_umx( k_order, kt, zdt, zudy, zvdx, zcu_box, zcv_box, pe_s(:,:,1,jl) ) ! Snow heat content 125 126 IF ( nn_pnd_scheme > 0 ) THEN 126 CALL adv_umx( k t, zdt, zudy, zvdx, zcu_box, zcv_box, pa_ip(:,:,jl) ) ! Melt pond fraction127 CALL adv_umx( k t, zdt, zudy, zvdx, zcu_box, zcv_box, pv_ip(:,:,jl) ) ! Melt pond volume127 CALL adv_umx( k_order, kt, zdt, zudy, zvdx, zcu_box, zcv_box, pa_ip(:,:,jl) ) ! Melt pond fraction 128 CALL adv_umx( k_order, kt, zdt, zudy, zvdx, zcu_box, zcv_box, pv_ip(:,:,jl) ) ! Melt pond volume 128 129 ENDIF 129 130 END DO … … 134 135 END SUBROUTINE ice_adv_umx 135 136 136 SUBROUTINE adv_umx( k t, pdt, puc, pvc, pubox, pvbox, ptc )137 SUBROUTINE adv_umx( k_order, kt, pdt, puc, pvc, pubox, pvbox, ptc ) 137 138 !!---------------------------------------------------------------------- 138 139 !! *** ROUTINE adv_umx *** … … 147 148 !! ** Action : - pt the after advective tracer 148 149 !!---------------------------------------------------------------------- 150 INTEGER , INTENT(in ) :: k_order ! order of the ULTIMATE scheme 149 151 INTEGER , INTENT(in ) :: kt ! number of iteration 150 152 REAL(wp) , INTENT(in ) :: pdt ! tracer time-step … … 189 191 ! High order (_ho) fluxes 190 192 ! ----------------------- 191 SELECT CASE( nn_UMx)193 SELECT CASE( k_order ) 192 194 CASE ( 20 ) ! centered second order 193 195 DO jj = 2, jpjm1 … … 199 201 ! 200 202 CASE ( 1:5 ) ! 1st to 5th order ULTIMATE-MACHO scheme 201 CALL macho( k t, nn_UMx, pdt, ptc, puc, pvc, pubox, pvbox, zt_u, zt_v )203 CALL macho( k_order, kt, pdt, ptc, puc, pvc, pubox, pvbox, zt_u, zt_v ) 202 204 ! 203 205 DO jj = 2, jpjm1 … … 240 242 241 243 242 SUBROUTINE macho( k t, k_order, pdt, ptc, puc, pvc, pubox, pvbox, pt_u, pt_v )244 SUBROUTINE macho( k_order, kt, pdt, ptc, puc, pvc, pubox, pvbox, pt_u, pt_v ) 243 245 !!--------------------------------------------------------------------- 244 246 !! *** ROUTINE ultimate_x *** … … 251 253 !! Reference : Leonard, B.P., 1991, Comput. Methods Appl. Mech. Eng., 88, 17-74. 252 254 !!---------------------------------------------------------------------- 255 INTEGER , INTENT(in ) :: k_order ! order of the ULTIMATE scheme 253 256 INTEGER , INTENT(in ) :: kt ! number of iteration 254 INTEGER , INTENT(in ) :: k_order ! order of the ULTIMATE scheme255 257 REAL(wp) , INTENT(in ) :: pdt ! tracer time-step 256 258 REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: ptc ! tracer fields -
branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icedia.F90
r8514 r8518 45 45 !!---------------------------------------------------------------------- 46 46 CONTAINS 47 48 INTEGER FUNCTION ice_dia_alloc() 49 !!---------------------------------------------------------------------! 50 !! *** ROUTINE ice_rdgrft_alloc *** 51 !!---------------------------------------------------------------------! 52 ALLOCATE( vol_loc_ini(jpi,jpj), sal_loc_ini(jpi,jpj), tem_loc_ini(jpi,jpj), STAT=ice_dia_alloc ) 53 54 IF( lk_mpp ) CALL mpp_sum ( ice_dia_alloc ) 55 IF( ice_dia_alloc /= 0 ) CALL ctl_warn( 'ice_dia_alloc: failed to allocate arrays' ) 56 ! 57 END FUNCTION ice_dia_alloc 47 58 48 59 SUBROUTINE ice_dia( kt ) … … 188 199 ! 189 200 IF( ln_icediahsb ) THEN 190 ALLOCATE( vol_loc_ini(jpi,jpj), sal_loc_ini(jpi,jpj), tem_loc_ini(jpi,jpj), STAT=ierror ) 191 IF( ierror > 0 ) THEN 192 CALL ctl_stop( 'ice_dia: unable to allocate vol_loc_ini' ) 193 RETURN 194 ENDIF 195 ! 196 CALL ice_dia_rst( 'READ' ) !* read or initialize all required files 201 IF( ice_dia_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'ice_dia_init : unable to allocate arrays' ) ! allocate tke arrays 202 CALL ice_dia_rst( 'READ' ) ! read or initialize all required files 197 203 ENDIF 198 204 ! … … 220 226 CALL iom_get( numrir, 'kt_ice' , ziter ) 221 227 IF(lwp) WRITE(numout,*) 222 IF(lwp) WRITE(numout,*) ' 223 IF(lwp) WRITE(numout,*) '~~~~~~~ '228 IF(lwp) WRITE(numout,*) 'ice_dia_rst read at time step = ', ziter 229 IF(lwp) WRITE(numout,*) '~~~~~~~~~~' 224 230 CALL iom_get( numrir, 'frc_voltop' , frc_voltop ) 225 231 CALL iom_get( numrir, 'frc_volbot' , frc_volbot ) … … 252 258 IF( iter == nitrst ) THEN 253 259 IF(lwp) WRITE(numout,*) 254 IF(lwp) WRITE(numout,*) ' 255 IF(lwp) WRITE(numout,*) '~~~~~~~ '260 IF(lwp) WRITE(numout,*) 'ice_dia_rst write at time step = ', kt 261 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~' 256 262 ENDIF 257 263 ! -
branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icedyn.F90
r8517 r8518 42 42 ! 43 43 ! ** namelist (namdyn) ** 44 REAL(wp) :: rn_uice ! prescribed u-vel (case np_dynADV) 45 REAL(wp) :: rn_vice ! prescribed v-vel (case np_dynADV) 44 LOGICAL :: ln_dynFULL ! full ice dynamics (rheology + advection + ridging/rafting + correction) 45 LOGICAL :: ln_dynRHGADV ! no ridge/raft & no corrections (rheology + advection) 46 LOGICAL :: ln_dynADV ! only advection w prescribed vel.(rn_uvice + advection) 47 REAL(wp) :: rn_uice ! prescribed u-vel (case np_dynADV) 48 REAL(wp) :: rn_vice ! prescribed v-vel (case np_dynADV) 46 49 47 50 !! * Substitutions -
branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/iceistate.F90
r8515 r8518 17 17 !! ice_istate_init : initialization of ice state and namelist read 18 18 !!---------------------------------------------------------------------- 19 USE par_oce ! ocean parameters 19 20 USE phycst ! physical constant 20 21 USE oce ! dynamics and tracers variables 21 22 USE dom_oce ! ocean domain 22 USE sbc_oce , ONLY : sst_m, sss_m 23 USE sbc_ice , ONLY : tn_ice 23 USE sbc_oce , ONLY : sst_m, sss_m, ln_ice_embd 24 USE sbc_ice , ONLY : tn_ice, snwice_mass, snwice_mass_b 24 25 USE eosbn2 ! equation of state 26 USE domvvl ! Variable volume 25 27 USE ice ! sea-ice variables 26 USE par_oce ! ocean parameters27 28 USE icevar ! ice_var_salprof 28 29 ! … … 91 92 !! where there is no ice (clem: I do not know why, is it mandatory?) 92 93 !!-------------------------------------------------------------------- 93 INTEGER :: ji, jj, jk, jl ! dummy loop indices 94 REAL(wp) :: ztmelts, zdh 95 INTEGER :: i_hemis, i_fill, jl0 96 REAL(wp) :: zarg, zV, zconv, zdv 97 REAL(wp), DIMENSION(jpi,jpj) :: zswitch ! ice indicator 98 REAL(wp), DIMENSION(jpi,jpj) :: zht_i_ini, zat_i_ini, zvt_i_ini !data from namelist or nc file 99 REAL(wp), DIMENSION(jpi,jpj) :: zts_u_ini, zht_s_ini, zsm_i_ini, ztm_i_ini !data from namelist or nc file 100 REAL(wp), DIMENSION(jpi,jpj,jpl) :: zh_i_ini, za_i_ini !data by cattegories to fill 101 INTEGER , DIMENSION(4) :: itest 94 INTEGER :: ji, jj, jk, jl ! dummy loop indices 95 REAL(wp) :: ztmelts, zdh 96 INTEGER :: i_hemis, i_fill, jl0 97 REAL(wp) :: zarg, zV, zconv, zdv 98 INTEGER , DIMENSION(4) :: itest 99 REAL(wp), DIMENSION(jpi,jpj) :: z2d 100 REAL(wp), DIMENSION(jpi,jpj) :: zswitch ! ice indicator 101 REAL(wp), DIMENSION(jpi,jpj) :: zht_i_ini, zat_i_ini, zvt_i_ini !data from namelist or nc file 102 REAL(wp), DIMENSION(jpi,jpj) :: zts_u_ini, zht_s_ini, zsm_i_ini, ztm_i_ini !data from namelist or nc file 103 REAL(wp), DIMENSION(jpi,jpj,jpl) :: zh_i_ini, za_i_ini !data by cattegories to fill 102 104 !-------------------------------------------------------------------- 103 105 … … 119 121 CALL eos_fzp( sss_m(:,:), t_bo(:,:) ) 120 122 t_bo(:,:) = ( t_bo(:,:) + rt0 ) * tmask(:,:,1) 121 122 123 123 124 !-------------------------------------------------------------------- … … 377 378 ! Melt pond volume and fraction 378 379 IF ( ln_pnd ) THEN 379 380 380 DO jl = 1, jpl 381 382 381 a_ip_frac(:,:,jl) = 0.2 * zswitch(:,:) 383 382 h_ip (:,:,jl) = 0.05 * zswitch(:,:) 384 383 a_ip(:,:,jl) = a_ip_frac(:,:,jl) * a_i (:,:,jl) 385 384 v_ip(:,:,jl) = h_ip (:,:,jl) * a_ip(:,:,jl) 386 387 END DO 388 385 END DO 389 386 ELSE 390 391 387 a_ip(:,:,:) = 0._wp 392 388 v_ip(:,:,:) = 0._wp 393 389 a_ip_frac(:,:,:) = 0._wp 394 390 h_ip (:,:,:) = 0._wp 395 396 391 ENDIF 397 392 ! END MV MP 2016 … … 437 432 u_ice (:,:) = 0._wp 438 433 v_ice (:,:) = 0._wp 439 stress1_i(:,:) = 0._wp 440 stress2_i(:,:) = 0._wp 441 stress12_i(:,:) = 0._wp 442 443 !-------------------------------------------------------------------- 444 ! 5) Moments for advection 445 !-------------------------------------------------------------------- 446 447 sxopw (:,:) = 0._wp 448 syopw (:,:) = 0._wp 449 sxxopw(:,:) = 0._wp 450 syyopw(:,:) = 0._wp 451 sxyopw(:,:) = 0._wp 452 453 sxice (:,:,:) = 0._wp ; sxsn (:,:,:) = 0._wp ; sxa (:,:,:) = 0._wp 454 syice (:,:,:) = 0._wp ; sysn (:,:,:) = 0._wp ; sya (:,:,:) = 0._wp 455 sxxice(:,:,:) = 0._wp ; sxxsn(:,:,:) = 0._wp ; sxxa (:,:,:) = 0._wp 456 syyice(:,:,:) = 0._wp ; syysn(:,:,:) = 0._wp ; syya (:,:,:) = 0._wp 457 sxyice(:,:,:) = 0._wp ; sxysn(:,:,:) = 0._wp ; sxya (:,:,:) = 0._wp 458 459 sxc0 (:,:,:) = 0._wp ; sxe (:,:,:,:)= 0._wp 460 syc0 (:,:,:) = 0._wp ; sye (:,:,:,:)= 0._wp 461 sxxc0 (:,:,:) = 0._wp ; sxxe (:,:,:,:)= 0._wp 462 syyc0 (:,:,:) = 0._wp ; syye (:,:,:,:)= 0._wp 463 sxyc0 (:,:,:) = 0._wp ; sxye (:,:,:,:)= 0._wp 464 465 sxsal (:,:,:) = 0._wp 466 sysal (:,:,:) = 0._wp 467 sxxsal (:,:,:) = 0._wp 468 syysal (:,:,:) = 0._wp 469 sxysal (:,:,:) = 0._wp 470 471 sxage (:,:,:) = 0._wp 472 syage (:,:,:) = 0._wp 473 sxxage (:,:,:) = 0._wp 474 syyage (:,:,:) = 0._wp 475 sxyage (:,:,:) = 0._wp 476 434 ! 435 !-------------------------------------------------------------------- 436 ! Snow-ice mass (case ice is fully embedded) | 437 !-------------------------------------------------------------------- 438 snwice_mass (:,:) = tmask(:,:,1) * ( rhosn * vt_s(:,:) + rhoic * vt_i(:,:) ) ! snow+ice mass 439 snwice_mass_b(:,:) = snwice_mass(:,:) 440 ! 441 IF( ln_ice_embd ) THEN ! embedded sea-ice: deplete the initial ssh below sea-ice area 442 443 sshn(:,:) = sshn(:,:) - snwice_mass(:,:) * r1_rau0 444 sshb(:,:) = sshb(:,:) - snwice_mass(:,:) * r1_rau0 445 446 IF( .NOT.ln_linssh ) THEN 447 448 WHERE( ht_0(:,:) > 0 ) ; z2d(:,:) = 1._wp + sshn(:,:)*tmask(:,:,1) / ht_0(:,:) 449 ELSEWHERE ; z2d(:,:) = 1._wp ; END WHERE 450 451 DO jk = 1,jpkm1 ! adjust initial vertical scale factors 452 e3t_n(:,:,jk) = e3t_0(:,:,jk) * z2d(:,:) 453 e3t_b(:,:,jk) = e3t_n(:,:,jk) 454 e3t_a(:,:,jk) = e3t_n(:,:,jk) 455 END DO 456 457 ! Reconstruction of all vertical scale factors at now and before time-steps 458 ! ========================================================================= 459 ! Horizontal scale factor interpolations 460 ! -------------------------------------- 461 CALL dom_vvl_interpol( e3t_b(:,:,:), e3u_b(:,:,:), 'U' ) 462 CALL dom_vvl_interpol( e3t_b(:,:,:), e3v_b(:,:,:), 'V' ) 463 CALL dom_vvl_interpol( e3t_n(:,:,:), e3u_n(:,:,:), 'U' ) 464 CALL dom_vvl_interpol( e3t_n(:,:,:), e3v_n(:,:,:), 'V' ) 465 CALL dom_vvl_interpol( e3u_n(:,:,:), e3f_n(:,:,:), 'F' ) 466 ! Vertical scale factor interpolations 467 ! ------------------------------------ 468 CALL dom_vvl_interpol( e3t_n(:,:,:), e3w_n (:,:,:), 'W' ) 469 CALL dom_vvl_interpol( e3u_n(:,:,:), e3uw_n(:,:,:), 'UW' ) 470 CALL dom_vvl_interpol( e3v_n(:,:,:), e3vw_n(:,:,:), 'VW' ) 471 CALL dom_vvl_interpol( e3u_b(:,:,:), e3uw_b(:,:,:), 'UW' ) 472 CALL dom_vvl_interpol( e3v_b(:,:,:), e3vw_b(:,:,:), 'VW' ) 473 ! t- and w- points depth 474 ! ---------------------- 475 !!gm not sure of that.... 476 gdept_n(:,:,1) = 0.5_wp * e3w_n(:,:,1) 477 gdepw_n(:,:,1) = 0.0_wp 478 gde3w_n(:,:,1) = gdept_n(:,:,1) - sshn(:,:) 479 DO jk = 2, jpk 480 gdept_n(:,:,jk) = gdept_n(:,:,jk-1) + e3w_n(:,:,jk ) 481 gdepw_n(:,:,jk) = gdepw_n(:,:,jk-1) + e3t_n(:,:,jk-1) 482 gde3w_n(:,:,jk) = gdept_n(:,:,jk ) - sshn (:,:) 483 END DO 484 ENDIF 485 ENDIF 486 477 487 !------------------------------------ 478 488 ! 6) store fields at before time-step … … 488 498 u_ice_b(:,:) = u_ice(:,:) 489 499 v_ice_b(:,:) = v_ice(:,:) 490 491 ! MV MP 2016492 IF ( nn_pnd_scheme > 0 ) THEN493 sxap (:,:,:) = 0._wp ; sxvp (:,:,:) = 0._wp494 syap (:,:,:) = 0._wp ; syvp (:,:,:) = 0._wp495 sxxap (:,:,:) = 0._wp ; sxxvp (:,:,:) = 0._wp496 syyap (:,:,:) = 0._wp ; syyvp (:,:,:) = 0._wp497 sxyap (:,:,:) = 0._wp ; sxyvp (:,:,:) = 0._wp498 ENDIF499 ! END MV MP 2016500 500 501 501 !!!clem -
branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icerdgrft.F90
r8517 r8518 34 34 PUBLIC ice_rdgrft_strength ! called by icerhg_evp 35 35 PUBLIC ice_rdgrft_init ! called by ice_stp 36 PUBLIC ice_rdgrft_alloc ! called by ice_init37 36 38 37 ! Variables shared among ridging subroutines 39 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: asum ! sum of total ice and open water area40 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: aksum ! ratio of area removed to area ridged41 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: athorn ! participation function; fraction of ridging/closing associated w/ category n42 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: hrmin ! minimum ridge thickness43 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: hrmax ! maximum ridge thickness44 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: hraft ! thickness of rafted ice45 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: krdg ! thickness of ridging ice / mean ridge thickness46 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: aridge ! participating ice ridging47 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: araft ! participating ice rafting38 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: asum ! sum of total ice and open water area 39 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: aksum ! ratio of area removed to area ridged 40 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: athorn ! participation function; fraction of ridging/closing associated w/ category n 41 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: hrmin ! minimum ridge thickness 42 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: hrmax ! maximum ridge thickness 43 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: hraft ! thickness of rafted ice 44 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: krdg ! thickness of ridging ice / mean ridge thickness 45 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: aridge ! participating ice ridging 46 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: araft ! participating ice rafting 48 47 ! 49 48 REAL(wp), PARAMETER :: krdgmin = 1.1_wp ! min ridge thickness multiplier … … 52 51 ! 53 52 ! ** namelist (namice_rdgrft) ** 53 LOGICAL :: ln_str_H79 ! ice strength parameterization (Hibler79) 54 REAL(wp) :: rn_pstar ! determines ice strength, Hibler JPO79 55 REAL(wp) :: rn_crhg ! determines changes in ice strength 56 LOGICAL :: ln_str_R75 ! ice strength parameterization (Rothrock75) 57 REAL(wp) :: rn_perdg ! ridging work divided by pot. energy change in ridging 54 58 REAL(wp) :: rn_csrdg ! fraction of shearing energy contributing to ridging 55 59 LOGICAL :: ln_partf_lin ! participation function linear (Thorndike et al. (1975)) … … 82 86 & hrmin(jpi,jpj,jpl) , hraft (jpi,jpj,jpl) , aridge(jpi,jpj,jpl) , & 83 87 & hrmax(jpi,jpj,jpl) , krdg (jpi,jpj,jpl) , araft (jpi,jpj,jpl) , STAT=ice_rdgrft_alloc ) 84 ! 88 89 IF( lk_mpp ) CALL mpp_sum ( ice_rdgrft_alloc ) 85 90 IF( ice_rdgrft_alloc /= 0 ) CALL ctl_warn( 'ice_rdgrft_alloc: failed to allocate arrays' ) 86 91 ! … … 909 914 CALL ctl_stop( 'ice_rdgrft_init: choose one and only one participation function (ln_partf_lin or ln_partf_exp)' ) 910 915 ENDIF 916 ! ! allocate tke arrays 917 IF( ice_rdgrft_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'ice_rdgrft_init : unable to allocate arrays' ) 911 918 ! 912 919 END SUBROUTINE ice_rdgrft_init -
branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icerhg.F90
r8517 r8518 37 37 INTEGER, PARAMETER :: np_rhgEVP = 1 ! EVP rheology 38 38 !! INTEGER, PARAMETER :: np_rhgEAP = 2 ! EAP rheology 39 39 40 ! ** namelist (namrhg) ** 41 LOGICAL :: ln_rhg_EVP ! EVP rheology 42 ! 40 43 !! * Substitutions 41 44 # include "vectopt_loop_substitute.h90" … … 77 80 ! !------------------------! 78 81 CALL ice_rhg_evp( kt, stress1_i, stress2_i, stress12_i, u_ice, v_ice, shear_i, divu_i, delta_i ) 79 82 ! 80 83 END SELECT 84 ! 85 IF( lrst_ice ) THEN !* write EVP fields in the restart file 86 IF( ln_rhg_EVP ) CALL rhg_evp_rst( 'WRITE', kt ) 87 ENDIF 81 88 ! 82 89 ! controls … … 86 93 ! 87 94 END SUBROUTINE ice_rhg 88 89 95 90 96 SUBROUTINE ice_rhg_init … … 131 137 !! IF( ln_rhg_EAP ) THEN ; ioptio = ioptio + 1 ; nice_rhg = np_rhgEAP ; ENDIF 132 138 IF( ioptio /= 1 ) CALL ctl_stop( 'ice_rhg_init: choose one and only one ice rheology' ) 133 134 ! ! allocate tke arrays 135 !!clem example IF( zdf_tke_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'zdf_tke_init : unable to allocate arrays' ) 139 ! 140 IF( ln_rhg_EVP ) CALL rhg_evp_rst( 'READ' ) !* read or initialize all required files 136 141 ! 137 142 END SUBROUTINE ice_rhg_init -
branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icerhg_evp.F90
r8517 r8518 41 41 PRIVATE 42 42 43 PUBLIC ice_rhg_evp ! routine called by icerhg.F90 43 PUBLIC ice_rhg_evp ! called by icerhg.F90 44 PUBLIC rhg_evp_rst ! called by icerhg.F90 44 45 45 46 !! * Substitutions … … 108 109 !!------------------------------------------------------------------- 109 110 INTEGER, INTENT(in) :: kt ! time step 110 REAL(wp), DIMENSION( jpi,jpj), INTENT(inout) :: pstress1_i, pstress2_i, pstress12_i111 REAL(wp), DIMENSION( jpi,jpj), INTENT( out) :: pu_ice, pv_ice, pshear_i, pdivu_i, pdelta_i111 REAL(wp), DIMENSION(:,:), INTENT(inout) :: pstress1_i, pstress2_i, pstress12_i 112 REAL(wp), DIMENSION(:,:), INTENT( out) :: pu_ice, pv_ice, pshear_i, pdivu_i, pdelta_i 112 113 !! 113 114 INTEGER :: ji, jj ! dummy loop indices … … 825 826 END SUBROUTINE ice_rhg_evp 826 827 828 SUBROUTINE rhg_evp_rst( cdrw, kt ) 829 !!--------------------------------------------------------------------- 830 !! *** ROUTINE rhg_evp_rst *** 831 !! 832 !! ** Purpose : Read or write RHG file in restart file 833 !! 834 !! ** Method : use of IOM library 835 !!---------------------------------------------------------------------- 836 CHARACTER(len=*) , INTENT(in) :: cdrw ! "READ"/"WRITE" flag 837 INTEGER, OPTIONAL, INTENT(in) :: kt ! ice time-step 838 ! 839 INTEGER :: iter ! local integer 840 INTEGER :: id1, id2, id3 ! local integers 841 !!---------------------------------------------------------------------- 842 ! 843 IF( TRIM(cdrw) == 'READ' ) THEN ! Read/initialize 844 ! ! --------------- 845 IF( ln_rstart ) THEN !* Read the restart file 846 ! 847 id1 = iom_varid( numrir, 'stress1_i' , ldstop = .FALSE. ) 848 id2 = iom_varid( numrir, 'stress2_i' , ldstop = .FALSE. ) 849 id3 = iom_varid( numrir, 'stress12_i', ldstop = .FALSE. ) 850 ! 851 IF( MIN( id1, id2, id3 ) > 0 ) THEN ! fields exist 852 CALL iom_get( numrir, jpdom_autoglo, 'stress1_i' , stress1_i ) 853 CALL iom_get( numrir, jpdom_autoglo, 'stress2_i' , stress2_i ) 854 CALL iom_get( numrir, jpdom_autoglo, 'stress12_i', stress12_i ) 855 ELSE ! start rheology from rest 856 IF(lwp) WRITE(numout,*) ' ==>> previous run without rheology, set stresses to 0' 857 stress1_i (:,:) = 0._wp 858 stress2_i (:,:) = 0._wp 859 stress12_i(:,:) = 0._wp 860 ENDIF 861 ELSE !* Start from rest 862 IF(lwp) WRITE(numout,*) ' ==>> start from rest: set stresses to 0' 863 stress1_i (:,:) = 0._wp 864 stress2_i (:,:) = 0._wp 865 stress12_i(:,:) = 0._wp 866 ENDIF 867 ! 868 ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN ! Create restart file 869 ! ! ------------------- 870 IF(lwp) WRITE(numout,*) '---- rhg-rst ----' 871 iter = kt + nn_fsbc - 1 ! ice restarts are written at kt == nitrst - nn_fsbc + 1 872 ! 873 CALL iom_rstput( iter, nitrst, numriw, 'stress1_i' , stress1_i ) 874 CALL iom_rstput( iter, nitrst, numriw, 'stress2_i' , stress2_i ) 875 CALL iom_rstput( iter, nitrst, numriw, 'stress12_i', stress12_i ) 876 ! 877 ENDIF 878 ! 879 END SUBROUTINE rhg_evp_rst 880 827 881 #else 828 882 !!---------------------------------------------------------------------- -
branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icerst.F90
r8514 r8518 17 17 !!---------------------------------------------------------------------- 18 18 USE ice ! sea-ice variables 19 USE sbc_ice , ONLY : snwice_mass, snwice_mass_b20 19 USE dom_oce ! ocean domain 21 20 USE sbc_oce , ONLY : nn_fsbc … … 30 29 PRIVATE 31 30 32 PUBLIC ice_rst_opn ! routine called by icestep.F90 33 PUBLIC ice_rst_write ! routine called by icestep.F90 34 PUBLIC ice_rst_read ! routine called by ice_init 35 36 LOGICAL, PUBLIC :: lrst_ice !: logical to control the ice restart write 37 INTEGER, PUBLIC :: numrir, numriw !: logical unit for ice restart (read and write) 31 PUBLIC ice_rst_opn ! called by icestp 32 PUBLIC ice_rst_write ! called by icestp 33 PUBLIC ice_rst_read ! called by ice_init 38 34 39 35 !!---------------------------------------------------------------------- … … 79 75 WRITE(numout,*) ' open ice restart NetCDF file: ',TRIM(clpath)//clname 80 76 END SELECT 81 IF( kt == nitrst - 2*nn_fsbc + 1 ) THEN 82 WRITE(numout,*) ' kt = nitrst - 2*nn_fsbc + 1 = ', kt,' date= ', ndastp 83 ELSE ; WRITE(numout,*) ' kt = ' , kt,' date= ', ndastp 77 IF( kt == nitrst - 2*nn_fsbc + 1 ) THEN 78 WRITE(numout,*) ' kt = nitrst - 2*nn_fsbc + 1 = ', kt,' date= ', ndastp 79 ELSE 80 WRITE(numout,*) ' kt = ' , kt,' date= ', ndastp 84 81 ENDIF 85 82 ENDIF … … 102 99 INTEGER, INTENT(in) :: kt ! number of iteration 103 100 !! 104 INTEGER :: j i, jj, jk ,jl ! dummy loop indices101 INTEGER :: jk ,jl ! dummy loop indices 105 102 INTEGER :: iter 106 103 CHARACTER(len=25) :: znam … … 127 124 !!gm "just" ask Sebatien 128 125 129 130 126 ! Prognostic variables 131 127 DO jl = 1, jpl … … 133 129 znam = 'v_i'//'_htc'//zchar 134 130 z2d(:,:) = v_i(:,:,jl) 135 CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 131 CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) ! v_i 136 132 znam = 'v_s'//'_htc'//zchar 137 133 z2d(:,:) = v_s(:,:,jl) 138 CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 134 CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) ! v_s 139 135 znam = 'smv_i'//'_htc'//zchar 140 136 z2d(:,:) = smv_i(:,:,jl) 141 CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 137 CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) ! smv_i 142 138 znam = 'oa_i'//'_htc'//zchar 143 139 z2d(:,:) = oa_i(:,:,jl) 144 CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 140 CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) ! oa_i 145 141 znam = 'a_i'//'_htc'//zchar 146 142 z2d(:,:) = a_i(:,:,jl) 147 CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 143 CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) ! a_i 148 144 znam = 't_su'//'_htc'//zchar 149 145 z2d(:,:) = t_su(:,:,jl) 150 CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 146 CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) ! t_su 151 147 END DO 152 148 … … 157 153 znam = 'a_ip'//'_htc'//zchar 158 154 z2d(:,:) = a_ip(:,:,jl) 159 CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 155 CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) ! a_ip 160 156 znam = 'v_ip'//'_htc'//zchar 161 157 z2d(:,:) = v_ip(:,:,jl) 162 CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 158 CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) ! v_ip 163 159 END DO 164 160 ENDIF … … 169 165 znam = 'tempt_sl1'//'_htc'//zchar 170 166 z2d(:,:) = e_s(:,:,1,jl) 171 CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 167 CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) ! e_s 172 168 DO jk = 1, nlay_i 173 169 WRITE(zchar1,'(I2.2)') jk 174 170 znam = 'tempt'//'_il'//zchar1//'_htc'//zchar 175 171 z2d(:,:) = e_i(:,:,jk,jl) 176 CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 177 END DO 178 END DO 179 180 CALL iom_rstput( iter, nitrst, numriw, 'u_ice' , u_ice ) 181 CALL iom_rstput( iter, nitrst, numriw, 'v_ice' , v_ice ) 182 CALL iom_rstput( iter, nitrst, numriw, 'stress1_i' , stress1_i ) 183 CALL iom_rstput( iter, nitrst, numriw, 'stress2_i' , stress2_i ) 184 CALL iom_rstput( iter, nitrst, numriw, 'stress12_i' , stress12_i ) 185 CALL iom_rstput( iter, nitrst, numriw, 'snwice_mass' , snwice_mass ) 186 CALL iom_rstput( iter, nitrst, numriw, 'snwice_mass_b', snwice_mass_b ) 187 188 ! In case Prather scheme is used for advection, write second order moments 189 ! ------------------------------------------------------------------------ 190 IF( ln_adv_Pra ) THEN 191 192 DO jl = 1, jpl 193 WRITE(zchar,'(I2.2)') jl 194 znam = 'sxice'//'_htc'//zchar 195 z2d(:,:) = sxice(:,:,jl) 196 CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 197 znam = 'syice'//'_htc'//zchar 198 z2d(:,:) = syice(:,:,jl) 199 CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 200 znam = 'sxxice'//'_htc'//zchar 201 z2d(:,:) = sxxice(:,:,jl) 202 CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 203 znam = 'syyice'//'_htc'//zchar 204 z2d(:,:) = syyice(:,:,jl) 205 CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 206 znam = 'sxyice'//'_htc'//zchar 207 z2d(:,:) = sxyice(:,:,jl) 208 CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 209 znam = 'sxsn'//'_htc'//zchar 210 z2d(:,:) = sxsn(:,:,jl) 211 CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 212 znam = 'sysn'//'_htc'//zchar 213 z2d(:,:) = sysn(:,:,jl) 214 CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 215 znam = 'sxxsn'//'_htc'//zchar 216 z2d(:,:) = sxxsn(:,:,jl) 217 CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 218 znam = 'syysn'//'_htc'//zchar 219 z2d(:,:) = syysn(:,:,jl) 220 CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 221 znam = 'sxysn'//'_htc'//zchar 222 z2d(:,:) = sxysn(:,:,jl) 223 CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 224 znam = 'sxa'//'_htc'//zchar 225 z2d(:,:) = sxa(:,:,jl) 226 CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 227 znam = 'sya'//'_htc'//zchar 228 z2d(:,:) = sya(:,:,jl) 229 CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 230 znam = 'sxxa'//'_htc'//zchar 231 z2d(:,:) = sxxa(:,:,jl) 232 CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 233 znam = 'syya'//'_htc'//zchar 234 z2d(:,:) = syya(:,:,jl) 235 CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 236 znam = 'sxya'//'_htc'//zchar 237 z2d(:,:) = sxya(:,:,jl) 238 CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 239 znam = 'sxc0'//'_htc'//zchar 240 z2d(:,:) = sxc0(:,:,jl) 241 CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 242 znam = 'syc0'//'_htc'//zchar 243 z2d(:,:) = syc0(:,:,jl) 244 CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 245 znam = 'sxxc0'//'_htc'//zchar 246 z2d(:,:) = sxxc0(:,:,jl) 247 CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 248 znam = 'syyc0'//'_htc'//zchar 249 z2d(:,:) = syyc0(:,:,jl) 250 CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 251 znam = 'sxyc0'//'_htc'//zchar 252 z2d(:,:) = sxyc0(:,:,jl) 253 CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 254 znam = 'sxsal'//'_htc'//zchar 255 z2d(:,:) = sxsal(:,:,jl) 256 CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 257 znam = 'sysal'//'_htc'//zchar 258 z2d(:,:) = sysal(:,:,jl) 259 CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 260 znam = 'sxxsal'//'_htc'//zchar 261 z2d(:,:) = sxxsal(:,:,jl) 262 CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 263 znam = 'syysal'//'_htc'//zchar 264 z2d(:,:) = syysal(:,:,jl) 265 CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 266 znam = 'sxysal'//'_htc'//zchar 267 z2d(:,:) = sxysal(:,:,jl) 268 CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 269 znam = 'sxage'//'_htc'//zchar 270 z2d(:,:) = sxage(:,:,jl) 271 CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 272 znam = 'syage'//'_htc'//zchar 273 z2d(:,:) = syage(:,:,jl) 274 CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 275 znam = 'sxxage'//'_htc'//zchar 276 z2d(:,:) = sxxage(:,:,jl) 277 CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 278 znam = 'syyage'//'_htc'//zchar 279 z2d(:,:) = syyage(:,:,jl) 280 CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 281 znam = 'sxyage'//'_htc'//zchar 282 z2d(:,:) = sxyage(:,:,jl) 283 CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 284 END DO 285 286 CALL iom_rstput( iter, nitrst, numriw, 'sxopw ' , sxopw ) 287 CALL iom_rstput( iter, nitrst, numriw, 'syopw ' , syopw ) 288 CALL iom_rstput( iter, nitrst, numriw, 'sxxopw' , sxxopw ) 289 CALL iom_rstput( iter, nitrst, numriw, 'syyopw' , syyopw ) 290 CALL iom_rstput( iter, nitrst, numriw, 'sxyopw' , sxyopw ) 291 292 DO jl = 1, jpl 293 WRITE(zchar,'(I2.2)') jl 294 DO jk = 1, nlay_i 295 WRITE(zchar1,'(I2.2)') jk 296 znam = 'sxe'//'_il'//zchar1//'_htc'//zchar 297 z2d(:,:) = sxe(:,:,jk,jl) 298 CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 299 znam = 'sye'//'_il'//zchar1//'_htc'//zchar 300 z2d(:,:) = sye(:,:,jk,jl) 301 CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 302 znam = 'sxxe'//'_il'//zchar1//'_htc'//zchar 303 z2d(:,:) = sxxe(:,:,jk,jl) 304 CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 305 znam = 'syye'//'_il'//zchar1//'_htc'//zchar 306 z2d(:,:) = syye(:,:,jk,jl) 307 CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 308 znam = 'sxye'//'_il'//zchar1//'_htc'//zchar 309 z2d(:,:) = sxye(:,:,jk,jl) 310 CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 311 END DO 312 END DO 313 ! MV MP 2016 314 IF ( nn_pnd_scheme > 0 ) THEN 315 DO jl = 1, jpl 316 WRITE(zchar,'(I2.2)') jl 317 znam = 'sxap'//'_htc'//zchar 318 z2d(:,:) = sxap(:,:,jl) 319 CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 320 znam = 'syap'//'_htc'//zchar 321 z2d(:,:) = syap(:,:,jl) 322 CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 323 znam = 'sxxap'//'_htc'//zchar 324 z2d(:,:) = sxxap(:,:,jl) 325 CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 326 znam = 'syyap'//'_htc'//zchar 327 z2d(:,:) = syyap(:,:,jl) 328 CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 329 znam = 'sxyap'//'_htc'//zchar 330 z2d(:,:) = sxyap(:,:,jl) 331 CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 332 333 znam = 'sxvp'//'_htc'//zchar 334 z2d(:,:) = sxvp(:,:,jl) 335 CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 336 znam = 'syvp'//'_htc'//zchar 337 z2d(:,:) = syvp(:,:,jl) 338 CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 339 znam = 'sxxvp'//'_htc'//zchar 340 z2d(:,:) = sxxvp(:,:,jl) 341 CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 342 znam = 'syyvp'//'_htc'//zchar 343 z2d(:,:) = syyvp(:,:,jl) 344 CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 345 znam = 'sxyvp'//'_htc'//zchar 346 z2d(:,:) = sxyvp(:,:,jl) 347 CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 348 END DO 349 ENDIF 350 351 ENDIF 172 CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) ! e_i 173 END DO 174 END DO 175 176 CALL iom_rstput( iter, nitrst, numriw, 'u_ice', u_ice ) ! u_ice 177 CALL iom_rstput( iter, nitrst, numriw, 'v_ice', v_ice ) ! v_ice 352 178 353 179 ! close restart file … … 368 194 !! ** purpose : read of sea-ice variable restart in a netcdf file 369 195 !!---------------------------------------------------------------------- 370 INTEGER :: j i, jj, jk, jl196 INTEGER :: jk, jl 371 197 REAL(wp) :: zfice, ziter 372 198 REAL(wp), DIMENSION(jpi,jpj) :: z2d … … 379 205 IF(lwp) THEN 380 206 WRITE(numout,*) 381 WRITE(numout,*) 'ice_rst_read 382 WRITE(numout,*) '~~~~~~~~~~~~ ~'207 WRITE(numout,*) 'ice_rst_read: read ice NetCDF restart file' 208 WRITE(numout,*) '~~~~~~~~~~~~' 383 209 ENDIF 384 210 … … 390 216 IF(lwp) WRITE(numout,*) ' in any case we force it to nit000 - 1 : ', nit000 - 1 391 217 392 !Control of date 393 218 ! Control of date 394 219 IF( ( nit000 - NINT(ziter) ) /= 1 .AND. ABS( nrstdt ) == 1 ) & 395 220 & CALL ctl_stop( 'ice_rst_read ===>>>> : problem with nit000 in ice restart', & … … 451 276 END DO 452 277 453 CALL iom_get( numrir, jpdom_autoglo, 'u_ice' , u_ice ) 454 CALL iom_get( numrir, jpdom_autoglo, 'v_ice' , v_ice ) 455 CALL iom_get( numrir, jpdom_autoglo, 'stress1_i' , stress1_i ) 456 CALL iom_get( numrir, jpdom_autoglo, 'stress2_i' , stress2_i ) 457 CALL iom_get( numrir, jpdom_autoglo, 'stress12_i', stress12_i ) 458 CALL iom_get( numrir, jpdom_autoglo, 'snwice_mass' , snwice_mass ) 459 CALL iom_get( numrir, jpdom_autoglo, 'snwice_mass_b', snwice_mass_b ) 460 461 ! In case Prather scheme is used for advection, read second order moments 462 ! ------------------------------------------------------------------------ 463 IF( ln_adv_Pra ) THEN 464 465 DO jl = 1, jpl 466 WRITE(zchar,'(I2.2)') jl 467 znam = 'sxice'//'_htc'//zchar 468 CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 469 sxice(:,:,jl) = z2d(:,:) 470 znam = 'syice'//'_htc'//zchar 471 CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 472 syice(:,:,jl) = z2d(:,:) 473 znam = 'sxxice'//'_htc'//zchar 474 CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 475 sxxice(:,:,jl) = z2d(:,:) 476 znam = 'syyice'//'_htc'//zchar 477 CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 478 syyice(:,:,jl) = z2d(:,:) 479 znam = 'sxyice'//'_htc'//zchar 480 CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 481 sxyice(:,:,jl) = z2d(:,:) 482 znam = 'sxsn'//'_htc'//zchar 483 CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 484 sxsn(:,:,jl) = z2d(:,:) 485 znam = 'sysn'//'_htc'//zchar 486 CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 487 sysn(:,:,jl) = z2d(:,:) 488 znam = 'sxxsn'//'_htc'//zchar 489 CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 490 sxxsn(:,:,jl) = z2d(:,:) 491 znam = 'syysn'//'_htc'//zchar 492 CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 493 syysn(:,:,jl) = z2d(:,:) 494 znam = 'sxysn'//'_htc'//zchar 495 CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 496 sxysn(:,:,jl) = z2d(:,:) 497 znam = 'sxa'//'_htc'//zchar 498 CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 499 sxa(:,:,jl) = z2d(:,:) 500 znam = 'sya'//'_htc'//zchar 501 CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 502 sya(:,:,jl) = z2d(:,:) 503 znam = 'sxxa'//'_htc'//zchar 504 CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 505 sxxa(:,:,jl) = z2d(:,:) 506 znam = 'syya'//'_htc'//zchar 507 CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 508 syya(:,:,jl) = z2d(:,:) 509 znam = 'sxya'//'_htc'//zchar 510 CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 511 sxya(:,:,jl) = z2d(:,:) 512 znam = 'sxc0'//'_htc'//zchar 513 CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 514 sxc0(:,:,jl) = z2d(:,:) 515 znam = 'syc0'//'_htc'//zchar 516 CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 517 syc0(:,:,jl) = z2d(:,:) 518 znam = 'sxxc0'//'_htc'//zchar 519 CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 520 sxxc0(:,:,jl) = z2d(:,:) 521 znam = 'syyc0'//'_htc'//zchar 522 CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 523 syyc0(:,:,jl) = z2d(:,:) 524 znam = 'sxyc0'//'_htc'//zchar 525 CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 526 sxyc0(:,:,jl) = z2d(:,:) 527 znam = 'sxsal'//'_htc'//zchar 528 CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 529 sxsal(:,:,jl) = z2d(:,:) 530 znam = 'sysal'//'_htc'//zchar 531 CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 532 sysal(:,:,jl) = z2d(:,:) 533 znam = 'sxxsal'//'_htc'//zchar 534 CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 535 sxxsal(:,:,jl) = z2d(:,:) 536 znam = 'syysal'//'_htc'//zchar 537 CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 538 syysal(:,:,jl) = z2d(:,:) 539 znam = 'sxysal'//'_htc'//zchar 540 CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 541 sxysal(:,:,jl) = z2d(:,:) 542 znam = 'sxage'//'_htc'//zchar 543 CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 544 sxage(:,:,jl) = z2d(:,:) 545 znam = 'syage'//'_htc'//zchar 546 CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 547 syage(:,:,jl) = z2d(:,:) 548 znam = 'sxxage'//'_htc'//zchar 549 CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 550 sxxage(:,:,jl) = z2d(:,:) 551 znam = 'syyage'//'_htc'//zchar 552 CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 553 syyage(:,:,jl) = z2d(:,:) 554 znam = 'sxyage'//'_htc'//zchar 555 CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 556 sxyage(:,:,jl)= z2d(:,:) 557 END DO 558 ! MV MP 2016 559 IF ( nn_pnd_scheme > 0 ) THEN 560 DO jl = 1, jpl 561 WRITE(zchar,'(I2.2)') jl 562 znam = 'sxap'//'_htc'//zchar 563 CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 564 sxap(:,:,jl) = z2d(:,:) 565 znam = 'syap'//'_htc'//zchar 566 CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 567 syap(:,:,jl) = z2d(:,:) 568 znam = 'sxxap'//'_htc'//zchar 569 CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 570 sxxap(:,:,jl) = z2d(:,:) 571 znam = 'syyap'//'_htc'//zchar 572 CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 573 syyap(:,:,jl) = z2d(:,:) 574 znam = 'sxyap'//'_htc'//zchar 575 CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 576 sxyap(:,:,jl) = z2d(:,:) 577 578 znam = 'sxvp'//'_htc'//zchar 579 CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 580 sxvp(:,:,jl) = z2d(:,:) 581 znam = 'syvp'//'_htc'//zchar 582 CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 583 syvp(:,:,jl) = z2d(:,:) 584 znam = 'sxxvp'//'_htc'//zchar 585 CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 586 sxxvp(:,:,jl) = z2d(:,:) 587 znam = 'syyvp'//'_htc'//zchar 588 CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 589 syyvp(:,:,jl) = z2d(:,:) 590 znam = 'sxyvp'//'_htc'//zchar 591 CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 592 sxyvp(:,:,jl) = z2d(:,:) 593 END DO 594 ENDIF 595 ! END MV MP 2016 596 597 CALL iom_get( numrir, jpdom_autoglo, 'sxopw ' , sxopw ) 598 CALL iom_get( numrir, jpdom_autoglo, 'syopw ' , syopw ) 599 CALL iom_get( numrir, jpdom_autoglo, 'sxxopw' , sxxopw ) 600 CALL iom_get( numrir, jpdom_autoglo, 'syyopw' , syyopw ) 601 CALL iom_get( numrir, jpdom_autoglo, 'sxyopw' , sxyopw ) 602 603 DO jl = 1, jpl 604 WRITE(zchar,'(I2.2)') jl 605 DO jk = 1, nlay_i 606 WRITE(zchar1,'(I2.2)') jk 607 znam = 'sxe'//'_il'//zchar1//'_htc'//zchar 608 CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 609 sxe(:,:,jk,jl) = z2d(:,:) 610 znam = 'sye'//'_il'//zchar1//'_htc'//zchar 611 CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 612 sye(:,:,jk,jl) = z2d(:,:) 613 znam = 'sxxe'//'_il'//zchar1//'_htc'//zchar 614 CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 615 sxxe(:,:,jk,jl) = z2d(:,:) 616 znam = 'syye'//'_il'//zchar1//'_htc'//zchar 617 CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 618 syye(:,:,jk,jl) = z2d(:,:) 619 znam = 'sxye'//'_il'//zchar1//'_htc'//zchar 620 CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 621 sxye(:,:,jk,jl) = z2d(:,:) 622 END DO 623 END DO 624 ! 625 END IF 626 627 ! clem: I do not understand why the following IF is needed 628 ! I suspect something inconsistent in the main code with option nn_icesal=1 629 IF( nn_icesal == 1 ) THEN 630 DO jl = 1, jpl 631 sm_i(:,:,jl) = rn_icesal 632 DO jk = 1, nlay_i 633 s_i(:,:,jk,jl) = rn_icesal 634 END DO 635 END DO 636 ENDIF 637 ! 638 !CALL iom_close( numrir ) !clem: closed in icestp.F90 639 ! 278 CALL iom_get( numrir, jpdom_autoglo, 'u_ice', u_ice ) 279 CALL iom_get( numrir, jpdom_autoglo, 'v_ice', v_ice ) 280 640 281 END SUBROUTINE ice_rst_read 641 282 -
branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icestp.F90
r8517 r8518 206 206 CALL ice_wri( 1 ) ! -- Ice outputs 207 207 ! 208 IF( kt == nit000 .AND. ln_rstart ) & !!gm I don't understand the ln_rstart, if needed, add a comment, please209 & CALL iom_close( numrir ) ! close input ice restart file210 208 ! 211 209 IF( lrst_ice ) CALL ice_rst_write( kt ) ! -- Ice restart file … … 251 249 ierr = ierr + sbc_ice_alloc () ! surface forcing 252 250 ierr = ierr + ice1D_alloc () ! thermodynamics 253 ierr = ierr + ice_rdgrft_alloc () ! ridging/rafting254 251 ! 255 252 IF( lk_mpp ) CALL mpp_sum( ierr ) … … 258 255 CALL ice_itd_init ! ice thickness distribution initialization 259 256 ! 260 IF( ln_icedyn ) THEN261 CALL ice_dyn_init ! set ice dynamics parameters262 CALL ice_rdgrft_init ! set ice ridging/rafting parameters263 CALL ice_rhg_init ! set ice rheology parameters264 CALL ice_adv_init ! set ice advection parameters265 ENDIF266 267 IF( ln_icethd ) THEN268 CALL ice_thd_init ! set ice thermodynics parameters269 CALL ice_thd_sal_init ! set ice salinity parameters270 ENDIF271 272 257 ! MV MP 2016 273 CALL lim_mp_init ! set melt ponds parameters 258 CALL lim_mp_init ! set melt ponds parameters (clem: important to be located here) 274 259 ! END MV MP 2016 275 260 ! ! Initial sea-ice state … … 283 268 CALL ice_var_glo2eqv 284 269 ! 270 IF( ln_icedyn ) THEN 271 CALL ice_dyn_init ! set ice dynamics parameters 272 CALL ice_rdgrft_init ! set ice ridging/rafting parameters 273 CALL ice_rhg_init ! set ice rheology parameters 274 CALL ice_adv_init ! set ice advection parameters 275 ENDIF 276 277 IF( ln_icethd ) THEN 278 CALL ice_thd_init ! set ice thermodynics parameters 279 CALL ice_thd_sal_init ! set ice salinity parameters 280 ENDIF 281 ! 285 282 CALL ice_update_init ! ice surface boundary condition 286 283 ! … … 296 293 ELSEWHERE ; rn_amax_2d(:,:) = rn_amax_s ! SH 297 294 END WHERE 295 296 IF( ln_rstart ) CALL iom_close( numrir ) ! close input ice restart file 298 297 ! 299 298 END SUBROUTINE ice_init … … 429 428 hfx_res(:,:) = 0._wp ; hfx_sub(:,:) = 0._wp 430 429 hfx_spr(:,:) = 0._wp ; hfx_dif(:,:) = 0._wp 431 hfx_err (:,:) = 0._wp ; hfx_err_rem(:,:) = 0._wp430 hfx_err_rem(:,:) = 0._wp 432 431 hfx_err_dif(:,:) = 0._wp 433 432 wfx_err_sub(:,:) = 0._wp -
branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd.F90
r8517 r8518 413 413 CALL tab_2d_1d( nidx, idxice(1:nidx), hfx_snw_1d (1:nidx), hfx_snw ) 414 414 CALL tab_2d_1d( nidx, idxice(1:nidx), hfx_sub_1d (1:nidx), hfx_sub ) 415 CALL tab_2d_1d( nidx, idxice(1:nidx), hfx_err_1d (1:nidx), hfx_err )416 415 CALL tab_2d_1d( nidx, idxice(1:nidx), hfx_res_1d (1:nidx), hfx_res ) 417 416 CALL tab_2d_1d( nidx, idxice(1:nidx), hfx_err_dif_1d(1:nidx), hfx_err_dif ) … … 499 498 CALL tab_1d_2d( nidx, idxice(1:nidx), hfx_snw_1d (1:nidx), hfx_snw ) 500 499 CALL tab_1d_2d( nidx, idxice(1:nidx), hfx_sub_1d (1:nidx), hfx_sub ) 501 CALL tab_1d_2d( nidx, idxice(1:nidx), hfx_err_1d (1:nidx), hfx_err )502 500 CALL tab_1d_2d( nidx, idxice(1:nidx), hfx_res_1d (1:nidx), hfx_res ) 503 501 CALL tab_1d_2d( nidx, idxice(1:nidx), hfx_err_dif_1d(1:nidx), hfx_err_dif ) … … 585 583 ENDIF 586 584 ! 587 IF ( rn_hinew < rn_himin ) CALL ctl_stop( ' STOP', 'ice_thd_init : rn_hinew should be >= rn_himin' )585 IF ( rn_hinew < rn_himin ) CALL ctl_stop( 'ice_thd_init : rn_hinew should be >= rn_himin' ) 588 586 ! 589 587 IF(lwp) WRITE(numout,*) -
branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd_dif.F90
r8514 r8518 94 94 REAL(wp) :: ztsu_err = 1.e-5_wp ! range around which t_su is considered at 0C 95 95 REAL(wp) :: ztmelt_i ! ice melting temperature 96 REAL(wp) :: z hsu96 REAL(wp) :: z1_hsu 97 97 REAL(wp) :: zdti_max ! current maximal error on temperature 98 98 REAL(wp) :: zdti_bnd = 1.e-4_wp ! maximal authorized error on temperature 99 REAL(wp) :: zcpi ! Ice specific heat 100 REAL(wp) :: zhi ! 101 REAL(wp) :: zhfx_err, zdq ! diag errors on heat 99 102 100 103 REAL(wp), DIMENSION(jpij) :: isnow ! switch for presence (1) or absence (0) of snow 101 REAL(wp), DIMENSION(jpij) :: ztsub ! old surface temperature (before the iterative procedure ) 102 REAL(wp), DIMENSION(jpij) :: ztsubit ! surface temperature at previous iteration 104 REAL(wp), DIMENSION(jpij) :: ztsub ! surface temperature at previous iteration 103 105 REAL(wp), DIMENSION(jpij) :: zh_i ! ice layer thickness 104 106 REAL(wp), DIMENSION(jpij) :: zh_s ! snow layer thickness … … 106 108 REAL(wp), DIMENSION(jpij) :: zqns_ice_b ! solar radiation absorbed at the surface 107 109 REAL(wp), DIMENSION(jpij) :: zf ! surface flux function 108 REAL(wp), DIMENSION(jpij) :: dzf! derivative of the surface flux function110 REAL(wp), DIMENSION(jpij) :: zdqns_ice_b ! derivative of the surface flux function 109 111 REAL(wp), DIMENSION(jpij) :: zdti ! current error on temperature 110 REAL(wp), DIMENSION(jpij) :: zdifcase ! case of the equation resolution (1->4)111 112 REAL(wp), DIMENSION(jpij) :: zftrice ! solar radiation transmitted through the ice 112 REAL(wp), DIMENSION(jpij) :: zihic113 113 114 REAL(wp), DIMENSION(jpij,nlay_i) :: z_i ! Vertical cotes of the layers in the ice 115 REAL(wp), DIMENSION(jpij,nlay_s) :: z_s ! Vertical cotes of the layers in the snow 116 REAL(wp), DIMENSION(jpij,nlay_i) :: ztib ! Old temperature in the ice 117 REAL(wp), DIMENSION(jpij,nlay_s) :: ztsb ! Temporary temperature in the snow 118 REAL(wp), DIMENSION(jpij,nlay_i) :: ztitemp ! Temporary temperature in the ice to check the convergence 119 REAL(wp), DIMENSION(jpij,nlay_s) :: ztstemp ! Temporary temperature in the snow to check the convergence 114 120 REAL(wp), DIMENSION(jpij,0:nlay_i) :: ztcond_i ! Ice thermal conductivity 115 121 REAL(wp), DIMENSION(jpij,0:nlay_i) :: zradtr_i ! Radiation transmitted through the ice 116 122 REAL(wp), DIMENSION(jpij,0:nlay_i) :: zradab_i ! Radiation absorbed in the ice 117 123 REAL(wp), DIMENSION(jpij,0:nlay_i) :: zkappa_i ! Kappa factor in the ice 118 REAL(wp), DIMENSION(jpij,0:nlay_i) :: ztib ! Old temperature in the ice119 124 REAL(wp), DIMENSION(jpij,0:nlay_i) :: zeta_i ! Eta factor in the ice 120 REAL(wp), DIMENSION(jpij,0:nlay_i) :: ztitemp ! Temporary temperature in the ice to check the convergence121 REAL(wp), DIMENSION(jpij,0:nlay_i) :: zspeche_i ! Ice specific heat122 REAL(wp), DIMENSION(jpij,0:nlay_i) :: z_i ! Vertical cotes of the layers in the ice123 125 REAL(wp), DIMENSION(jpij,0:nlay_s) :: zradtr_s ! Radiation transmited through the snow 124 126 REAL(wp), DIMENSION(jpij,0:nlay_s) :: zradab_s ! Radiation absorbed in the snow 125 127 REAL(wp), DIMENSION(jpij,0:nlay_s) :: zkappa_s ! Kappa factor in the snow 126 128 REAL(wp), DIMENSION(jpij,0:nlay_s) :: zeta_s ! Eta factor in the snow 127 REAL(wp), DIMENSION(jpij,0:nlay_s) :: ztstemp ! Temporary temperature in the snow to check the convergence128 REAL(wp), DIMENSION(jpij,0:nlay_s) :: ztsb ! Temporary temperature in the snow129 REAL(wp), DIMENSION(jpij,0:nlay_s) :: z_s ! Vertical cotes of the layers in the snow130 129 REAL(wp), DIMENSION(jpij,nlay_i+3) :: zindterm ! 'Ind'ependent term 131 130 REAL(wp), DIMENSION(jpij,nlay_i+3) :: zindtbis ! Temporary 'ind'ependent term 132 131 REAL(wp), DIMENSION(jpij,nlay_i+3) :: zdiagbis ! Temporary 'dia'gonal term 133 132 REAL(wp), DIMENSION(jpij,nlay_i+3,3) :: ztrid ! Tridiagonal system terms 134 REAL(wp), DIMENSION(jpij) :: z dq, zq_ini, zhfx_err! diag errors on heat133 REAL(wp), DIMENSION(jpij) :: zq_ini ! diag errors on heat 135 134 REAL(wp), DIMENSION(jpij) :: zghe ! G(he), th. conduct enhancement factor, mono-cat 136 135 … … 148 147 149 148 ! --- diag error on heat diffusion - PART 1 --- ! 150 zdq(:) = 0._wp ; zq_ini(:) = 0._wp151 149 DO ji = 1, nidx 152 150 zq_ini(ji) = ( SUM( e_i_1d(ji,1:nlay_i) ) * ht_i_1d(ji) * r1_nlay_i + & … … 167 165 ! Ice / snow layers 168 166 !-------------------- 169 170 z_s(:,0) = 0._wp ! vert. coord. of the up. lim. of the 1st snow layer 171 z_i(:,0) = 0._wp ! vert. coord. of the up. lim. of the 1st ice layer 172 173 DO jk = 1, nlay_s ! vert. coord of the up. lim. of the layer-th snow layer 174 DO ji = 1 , nidx 175 z_s(ji,jk) = z_s(ji,jk-1) + ht_s_1d(ji) * r1_nlay_s 176 END DO 177 END DO 178 179 DO jk = 1, nlay_i ! vert. coord of the up. lim. of the layer-th ice layer 180 DO ji = 1 , nidx 181 z_i(ji,jk) = z_i(ji,jk-1) + ht_i_1d(ji) * r1_nlay_i 167 z_s(1:nidx,1) = zh_s(1:nidx) 168 z_i(1:nidx,1) = zh_i(1:nidx) 169 DO jk = 2, nlay_s ! vert. coord of the up. lim. of the layer-th snow layer 170 DO ji = 1 , nidx 171 z_s(ji,jk) = z_s(ji,jk-1) + zh_s(ji) 172 END DO 173 END DO 174 175 DO jk = 2, nlay_i ! vert. coord of the up. lim. of the layer-th ice layer 176 DO ji = 1 , nidx 177 z_i(ji,jk) = z_i(ji,jk-1) + zh_i(ji) 182 178 END DO 183 179 END DO … … 187 183 !------------------------------------------------------------------------------| 188 184 ! 189 !------------------- 190 ! Computation of i0 191 !------------------- 192 ! i0 describes the fraction of solar radiation which does not contribute 193 ! to the surface energy budget but rather penetrates inside the ice. 194 ! We assume that no radiation is transmitted through the snow 195 ! If there is no no snow 196 ! zfsw = (1-i0).qsr_ice is absorbed at the surface 197 ! zftrice = io.qsr_ice is below the surface 198 ! ftr_ice = io.qsr_ice.exp(-k(h_i)) transmitted below the ice 199 ! fr1_i0_1d = i0 for a thin ice cover, fr1_i0_2d = i0 for a thick ice cover 200 zhsu = 0.1_wp ! threshold for the computation of i0 185 z1_hsu = 1._wp / 0.1_wp ! threshold for the computation of i0 201 186 DO ji = 1 , nidx 202 ! switches 203 isnow(ji) = 1._wp - MAX( 0._wp , SIGN( 1._wp , - ht_s_1d(ji) ) ) 204 ! hs > 0, isnow = 1 205 zihic(ji) = MAX( 0._wp , 1._wp - ( ht_i_1d(ji) / zhsu ) ) 206 207 i0(ji) = ( 1._wp - isnow(ji) ) * ( fr1_i0_1d(ji) + zihic(ji) * fr2_i0_1d(ji) ) 208 END DO 209 210 !------------------------------------------------------- 211 ! Solar radiation absorbed / transmitted at the surface 212 ! Derivative of the non solar flux 213 !------------------------------------------------------- 214 DO ji = 1 , nidx 215 zfsw (ji) = qsr_ice_1d(ji) * ( 1 - i0(ji) ) ! Shortwave radiation absorbed at surface 216 zftrice(ji) = qsr_ice_1d(ji) * i0(ji) ! Solar radiation transmitted below the surface layer 217 dzf (ji) = dqns_ice_1d(ji) ! derivative of incoming nonsolar flux 218 zqns_ice_b(ji) = qns_ice_1d(ji) ! store previous qns_ice_1d value 187 !------------------- 188 ! Computation of i0 189 !------------------- 190 ! i0 describes the fraction of solar radiation which does not contribute 191 ! to the surface energy budget but rather penetrates inside the ice. 192 ! We assume that no radiation is transmitted through the snow 193 ! If there is no no snow 194 ! zfsw = (1-i0).qsr_ice is absorbed at the surface 195 ! zftrice = io.qsr_ice is below the surface 196 ! ftr_ice = io.qsr_ice.exp(-k(h_i)) transmitted below the ice 197 ! fr1_i0_1d = i0 for a thin ice cover, fr1_i0_2d = i0 for a thick ice cover 198 zhi = MAX( 0._wp , 1._wp - ( ht_i_1d(ji) * z1_hsu ) ) 199 i0(ji) = ( 1._wp - isnow(ji) ) * ( fr1_i0_1d(ji) + zhi * fr2_i0_1d(ji) ) 200 201 !------------------------------------------------------- 202 ! Solar radiation absorbed / transmitted at the surface 203 ! Derivative of the non solar flux 204 !------------------------------------------------------- 205 zfsw (ji) = qsr_ice_1d(ji) * ( 1 - i0(ji) ) ! Shortwave radiation absorbed at surface 206 zftrice(ji) = qsr_ice_1d(ji) * i0(ji) ! Solar radiation transmitted below the surface layer 207 zdqns_ice_b(ji) = dqns_ice_1d(ji) ! derivative of incoming nonsolar flux 208 zqns_ice_b (ji) = qns_ice_1d(ji) ! store previous qns_ice_1d value 219 209 END DO 220 210 … … 222 212 ! Transmission - absorption of solar radiation in the ice 223 213 !--------------------------------------------------------- 224 225 DO ji = 1, nidx ! snow initialization 226 zradtr_s(ji,0) = zftrice(ji) ! radiation penetrating through snow 227 END DO 228 229 DO jk = 1, nlay_s ! Radiation through snow 214 zradtr_s(1:nidx,0) = zftrice(1:nidx) 215 DO jk = 1, nlay_s 230 216 DO ji = 1, nidx 231 217 ! ! radiation transmitted below the layer-th snow layer 232 zradtr_s(ji,jk) = zradtr_s(ji,0) * EXP( - zraext_s * ( MAX ( 0._wp , z_s(ji,jk) )) )218 zradtr_s(ji,jk) = zradtr_s(ji,0) * EXP( - zraext_s * z_s(ji,jk) ) 233 219 ! ! radiation absorbed by the layer-th snow layer 234 220 zradab_s(ji,jk) = zradtr_s(ji,jk-1) - zradtr_s(ji,jk) 235 221 END DO 236 222 END DO 237 238 DO ji = 1, nidx ! ice initialization 239 zradtr_i(ji,0) = zradtr_s(ji,nlay_s) * isnow(ji) + zftrice(ji) * ( 1._wp - isnow(ji) ) 240 END DO 241 242 DO jk = 1, nlay_i ! Radiation through ice 223 224 zradtr_i(1:nidx,0) = zradtr_s(1:nidx,nlay_s) * isnow(1:nidx) + zftrice(1:nidx) * ( 1._wp - isnow(1:nidx) ) 225 DO jk = 1, nlay_i 243 226 DO ji = 1, nidx 244 227 ! ! radiation transmitted below the layer-th ice layer 245 zradtr_i(ji,jk) = zradtr_i(ji,0) * EXP( - rn_kappa_i * ( MAX ( 0._wp , z_i(ji,jk) )) )228 zradtr_i(ji,jk) = zradtr_i(ji,0) * EXP( - rn_kappa_i * z_i(ji,jk) ) 246 229 ! ! radiation absorbed by the layer-th ice layer 247 230 zradab_i(ji,jk) = zradtr_i(ji,jk-1) - zradtr_i(ji,jk) … … 249 232 END DO 250 233 251 DO ji = 1, nidx ! Radiation transmitted below the ice 252 ftr_ice_1d(ji) = zradtr_i(ji,nlay_i) 253 END DO 234 ftr_ice_1d(1:nidx) = zradtr_i(1:nidx,nlay_i) ! record radiation transmitted below the ice 254 235 255 236 !------------------------------------------------------------------------------| … … 257 238 !------------------------------------------------------------------------------| 258 239 ! 259 DO ji = 1, nidx ! Old surface temperature 260 ztsub (ji) = t_su_1d(ji) ! temperature at the beg of iter pr. 261 ztsubit(ji) = t_su_1d(ji) ! temperature at the previous iter 262 t_su_1d(ji) = MIN( t_su_1d(ji), rt0 - ztsu_err ) ! necessary 263 zdti (ji) = 1000._wp ! initial value of error 264 END DO 265 266 DO jk = 1, nlay_s ! Old snow temperature 267 DO ji = 1 , nidx 268 ztsb(ji,jk) = t_s_1d(ji,jk) 269 END DO 270 END DO 271 272 DO jk = 1, nlay_i ! Old ice temperature 273 DO ji = 1 , nidx 274 ztib(ji,jk) = t_i_1d(ji,jk) 275 END DO 276 END DO 240 ztsub (1:nidx) = t_su_1d(1:nidx) ! temperature at the previous iter 241 t_su_1d(1:nidx) = MIN( t_su_1d(1:nidx), rt0 - ztsu_err ) ! necessary 242 zdti (1:nidx) = 1000._wp ! initial value of error 243 ztsb (1:nidx,:) = t_s_1d(1:nidx,:) ! Old snow temperature 244 ztib (1:nidx,:) = t_i_1d(1:nidx,:) ! Old ice temperature 277 245 278 246 iconv = 0 ! number of iterations … … 289 257 IF( ln_cndi_U64 ) THEN !-- Untersteiner (1964) formula: k = k0 + beta.S/T 290 258 DO ji = 1 , nidx 291 ztcond_i(ji,0) = rcdic + zbeta * s_i_1d(ji,1) / MIN( -epsi10, t_i_1d(ji,1) - rt0 ) 292 ztcond_i(ji,0) = MAX( ztcond_i(ji,0), zkimin ) 259 ztcond_i(ji,0) = MAX( zkimin, rcdic + zbeta * s_i_1d(ji,1) / MIN( -epsi10, t_i_1d(ji,1) - rt0 ) ) 293 260 END DO 294 261 DO jk = 1, nlay_i-1 295 262 DO ji = 1 , nidx 296 ztcond_i(ji,jk) = rcdic + zbeta * ( s_i_1d(ji,jk) + s_i_1d(ji,jk+1) ) / & 297 MIN(-2.0_wp * epsi10, t_i_1d(ji,jk) + t_i_1d(ji,jk+1) - 2.0_wp * rt0) 298 ztcond_i(ji,jk) = MAX( ztcond_i(ji,jk), zkimin ) 263 ztcond_i(ji,jk) = MAX( zkimin, rcdic + zbeta * ( s_i_1d(ji,jk) + s_i_1d(ji,jk+1) ) / & 264 & MIN(-2.0_wp * epsi10, t_i_1d(ji,jk) + t_i_1d(ji,jk+1) - 2.0_wp * rt0) ) 299 265 END DO 300 266 END DO … … 302 268 ELSEIF( ln_cndi_P07 ) THEN !-- Pringle et al formula: k = k0 + beta1.S/T - beta2.T 303 269 DO ji = 1 , nidx 304 ztcond_i(ji,0) = rcdic + 0.090_wp * s_i_1d(ji,1) / MIN( -epsi10, t_i_1d(ji,1) - rt0 ) & 305 & - 0.011_wp * ( t_i_1d(ji,1) - rt0 ) 306 ztcond_i(ji,0) = MAX( ztcond_i(ji,0), zkimin ) 270 ztcond_i(ji,0) = MAX( zkimin, rcdic + 0.090_wp * s_i_1d(ji,1) / MIN( -epsi10, t_i_1d(ji,1) - rt0 ) & 271 & - 0.011_wp * ( t_i_1d(ji,1) - rt0 ) ) 307 272 END DO 308 273 DO jk = 1, nlay_i-1 309 274 DO ji = 1 , nidx 310 ztcond_i(ji,jk) = rcdic +&275 ztcond_i(ji,jk) = MAX( zkimin, rcdic + & 311 276 & 0.09_wp * ( s_i_1d(ji,jk) + s_i_1d(ji,jk+1) ) & 312 277 & / MIN( -2._wp * epsi10, t_i_1d(ji,jk) + t_i_1d(ji,jk+1) - 2.0_wp * rt0 ) & 313 & - 0.0055_wp * ( t_i_1d(ji,jk) + t_i_1d(ji,jk+1) - 2.0 * rt0 ) 314 ztcond_i(ji,jk) = MAX( ztcond_i(ji,jk), zkimin ) 278 & - 0.0055_wp * ( t_i_1d(ji,jk) + t_i_1d(ji,jk+1) - 2.0 * rt0 ) ) 315 279 END DO 316 280 END DO 317 281 DO ji = 1 , nidx 318 ztcond_i(ji,nlay_i) = rcdic + 0.090_wp * s_i_1d(ji,nlay_i) / MIN( -epsi10, t_bo_1d(ji) - rt0 ) & 319 & - 0.011_wp * ( t_bo_1d(ji) - rt0 ) 320 ztcond_i(ji,nlay_i) = MAX( ztcond_i(ji,nlay_i), zkimin ) 282 ztcond_i(ji,nlay_i) = MAX( zkimin, rcdic + 0.090_wp * s_i_1d(ji,nlay_i) / MIN( -epsi10, t_bo_1d(ji) - rt0 ) & 283 & - 0.011_wp * ( t_bo_1d(ji) - rt0 ) ) 321 284 END DO 322 285 ENDIF … … 335 298 SELECT CASE ( nn_monocat ) 336 299 337 CASE (1,3) ! LIM3300 CASE (1,3) 338 301 339 302 zepsilon = 0.1_wp … … 404 367 DO jk = 1, nlay_i 405 368 DO ji = 1 , nidx 406 ztitemp(ji,jk) 407 z speche_i(ji,jk)= cpic + zgamma * s_i_1d(ji,jk) / MAX( ( t_i_1d(ji,jk) - rt0 ) * ( ztib(ji,jk) - rt0 ), epsi10 )408 zeta_i(ji,jk) = rdt_ice / MAX( rhoic * zspeche_i(ji,jk)* zh_i(ji), epsi10 )369 ztitemp(ji,jk) = t_i_1d(ji,jk) 370 zcpi = cpic + zgamma * s_i_1d(ji,jk) / MAX( ( t_i_1d(ji,jk) - rt0 ) * ( ztib(ji,jk) - rt0 ), epsi10 ) 371 zeta_i(ji,jk) = rdt_ice / MAX( rhoic * zcpi * zh_i(ji), epsi10 ) 409 372 END DO 410 373 END DO … … 425 388 DO ji = 1 , nidx 426 389 ! update of the non solar flux according to the update in T_su 427 qns_ice_1d(ji) = qns_ice_1d(ji) + dqns_ice_1d(ji) * ( t_su_1d(ji) - ztsub it(ji) )390 qns_ice_1d(ji) = qns_ice_1d(ji) + dqns_ice_1d(ji) * ( t_su_1d(ji) - ztsub(ji) ) 428 391 END DO 429 392 ENDIF … … 507 470 ! case 1 : no surface melting - snow present | 508 471 !------------------------------------------------------------------------------| 509 zdifcase(ji) = 1.0510 472 numeqmin(ji) = 1 511 473 numeqmax(ji) = nlay_i + nlay_s + 1 … … 513 475 !!surface equation 514 476 ztrid(ji,1,1) = 0.0 515 ztrid(ji,1,2) = dzf(ji) - zg1s * zkappa_s(ji,0)477 ztrid(ji,1,2) = zdqns_ice_b(ji) - zg1s * zkappa_s(ji,0) 516 478 ztrid(ji,1,3) = zg1s * zkappa_s(ji,0) 517 zindterm(ji,1) = dzf(ji) * t_su_1d(ji) - zf(ji)479 zindterm(ji,1) = zdqns_ice_b(ji) * t_su_1d(ji) - zf(ji) 518 480 519 481 !!first layer of snow equation … … 529 491 !------------------------------------------------------------------------------| 530 492 ! 531 zdifcase(ji) = 2.0532 493 numeqmin(ji) = 2 533 494 numeqmax(ji) = nlay_i + nlay_s + 1 … … 552 513 !------------------------------------------------------------------------------| 553 514 ! 554 zdifcase(ji) = 3.0555 515 numeqmin(ji) = nlay_s + 1 556 516 numeqmax(ji) = nlay_i + nlay_s + 1 … … 558 518 !!surface equation 559 519 ztrid(ji,numeqmin(ji),1) = 0.0 560 ztrid(ji,numeqmin(ji),2) = dzf(ji) - zkappa_i(ji,0)*zg1520 ztrid(ji,numeqmin(ji),2) = zdqns_ice_b(ji) - zkappa_i(ji,0)*zg1 561 521 ztrid(ji,numeqmin(ji),3) = zkappa_i(ji,0)*zg1 562 zindterm(ji,numeqmin(ji)) = dzf(ji)*t_su_1d(ji) - zf(ji)522 zindterm(ji,numeqmin(ji)) = zdqns_ice_b(ji)*t_su_1d(ji) - zf(ji) 563 523 564 524 !!first layer of ice equation … … 572 532 IF ( nlay_i == 1 ) THEN 573 533 ztrid(ji,numeqmin(ji),1) = 0.0 574 ztrid(ji,numeqmin(ji),2) = dzf(ji) - zkappa_i(ji,0) * 2.0534 ztrid(ji,numeqmin(ji),2) = zdqns_ice_b(ji) - zkappa_i(ji,0) * 2.0 575 535 ztrid(ji,numeqmin(ji),3) = zkappa_i(ji,0) * 2.0 576 536 ztrid(ji,numeqmin(ji)+1,1) = -zkappa_i(ji,0) * 2.0 * zeta_i(ji,1) … … 589 549 !------------------------------------------------------------------------------| 590 550 ! 591 zdifcase(ji) = 4.0592 551 numeqmin(ji) = nlay_s + 2 593 552 numeqmax(ji) = nlay_i + nlay_s + 1 … … 662 621 ! surface temperature 663 622 isnow(ji) = 1._wp - MAX( 0._wp , SIGN( 1._wp , -ht_s_1d(ji) ) ) 664 ztsub it(ji) = t_su_1d(ji)623 ztsub(ji) = t_su_1d(ji) 665 624 IF( t_su_1d(ji) < rt0 ) & 666 625 t_su_1d(ji) = ( zindtbis(ji,numeqmin(ji)) - ztrid(ji,numeqmin(ji),3) * & … … 676 635 DO ji = 1 , nidx 677 636 t_su_1d(ji) = MAX( MIN( t_su_1d(ji) , rt0 ) , 190._wp ) 678 zdti (ji) = ABS( t_su_1d(ji) - ztsub it(ji) )637 zdti (ji) = ABS( t_su_1d(ji) - ztsub(ji) ) 679 638 END DO 680 639 … … 754 713 755 714 ! --- diag conservation imbalance on heat diffusion - PART 2 --- ! 715 ! hfx_dif = Heat flux used to warm/cool ice in W.m-2 716 ! zhfx_err = correction on the diagnosed heat flux due to non-convergence of the algorithm used to solve heat equation 756 717 DO ji = 1, nidx 757 zdq(ji) = - zq_ini(ji) + ( SUM( e_i_1d(ji,1:nlay_i) ) * ht_i_1d(ji) * r1_nlay_i + & 758 & SUM( e_s_1d(ji,1:nlay_s) ) * ht_s_1d(ji) * r1_nlay_s ) 718 zdq = - zq_ini(ji) + ( SUM( e_i_1d(ji,1:nlay_i) ) * ht_i_1d(ji) * r1_nlay_i + & 719 & SUM( e_s_1d(ji,1:nlay_s) ) * ht_s_1d(ji) * r1_nlay_s ) 720 759 721 IF( t_su_1d(ji) < rt0 ) THEN ! case T_su < 0degC 760 zhfx_err (ji) = qns_ice_1d(ji) + qsr_ice_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) + zdq(ji) * r1_rdtice722 zhfx_err = ( qns_ice_1d(ji) + qsr_ice_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) + zdq * r1_rdtice ) * a_i_1d(ji) 761 723 ELSE ! case T_su = 0degC 762 zhfx_err (ji) = fc_su(ji) + i0(ji) * qsr_ice_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) + zdq(ji) * r1_rdtice724 zhfx_err = ( fc_su(ji) + i0(ji) * qsr_ice_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) + zdq * r1_rdtice ) * a_i_1d(ji) 763 725 ENDIF 764 hfx_ err_1d(ji) = hfx_err_1d(ji) + zhfx_err(ji)* a_i_1d(ji)726 hfx_dif_1d(ji) = hfx_dif_1d(ji) - zdq * r1_rdtice * a_i_1d(ji) 765 727 766 728 ! total heat that is sent to the ocean (i.e. not used in the heat diffusion equation) 767 hfx_err_dif_1d(ji) = hfx_err_dif_1d(ji) + zhfx_err(ji) * a_i_1d(ji) 768 END DO 769 770 !----------------------------------------- 771 ! Heat flux used to warm/cool ice in W.m-2 772 !----------------------------------------- 773 DO ji = 1, nidx 774 IF( t_su_1d(ji) < rt0 ) THEN ! case T_su < 0degC 775 hfx_dif_1d(ji) = hfx_dif_1d(ji) + & 776 & ( qns_ice_1d(ji) + qsr_ice_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) ) * a_i_1d(ji) 777 ELSE ! case T_su = 0degC 778 hfx_dif_1d(ji) = hfx_dif_1d(ji) + & 779 & ( fc_su(ji) + i0(ji) * qsr_ice_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) ) * a_i_1d(ji) 780 ENDIF 781 ! correction on the diagnosed heat flux due to non-convergence of the algorithm used to solve heat equation 782 hfx_dif_1d(ji) = hfx_dif_1d(ji) - zhfx_err(ji) * a_i_1d(ji) 729 hfx_err_dif_1d(ji) = hfx_err_dif_1d(ji) + zhfx_err 730 783 731 END DO 784 732 ! -
branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd_sal.F90
r8514 r8518 29 29 PUBLIC ice_thd_sal ! called by icethd module 30 30 PUBLIC ice_thd_sal_init ! called by ice_init 31 32 ! ** namelist (namsal) ** 33 LOGICAL :: ln_icedS ! activate gravity drainage and flushing (T) or not (F) 34 REAL(wp) :: rn_sal_gd ! restoring salinity for gravity drainage [PSU] 35 REAL(wp) :: rn_time_gd ! restoring time constant for gravity drainage (= 20 days) [s] 36 REAL(wp) :: rn_sal_fl ! restoring salinity for flushing [PSU] 37 REAL(wp) :: rn_time_fl ! restoring time constant for gravity drainage (= 10 days) [s] 31 38 32 39 !!---------------------------------------------------------------------- … … 115 122 !! 116 123 NAMELIST/namice_sal/ ln_icedS , nn_icesal , rn_icesal, rn_sal_gd, rn_time_gd, & 117 & rn_sal_fl, rn_time_fl, rn_simax , rn_simin124 & rn_sal_fl, rn_time_fl, rn_simax , rn_simin 118 125 !!------------------------------------------------------------------- 119 126 ! … … 143 150 ENDIF 144 151 ! 152 IF( ln_icedS .AND. nn_icesal == 1 ) THEN 153 ln_icedS = .FALSE. 154 CALL ctl_warn('ln_icedS is set to false since constant ice salinity is chosen (nn_icesal=1)') 155 ENDIF 156 ! 145 157 END SUBROUTINE ice_thd_sal_init 146 158 -
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 -
branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icewri.F90
r8500 r8518 192 192 CALL iom_put ('hfxsnw' , hfx_snw(:,:) ) ! 193 193 CALL iom_put ('hfxsub' , hfx_sub(:,:) ) ! 194 CALL iom_put ('hfxerr' , hfx_err (:,:)) !194 CALL iom_put ('hfxerr' , hfx_err_dif(:,:) ) ! 195 195 CALL iom_put ('hfxerr_rem' , hfx_err_rem(:,:) ) ! 196 196 -
branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/OPA_SRC/IOM/in_out_manager.F90
r8226 r8518 88 88 INTEGER :: nitrst !: time step at which restart file should be written 89 89 LOGICAL :: lrst_oce !: logical to control the oce restart write 90 LOGICAL :: lrst_ice !: logical to control the ice restart write 90 91 INTEGER :: numror = 0 !: logical unit for ocean restart (read). Init to 0 is needed for SAS (in daymod.F90) 92 INTEGER :: numrir !: logical unit for ice restart (read) 91 93 INTEGER :: numrow !: logical unit for ocean restart (write) 94 INTEGER :: numriw !: logical unit for ice restart (write) 92 95 INTEGER :: nrst_lst !: number of restart to output next 93 96
Note: See TracChangeset
for help on using the changeset viewer.