Changeset 3604 for trunk/NEMOGCM/NEMO/OPA_SRC/ASM/asminc.F90
- Timestamp:
- 2012-11-19T15:21:34+01:00 (11 years ago)
- File:
-
- 1 edited
-
trunk/NEMOGCM/NEMO/OPA_SRC/ASM/asminc.F90 (modified) (29 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/NEMOGCM/NEMO/OPA_SRC/ASM/asminc.F90
r3294 r3604 10 10 !! NEMO 3.3 ! 2010-05 (D. Lea) Update to work with NEMO v3.2 11 11 !! - ! 2010-05 (D. Lea) add calc_month_len routine based on day_init 12 !! 3.4 ! 2012-10 (A. Weaver and K. Mogensen) Fix for direct initialization 12 13 !!---------------------------------------------------------------------- 13 14 … … 20 21 !! dyn_asm_inc : Apply the dynamic (u and v) increments 21 22 !! ssh_asm_inc : Apply the SSH increment 23 !! seaice_asm_inc : Apply the seaice increment 22 24 !!---------------------------------------------------------------------- 23 25 USE wrk_nemo ! Memory Allocation … … 25 27 USE dom_oce ! Ocean space and time domain 26 28 USE oce ! Dynamics and active tracers defined in memory 27 USE divcur ! Horizontal divergence and relative vorticity28 29 USE ldfdyn_oce ! ocean dynamics: lateral physics 29 30 USE eosbn2 ! Equation of state - in situ and potential density … … 33 34 USE c1d ! 1D initialization 34 35 USE in_out_manager ! I/O manager 35 USE lib_mpp ! MPP library 36 USE lib_mpp ! MPP library 37 #if defined key_lim3 38 USE ice ! LIM3 39 #endif 40 #if defined key_lim2 41 USE ice_2 ! LIM2 42 #endif 43 USE sbc_oce ! Surface boundary condition variables. 36 44 37 45 IMPLICIT NONE … … 43 51 PUBLIC dyn_asm_inc !: Apply the dynamic (u and v) increments 44 52 PUBLIC ssh_asm_inc !: Apply the SSH increment 53 PUBLIC seaice_asm_inc !: Apply the seaice increment 45 54 46 55 #if defined key_asminc … … 50 59 #endif 51 60 LOGICAL, PUBLIC :: ln_bkgwri = .FALSE. !: No output of the background state fields 52 LOGICAL, PUBLIC :: ln_trjwri = .FALSE. !: No output of the state trajectory fields53 61 LOGICAL, PUBLIC :: ln_asmiau = .FALSE. !: No applying forcing with an assimilation increment 54 62 LOGICAL, PUBLIC :: ln_asmdin = .FALSE. !: No direct initialization … … 56 64 LOGICAL, PUBLIC :: ln_dyninc = .FALSE. !: No dynamics (u and v) assimilation increments 57 65 LOGICAL, PUBLIC :: ln_sshinc = .FALSE. !: No sea surface height assimilation increment 66 LOGICAL, PUBLIC :: ln_seaiceinc = .FALSE. !: No sea ice concentration increment 58 67 LOGICAL, PUBLIC :: ln_salfix = .FALSE. !: Apply minimum salinity check 68 LOGICAL, PUBLIC :: ln_temnofreeze = .FALSE. !: Don't allow the temperature to drop below freezing 59 69 INTEGER, PUBLIC :: nn_divdmp = 0 !: Apply divergence damping filter nn_divdmp times 60 70 … … 78 88 79 89 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: ssh_bkg, ssh_bkginc ! Background sea surface height and its increment 90 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: seaice_bkginc ! Increment to the background sea ice conc 80 91 81 92 !! * Substitutions … … 125 136 REAL(wp), POINTER, DIMENSION(:,:) :: hdiv 126 137 !! 127 NAMELIST/nam_asminc/ ln_bkgwri, ln_trjwri,&138 NAMELIST/nam_asminc/ ln_bkgwri, & 128 139 & ln_trainc, ln_dyninc, ln_sshinc, & 129 140 & ln_asmdin, ln_asmiau, & 130 141 & nitbkg, nitdin, nitiaustr, nitiaufin, niaufn, & 131 & nittrjfrq,ln_salfix, salfixmin, &142 & ln_salfix, salfixmin, & 132 143 & nn_divdmp 133 144 !!---------------------------------------------------------------------- … … 139 150 ! Set default values 140 151 ln_bkgwri = .FALSE. 141 ln_trjwri = .FALSE.142 152 ln_trainc = .FALSE. 143 153 ln_dyninc = .FALSE. 144 154 ln_sshinc = .FALSE. 155 ln_seaiceinc = .FALSE. 145 156 ln_asmdin = .FALSE. 146 157 ln_asmiau = .TRUE. 147 158 ln_salfix = .FALSE. 159 ln_temnofreeze = .FALSE. 148 160 salfixmin = -9999 149 161 nitbkg = 0 … … 152 164 nitiaufin = 150 ! = 10 days with ORCA2 153 165 niaufn = 0 154 nittrjfrq = 1155 166 156 167 REWIND ( numnam ) … … 164 175 WRITE(numout,*) ' Namelist namasm : set assimilation increment parameters' 165 176 WRITE(numout,*) ' Logical switch for writing out background state ln_bkgwri = ', ln_bkgwri 166 WRITE(numout,*) ' Logical switch for writing out state trajectory ln_trjwri = ', ln_trjwri167 177 WRITE(numout,*) ' Logical switch for applying tracer increments ln_trainc = ', ln_trainc 168 178 WRITE(numout,*) ' Logical switch for applying velocity increments ln_dyninc = ', ln_dyninc 169 179 WRITE(numout,*) ' Logical switch for applying SSH increments ln_sshinc = ', ln_sshinc 170 180 WRITE(numout,*) ' Logical switch for Direct Initialization (DI) ln_asmdin = ', ln_asmdin 181 WRITE(numout,*) ' Logical switch for applying sea ice increments ln_seaiceinc = ', ln_seaiceinc 171 182 WRITE(numout,*) ' Logical switch for Incremental Analysis Updating (IAU) ln_asmiau = ', ln_asmiau 172 183 WRITE(numout,*) ' Timestep of background in [0,nitend-nit000-1] nitbkg = ', nitbkg … … 175 186 WRITE(numout,*) ' Timestep of end of IAU interval in [0,nitend-nit000-1] nitiaufin = ', nitiaufin 176 187 WRITE(numout,*) ' Type of IAU weighting function niaufn = ', niaufn 177 WRITE(numout,*) ' Frequency of trajectory output for 4D-VAR nittrjfrq = ', nittrjfrq178 188 WRITE(numout,*) ' Logical switch for ensuring that the sa > salfixmin ln_salfix = ', ln_salfix 179 189 WRITE(numout,*) ' Minimum salinity after applying the increments salfixmin = ', salfixmin … … 213 223 WRITE(numout,*) ' nitiaustr_r = ', nitiaustr_r 214 224 WRITE(numout,*) ' nitiaufin_r = ', nitiaufin_r 215 WRITE(numout,*) ' nittrjfrq = ', nittrjfrq216 225 WRITE(numout,*) 217 226 WRITE(numout,*) ' Dates referenced to current cycle:' … … 235 244 236 245 IF ( ( ( .NOT. ln_asmdin ).AND.( .NOT. ln_asmiau ) ) & 237 .AND.( ( ln_trainc ).OR.( ln_dyninc ).OR.( ln_sshinc ) )) &238 & CALL ctl_stop( ' One or more of ln_trainc, ln_dyninc and ln_sshinc is set to .true.', &246 .AND.( ( ln_trainc ).OR.( ln_dyninc ).OR.( ln_sshinc ) .OR. ( ln_seaiceinc) )) & 247 & CALL ctl_stop( ' One or more of ln_trainc, ln_dyninc, ln_sshinc and ln_seaiceinc is set to .true.', & 239 248 & ' but ln_asmdin and ln_asmiau are both set to .false. :', & 240 249 & ' Inconsistent options') … … 248 257 & ' Type IAU weighting function is invalid') 249 258 250 IF ( ( .NOT. ln_trainc ).AND.( .NOT. ln_dyninc ).AND.( .NOT. ln_sshinc ) &259 IF ( ( .NOT. ln_trainc ).AND.( .NOT. ln_dyninc ).AND.( .NOT. ln_sshinc ).AND.( .NOT. ln_seaiceinc ) & 251 260 & ) & 252 & CALL ctl_warn( ' ln_trainc, ln_dyninc and ln_sshinc are set to .false. :', &261 & CALL ctl_warn( ' ln_trainc, ln_dyninc, ln_sshinc and ln_seaiceinc are set to .false. :', & 253 262 & ' The assimilation increments are not applied') 254 263 … … 353 362 ALLOCATE( v_bkginc(jpi,jpj,jpk) ) 354 363 ALLOCATE( ssh_bkginc(jpi,jpj) ) 364 ALLOCATE( seaice_bkginc(jpi,jpj)) 355 365 #if defined key_asminc 356 366 ALLOCATE( ssh_iau(jpi,jpj) ) … … 361 371 v_bkginc(:,:,:) = 0.0 362 372 ssh_bkginc(:,:) = 0.0 373 seaice_bkginc(:,:) = 0.0 363 374 #if defined key_asminc 364 375 ssh_iau(:,:) = 0.0 365 376 #endif 366 IF ( ( ln_trainc ).OR.( ln_dyninc ).OR.( ln_sshinc ) ) THEN377 IF ( ( ln_trainc ).OR.( ln_dyninc ).OR.( ln_sshinc ).OR.( ln_seaiceinc ) ) THEN 367 378 368 379 !-------------------------------------------------------------------- … … 429 440 ENDIF 430 441 442 IF ( ln_seaiceinc ) THEN 443 CALL iom_get( inum, jpdom_autoglo, 'bckinseaice', seaice_bkginc, 1 ) 444 ! Apply the masks 445 seaice_bkginc(:,:) = seaice_bkginc(:,:) * tmask(:,:,1) 446 ! Set missing increments to 0.0 rather than 1e+20 447 ! to allow for differences in masks 448 WHERE( ABS( seaice_bkginc(:,:) ) > 1.0e+10 ) seaice_bkginc(:,:) = 0.0 449 ENDIF 450 431 451 CALL iom_close( inum ) 432 452 … … 437 457 !----------------------------------------------------------------------- 438 458 439 440 459 IF ( ln_dyninc .AND. nn_divdmp > 0 ) THEN 441 460 442 CALL wrk_alloc(jpi,jpj,hdiv) 443 444 DO jt = 1, nn_divdmp 445 446 DO jk = 1, jpkm1 447 448 hdiv(:,:) = 0._wp 449 450 DO jj = 2, jpjm1 451 DO ji = fs_2, fs_jpim1 ! vector opt. 452 hdiv(ji,jj) = & 453 ( e2u(ji ,jj)*fse3u(ji ,jj,jk) * u_bkginc(ji ,jj,jk) & 454 - e2u(ji-1,jj)*fse3u(ji-1,jj,jk) * u_bkginc(ji-1,jj,jk) & 455 + e1v(ji,jj )*fse3v(ji,jj ,jk) * v_bkginc(ji,jj ,jk) & 456 - e1v(ji,jj-1)*fse3v(ji,jj-1,jk) * v_bkginc(ji,jj-1,jk) ) & 457 / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) ) 461 CALL wrk_alloc(jpi,jpj,hdiv) 462 463 DO jt = 1, nn_divdmp 464 465 DO jk = 1, jpkm1 466 467 hdiv(:,:) = 0._wp 468 469 DO jj = 2, jpjm1 470 DO ji = fs_2, fs_jpim1 ! vector opt. 471 hdiv(ji,jj) = & 472 ( e2u(ji ,jj ) * fse3u(ji ,jj ,jk) * u_bkginc(ji ,jj ,jk) & 473 - e2u(ji-1,jj ) * fse3u(ji-1,jj ,jk) * u_bkginc(ji-1,jj ,jk) & 474 + e1v(ji ,jj ) * fse3v(ji ,jj ,jk) * v_bkginc(ji ,jj ,jk) & 475 - e1v(ji ,jj-1) * fse3v(ji ,jj-1,jk) * v_bkginc(ji ,jj-1,jk) ) & 476 / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) ) 477 END DO 458 478 END DO 479 480 CALL lbc_lnk( hdiv, 'T', 1. ) ! lateral boundary cond. (no sign change) 481 482 DO jj = 2, jpjm1 483 DO ji = fs_2, fs_jpim1 ! vector opt. 484 u_bkginc(ji,jj,jk) = u_bkginc(ji,jj,jk) + 0.2_wp * ( e1t(ji+1,jj)*e2t(ji+1,jj) * hdiv(ji+1,jj) & 485 - e1t(ji ,jj)*e2t(ji ,jj) * hdiv(ji ,jj) ) & 486 / e1u(ji,jj) * umask(ji,jj,jk) 487 v_bkginc(ji,jj,jk) = v_bkginc(ji,jj,jk) + 0.2_wp * ( e1t(ji,jj+1)*e2t(ji,jj+1) * hdiv(ji,jj+1) & 488 - e1t(ji,jj )*e2t(ji,jj ) * hdiv(ji,jj ) ) & 489 / e2v(ji,jj) * vmask(ji,jj,jk) 490 END DO 491 END DO 492 459 493 END DO 460 494 461 CALL lbc_lnk( hdiv, 'T', 1. ) ! lateral boundary cond. (no sign change) 462 463 DO jj = 2, jpjm1 464 DO ji = fs_2, fs_jpim1 ! vector opt. 465 u_bkginc(ji,jj,jk) = u_bkginc(ji,jj,jk) + 0.2 * ( e1t(ji+1,jj)*e2t(ji+1,jj) * hdiv(ji+1,jj) & 466 - e1t(ji ,jj)*e2t(ji ,jj) * hdiv(ji ,jj) ) & 467 / e1u(ji,jj) * umask(ji,jj,jk) 468 v_bkginc(ji,jj,jk) = v_bkginc(ji,jj,jk) + 0.2 * ( e1t(ji,jj+1)*e2t(ji,jj+1) * hdiv(ji,jj+1) & 469 - e1t(ji,jj )*e2t(ji,jj ) * hdiv(ji,jj ) ) & 470 / e2v(ji,jj) * vmask(ji,jj,jk) 471 END DO 472 END DO 473 474 END DO 475 476 END DO 477 478 CALL wrk_dealloc(jpi,jpj,hdiv) 495 END DO 496 497 CALL wrk_dealloc(jpi,jpj,hdiv) 479 498 480 499 ENDIF … … 506 525 CALL iom_open( c_asmdin, inum ) 507 526 508 CALL iom_get( inum, ' zdate', zdate_bkg )527 CALL iom_get( inum, 'rdastp', zdate_bkg ) 509 528 510 529 IF(lwp) THEN … … 662 681 INTEGER :: it 663 682 REAL(wp) :: zincwgt ! IAU weight for current time step 664 !!---------------------------------------------------------------------- 683 REAL (wp), DIMENSION(jpi,jpj,jpk) :: fzptnz ! 3d freezing point values 684 !!---------------------------------------------------------------------- 685 686 ! freezing point calculation taken from oc_fz_pt (but calculated for all depths) 687 ! used to prevent the applied increments taking the temperature below the local freezing point 688 689 #if defined key_cice 690 fzptnz(:,:,:) = -1.8_wp 691 #else 692 DO jk = 1, jpk 693 DO jj = 1, jpj 694 DO ji = 1, jpk 695 fzptnz (ji,jj,jk) = ( -0.0575_wp + 1.710523e-3_wp * SQRT( tsn(ji,jj,jk,jp_sal) ) & 696 - 2.154996e-4_wp * tsn(ji,jj,jk,jp_sal) ) * tsn(ji,jj,jk,jp_sal) & 697 - 7.53e-4_wp * fsdepw(ji,jj,jk) ! (pressure in dbar) 698 END DO 699 END DO 700 END DO 701 #endif 665 702 666 703 IF ( ln_asmiau ) THEN … … 684 721 ! Update the tracer tendencies 685 722 DO jk = 1, jpkm1 686 tsa(:,:,jk,jp_tem) = tsa(:,:,jk,jp_tem) + t_bkginc(:,:,jk) * zincwgt 687 tsa(:,:,jk,jp_sal) = tsa(:,:,jk,jp_sal) + s_bkginc(:,:,jk) * zincwgt 723 IF (ln_temnofreeze) THEN 724 ! Do not apply negative increments if the temperature will fall below freezing 725 WHERE(t_bkginc(:,:,jk) > 0.0_wp .OR. & 726 & tsn(:,:,jk,jp_tem) + tsa(:,:,jk,jp_tem) + t_bkginc(:,:,jk) * wgtiau(it) > fzptnz(:,:,jk) ) 727 tsa(:,:,jk,jp_tem) = tsa(:,:,jk,jp_tem) + t_bkginc(:,:,jk) * zincwgt 728 END WHERE 729 ELSE 730 tsa(:,:,jk,jp_tem) = tsa(:,:,jk,jp_tem) + t_bkginc(:,:,jk) * zincwgt 731 ENDIF 732 IF (ln_salfix) THEN 733 ! Do not apply negative increments if the salinity will fall below a specified 734 ! minimum value salfixmin 735 WHERE(s_bkginc(:,:,jk) > 0.0_wp .OR. & 736 & tsn(:,:,jk,jp_sal) + tsa(:,:,jk,jp_sal) + s_bkginc(:,:,jk) * wgtiau(it) > salfixmin ) 737 tsa(:,:,jk,jp_sal) = tsa(:,:,jk,jp_sal) + s_bkginc(:,:,jk) * zincwgt 738 END WHERE 739 ELSE 740 tsa(:,:,jk,jp_sal) = tsa(:,:,jk,jp_sal) + s_bkginc(:,:,jk) * zincwgt 741 ENDIF 688 742 END DO 689 690 ! Salinity fix691 IF (ln_salfix) THEN692 DO jk = 1, jpkm1693 DO jj = 1, jpj694 DO ji= 1, jpi695 tsa(ji,jj,jk,jp_sal) = MAX( tsa(ji,jj,jk,jp_sal), salfixmin )696 END DO697 END DO698 END DO699 ENDIF700 743 701 744 ENDIF … … 718 761 719 762 ! Initialize the now fields with the background + increment 720 tsn(:,:,:,jp_tem) = t_bkg(:,:,:) + t_bkginc(:,:,:) 721 tsn(:,:,:,jp_sal) = s_bkg(:,:,:) + s_bkginc(:,:,:) 722 723 ! Optional salinity fix 763 IF (ln_temnofreeze) THEN 764 ! Do not apply negative increments if the temperature will fall below freezing 765 WHERE(t_bkginc(:,:,:) > 0.0_wp .OR. & 766 & tsn(:,:,:,jp_tem) + t_bkginc(:,:,:) > fzptnz(:,:,:) ) 767 tsn(:,:,:,jp_tem) = t_bkg(:,:,:) + t_bkginc(:,:,:) 768 END WHERE 769 ELSE 770 tsn(:,:,:,jp_tem) = t_bkg(:,:,:) + t_bkginc(:,:,:) 771 ENDIF 724 772 IF (ln_salfix) THEN 725 DO jk = 1, jpkm1 726 DO jj = 1, jpj 727 DO ji= 1, jpi 728 tsn(ji,jj,jk,jp_sal) = MAX( tsn(ji,jj,jk,jp_sal), salfixmin ) 729 END DO 730 END DO 731 END DO 773 ! Do not apply negative increments if the salinity will fall below a specified 774 ! minimum value salfixmin 775 WHERE(s_bkginc(:,:,:) > 0.0_wp .OR. & 776 & tsn(:,:,:,jp_sal) + s_bkginc(:,:,:) > salfixmin ) 777 tsn(:,:,:,jp_sal) = s_bkg(:,:,:) + s_bkginc(:,:,:) 778 END WHERE 779 ELSE 780 tsn(:,:,:,jp_sal) = s_bkg(:,:,:) + s_bkginc(:,:,:) 732 781 ENDIF 733 782 734 tsb(:,:,:,:) = tsn(:,:,:,:) 783 tsb(:,:,:,:) = tsn(:,:,:,:) ! Update before fields 735 784 736 785 CALL eos( tsb, rhd, rhop ) ! Before potential and in situ densities 737 786 738 787 IF( ln_zps .AND. .NOT. lk_c1d ) & 739 & CALL zps_hde( nit000, jpts, tsb, 740 & gtsu, gtsv, rhd, 788 & CALL zps_hde( nit000, jpts, tsb, & ! Partial steps: before horizontal derivative 789 & gtsu, gtsv, rhd, & ! of T, S, rd at the bottom ocean level 741 790 & gru , grv ) 791 792 #if defined key_zdfkpp 793 CALL eos( tsn, rhd ) ! Compute rhd 794 #endif 742 795 743 796 DEALLOCATE( t_bkginc ) … … 748 801 ! 749 802 ENDIF 803 ! Perhaps the following call should be in step 804 IF ( ln_seaiceinc ) CALL seaice_asm_inc ( kt ) ! apply sea ice concentration increment 750 805 ! 751 806 END SUBROUTINE tra_asm_inc … … 817 872 vb(:,:,:) = vn(:,:,:) 818 873 819 CALL div_cur( kt ) ! Compute divergence and curl for now fields820 821 rotb (:,:,:) = rotn (:,:,:) ! Update before fields822 hdivb(:,:,:) = hdivn(:,:,:)823 824 874 DEALLOCATE( u_bkg ) 825 875 DEALLOCATE( v_bkg ) … … 846 896 ! 847 897 INTEGER :: it 898 INTEGER :: jk 848 899 REAL(wp) :: zincwgt ! IAU weight for current time step 849 900 !!---------------------------------------------------------------------- … … 891 942 sshn(:,:) = ssh_bkg(:,:) + ssh_bkginc(:,:) 892 943 893 sshb(:,:) = sshn(:,:) ! Update before fields 944 ! Update before fields 945 sshb(:,:) = sshn(:,:) 946 947 IF( lk_vvl ) THEN 948 DO jk = 1, jpk 949 fse3t_b(:,:,jk) = fse3t_n(:,:,jk) 950 END DO 951 ENDIF 894 952 895 953 DEALLOCATE( ssh_bkg ) … … 902 960 END SUBROUTINE ssh_asm_inc 903 961 962 SUBROUTINE seaice_asm_inc( kt, kindic ) 963 !!---------------------------------------------------------------------- 964 !! *** ROUTINE seaice_asm_inc *** 965 !! 966 !! ** Purpose : Apply the sea ice assimilation increment. 967 !! 968 !! ** Method : Direct initialization or Incremental Analysis Updating. 969 !! 970 !! ** Action : 971 !! 972 !! History : 973 !! ! 07-2011 (D. Lea) Initial version based on ssh_asm_inc 974 !!---------------------------------------------------------------------- 975 976 IMPLICIT NONE 977 978 !! * Arguments 979 INTEGER, INTENT(IN) :: kt ! Current time step 980 INTEGER, OPTIONAL, INTENT(IN) :: kindic ! flag for disabling the deallocation 981 982 !! * Local declarations 983 INTEGER :: it 984 REAL(wp) :: zincwgt ! IAU weight for current time step 985 986 #if defined key_lim3 || defined key_lim2 987 REAL(wp), DIMENSION(jpi,jpj) :: zofrld, zohicif, zseaicendg, zhicifinc ! LIM 988 REAL(wp) :: zhicifmin=0.5_wp ! ice minimum depth in metres 989 990 #endif 991 992 993 IF ( ln_asmiau ) THEN 994 995 !-------------------------------------------------------------------- 996 ! Incremental Analysis Updating 997 !-------------------------------------------------------------------- 998 999 IF ( ( kt >= nitiaustr_r ).AND.( kt <= nitiaufin_r ) ) THEN 1000 1001 it = kt - nit000 + 1 1002 zincwgt = wgtiau(it) ! IAU weight for the current time step 1003 ! note this is not a tendency so should not be divided by rdt (as with the tracer and other increments) 1004 1005 IF(lwp) THEN 1006 WRITE(numout,*) 1007 WRITE(numout,*) 'seaice_asm_inc : sea ice conc IAU at time step = ', & 1008 & kt,' with IAU weight = ', wgtiau(it) 1009 WRITE(numout,*) '~~~~~~~~~~~~' 1010 ENDIF 1011 1012 #if defined key_lim3 || defined key_lim2 1013 1014 zofrld(:,:)=frld(:,:) 1015 zohicif(:,:)=hicif(:,:) 1016 1017 frld = MIN( MAX( frld(:,:) - seaice_bkginc(:,:) * zincwgt, 0.0_wp), 1.0_wp) 1018 pfrld = MIN( MAX( pfrld(:,:) - seaice_bkginc(:,:) * zincwgt, 0.0_wp), 1.0_wp) 1019 fr_i(:,:) = 1.0_wp - frld(:,:) ! adjust ice fraction 1020 1021 zseaicendg(:,:)=zofrld(:,:) - frld(:,:) ! find out actual sea ice nudge applied 1022 1023 ! Nudge sea ice depth to bring it up to a required minimum depth 1024 1025 WHERE( zseaicendg(:,:) > 0.0_wp .AND. hicif(:,:) < zhicifmin ) 1026 zhicifinc(:,:) = (zhicifmin - hicif(:,:)) * zincwgt 1027 ELSEWHERE 1028 zhicifinc(:,:) = 0.0_wp 1029 END WHERE 1030 1031 ! nudge ice depth 1032 hicif(:,:)=hicif(:,:) + zhicifinc(:,:) 1033 phicif(:,:)=phicif(:,:) + zhicifinc(:,:) 1034 1035 ! seaice salinity balancing (to add) 1036 1037 #endif 1038 1039 #if defined key_cice 1040 1041 ! Pass ice increment tendency into CICE 1042 ndaice_da(:,:) = seaice_bkginc(:,:) * zincwgt / rdt 1043 1044 #endif 1045 1046 IF ( kt == nitiaufin_r ) THEN 1047 DEALLOCATE( seaice_bkginc ) 1048 ENDIF 1049 1050 ELSE 1051 1052 #if defined key_cice 1053 1054 ! Zero ice increment tendency into CICE 1055 ndaice_da(:,:) = 0.0_wp 1056 1057 #endif 1058 1059 ENDIF 1060 1061 ELSEIF ( ln_asmdin ) THEN 1062 1063 !-------------------------------------------------------------------- 1064 ! Direct Initialization 1065 !-------------------------------------------------------------------- 1066 1067 IF ( kt == nitdin_r ) THEN 1068 1069 neuler = 0 ! Force Euler forward step 1070 1071 #if defined key_lim3 || defined key_lim2 1072 1073 zofrld(:,:)=frld(:,:) 1074 zohicif(:,:)=hicif(:,:) 1075 1076 ! Initialize the now fields the background + increment 1077 1078 frld(:,:) = MIN( MAX( frld(:,:) - seaice_bkginc(:,:), 0.0_wp), 1.0_wp) 1079 pfrld(:,:) = frld(:,:) 1080 fr_i(:,:) = 1.0_wp - frld(:,:) ! adjust ice fraction 1081 1082 zseaicendg(:,:)=zofrld(:,:) - frld(:,:) ! find out actual sea ice nudge applied 1083 1084 ! Nudge sea ice depth to bring it up to a required minimum depth 1085 1086 WHERE( zseaicendg(:,:) > 0.0_wp .AND. hicif(:,:) < zhicifmin ) 1087 zhicifinc(:,:) = (zhicifmin - hicif(:,:)) * zincwgt 1088 ELSEWHERE 1089 zhicifinc(:,:) = 0.0_wp 1090 END WHERE 1091 1092 ! nudge ice depth 1093 hicif(:,:)=hicif(:,:) + zhicifinc(:,:) 1094 phicif(:,:)=phicif(:,:) 1095 1096 ! seaice salinity balancing (to add) 1097 1098 #endif 1099 1100 #if defined key_cice 1101 1102 ! Pass ice increment tendency into CICE - is this correct? 1103 ndaice_da(:,:) = seaice_bkginc(:,:) / rdt 1104 1105 #endif 1106 IF ( .NOT. PRESENT(kindic) ) THEN 1107 DEALLOCATE( seaice_bkginc ) 1108 END IF 1109 1110 ELSE 1111 1112 #if defined key_cice 1113 1114 ! Zero ice increment tendency into CICE 1115 ndaice_da(:,:) = 0.0_wp 1116 1117 #endif 1118 1119 ENDIF 1120 1121 !#if defined key_lim3 || defined key_lim2 || defined key_cice 1122 ! 1123 ! IF (ln_seaicebal ) THEN 1124 ! !! balancing salinity increments 1125 ! !! simple case from limflx.F90 (doesn't include a mass flux) 1126 ! !! assumption is that as ice concentration is reduced or increased 1127 ! !! the snow and ice depths remain constant 1128 ! !! note that snow is being created where ice concentration is being increased 1129 ! !! - could be more sophisticated and 1130 ! !! not do this (but would need to alter h_snow) 1131 ! 1132 ! usave(:,:,:)=sb(:,:,:) ! use array as a temporary store 1133 ! 1134 ! DO jj = 1, jpj 1135 ! DO ji = 1, jpi 1136 ! ! calculate change in ice and snow mass per unit area 1137 ! ! positive values imply adding salt to the ocean (results from ice formation) 1138 ! ! fwf : ice formation and melting 1139 ! 1140 ! zfons = ( -nfresh_da(ji,jj)*soce + nfsalt_da(ji,jj) )*rdt 1141 ! 1142 ! ! change salinity down to mixed layer depth 1143 ! mld=hmld_kara(ji,jj) 1144 ! 1145 ! ! prevent small mld 1146 ! ! less than 10m can cause salinity instability 1147 ! IF (mld < 10) mld=10 1148 ! 1149 ! ! set to bottom of a level 1150 ! DO jk = jpk-1, 2, -1 1151 ! IF ((mld > gdepw(ji,jj,jk)) .and. (mld < gdepw(ji,jj,jk+1))) THEN 1152 ! mld=gdepw(ji,jj,jk+1) 1153 ! jkmax=jk 1154 ! ENDIF 1155 ! ENDDO 1156 ! 1157 ! ! avoid applying salinity balancing in shallow water or on land 1158 ! ! 1159 ! 1160 ! ! dsal_ocn (psu kg m^-2) / (kg m^-3 * m) 1161 ! 1162 ! dsal_ocn=0.0_wp 1163 ! sal_thresh=5.0_wp ! minimum salinity threshold for salinity balancing 1164 ! 1165 ! if (tmask(ji,jj,1) > 0 .AND. tmask(ji,jj,jkmax) > 0 ) & 1166 ! dsal_ocn = zfons / (rhop(ji,jj,1) * mld) 1167 ! 1168 ! ! put increments in for levels in the mixed layer 1169 ! ! but prevent salinity below a threshold value 1170 ! 1171 ! DO jk = 1, jkmax 1172 ! 1173 ! IF (dsal_ocn > 0.0_wp .or. sb(ji,jj,jk)+dsal_ocn > sal_thresh) THEN 1174 ! sb(ji,jj,jk) = sb(ji,jj,jk) + dsal_ocn 1175 ! sn(ji,jj,jk) = sn(ji,jj,jk) + dsal_ocn 1176 ! ENDIF 1177 ! 1178 ! ENDDO 1179 ! 1180 ! ! ! salt exchanges at the ice/ocean interface 1181 ! ! zpmess = zfons / rdt_ice ! rdt_ice is ice timestep 1182 ! ! 1183 ! !! Adjust fsalt. A +ve fsalt means adding salt to ocean 1184 ! !! fsalt(ji,jj) = fsalt(ji,jj) + zpmess ! adjust fsalt 1185 ! !! 1186 ! !! emps(ji,jj) = emps(ji,jj) + zpmess ! or adjust emps (see icestp1d) 1187 ! !! ! E-P (kg m-2 s-2) 1188 ! ! emp(ji,jj) = emp(ji,jj) + zpmess ! E-P (kg m-2 s-2) 1189 ! ENDDO !ji 1190 ! ENDDO !jj! 1191 ! 1192 ! ENDIF !ln_seaicebal 1193 ! 1194 !#endif 1195 1196 1197 ENDIF 1198 1199 END SUBROUTINE seaice_asm_inc 904 1200 !!====================================================================== 905 1201 END MODULE asminc
Note: See TracChangeset
for help on using the changeset viewer.
