PROGRAM sanity_check_LIM3 ! Checks if a LIM3 restart is physically consistent and outputs the updated ! version. Strongly based on limupdate.F90, and on sanity_check_LIM2 ! (c) Original from Chris König Beatty ! (c) Re-written by Francois Massonnet, November 2012, ! (c) Refreshed for NEMO3.6, Francois Massonnet, April 2016 ! francois.massonnet@uclouvain.be ! ! ! TODO: create subroutines to automaticate the extraction of NetCDF! ! ************************************************************************************************** ! PLAN of the program ! ! ! 0. Header ! 1. Check existence of NetCDF files, grab command line args ! 1) Ice analysis file ! 2) Ice forecast file ! 3) Oce analysis file ! 4) Oce forecast file ! 5) Mask file ! 6) Mesh file ! 2. Get dimensions, nb ice categories, nb ice/snow layers, allocate, get ice thickness bounds ! 3. Load ice+oce and mask+mesh variables from files ! 4. Ice analysis sanity check ! 1) Regularize ice concentration and adapt other variables ! accordingly ! 2) Manually update the snow and ice heat content according to volume ! changes ! 3) Rebin categories with thickness out of bounds (CALL itd_reb() ) ! Note that itd_reb() calls shift_ice if a shift needs be done ! 4) Several ice corrections ! For example, set volume to zero where concentration is zero. ! 5) Shrink ice (CALL shrink_ice() ) ! If sum(a_i) > 1 or < 0 or some a_i > 1 or < 0, then ice is shrunk ! or reset to zero. ! 6) Rebin categories with thickness out of bounds (CALL itd_reb() ) ! We do it once again as thickness may have changed ! during processes 2) and 3) ! 7) Final concentration check (CALL conc_check() ) ! The program stops if > 1 or < 0 occurs in sum or individual categories ! This can ultimately cause trouble in the code as a term (1-sum(a_i))^0.6 ! occurs in a routine => negative power 0.6 is very baaaad ! If ice concentration is just above 1 or just below zero (numerics) then ! the program resets to zero or one. ! 5. Oce analysis sanity check ! 1) Compute the difference between forecast and analyzed sea salinity ! If larger than XXX PSU, bound variations by that value ! Also update the sea surface salinity variable accordingly ! 2) Same for sea temperature; replace analysed temperature by forecast if ! less than -2°C ! 3) Limit ocean velocities to [-2,2] m/s ! 6. Record the files. The original file is copied and the variables are dumped ! in it ! 1) Ice analysis file ! 2) Oce analysis file !*************************************************************************************************** ! ! 0. Header ! --------- USE netcdf USE my_variables IMPLICIT NONE ! Dummy variables INTEGER :: jl ! ! 1. Grab Command Line Arguments ! ------------------------------ IF (iargc()==4) THEN CALL getarg(1, cfile_analysis_ice) CALL getarg(2, cfile_forecast_ice) CALL getarg(3, cfile_analysis_oce) CALL getarg(4, cfile_forecast_oce) ! Check if files exist ! 1) Ice analysis file INQUIRE(FILE=TRIM(cfile_analysis_ice), EXIST=l_ex) IF ( .NOT. l_ex ) THEN WRITE(*,*) "(sanity_check_lim3) File does not exist: " //& &TRIM(cfile_analysis_ice) STOP END IF ! 2) Ice forecast file INQUIRE(FILE=TRIM(cfile_forecast_ice), EXIST=l_ex) IF ( .NOT. l_ex ) THEN WRITE(*,*) "(sanity_check_lim3) File does not exist: " //& &TRIM(cfile_forecast_ice) STOP END IF ! 3) Oce analysis file INQUIRE(FILE=TRIM(cfile_analysis_oce), EXIST=l_ex) IF ( .NOT. l_ex ) THEN WRITE(*,*) "(sanity_check_lim3) File does not exist: " //& &TRIM(cfile_analysis_oce) STOP END IF ! 4) Oce forecast file INQUIRE(FILE=TRIM(cfile_forecast_oce), EXIST=l_ex) IF ( .NOT. l_ex ) THEN WRITE(*,*) "(sanity_check_lim3) File does not exist: " //& &TRIM(cfile_forecast_oce) STOP END IF ! 5) Mask file INQUIRE(FILE=TRIM(cmaskfile), EXIST=l_ex) IF ( .NOT. l_ex ) THEN WRITE(*,*) "(sanity_check_lim3) File does not exist: " //& &TRIM(cmaskfile) STOP END IF ! 6) Mesh file INQUIRE(FILE=TRIM(cmeshfile), EXIST=l_ex) IF ( .NOT. l_ex ) THEN WRITE(*,*) "(sanity_check_lim3) File does not exist: " //& &TRIM(cmeshfile) STOP END IF ! Everything went OK: WRITE(*,*) "(sanity_check_lim3) Starting ..." ELSE ! Write info WRITE(*,*) WRITE(*,*) " sanity_check_LIM3 needs arguments: " WRITE(*,*) " -analysis_file_ice " WRITE(*,*) " -forecast_file_ice " WRITE(*,*) " -analysis_file_oce " WRITE(*,*) " -forecast_file_oce " WRITE(*,*) " Checks NEMO-LIM3 ice and ocean analyses restarst (netcdf) file for sanity and fixes& &them if necessary." WRITE(*,*) WRITE(*,*) " Sanity means for now:" WRITE(*,*) " Strongly follow limupdate.F90" WRITE(*,*) " Files mask.nc and mesh_hgr.nc need to be in the current directory" WRITE(*,*) WRITE(*,*) " Hope to see you again soon." WRITE(*,*) WRITE(*,*) " Chris König Beatty " WRITE(*,*) " Francois Massonnet -- francois.massonnet@uclouvain.be" WRITE(*,*) " Last update: 2013" WRITE(*,*) " Last update: 2016 (to work with NEMO3.6)" STOP "(sanity_check): Stopped." END IF ! iargc ! 2. Get dimensions, and allocate the variables ! --------------------------------------------- CALL get_dimensions() ! Get number of ice categories CALL get_nb_cat() ! Get number of ice layers CALL get_nb_il() ! Get number of snow layers CALL get_nb_sl() ! Allocate variables CALL allocate_variables() ! Register ice thickness limits ! See Clement Rousset LIM3 paper or LIM3 doc: ! http://www.geosci-model-dev.net/8/2991/2015/gmd-8-2991-2015.pdf ! or the routine sbcice_lim.F90 ! ! ! ! zalpha = 0.05 ! exponent of the transform function rn_himean = 2.0 ! the expected mean thickness over the domain zhmax = 3.*rn_himean DO jl = 1, jpl znum = jpl * ( zhmax+1 )**zalpha zden = ( jpl - jl ) * ( zhmax+1 )**zalpha + jl hi_max(jl) = ( znum / zden )**(1./zalpha) - 1 END DO ! Set hi_max(jpl) to a big value to ensure that all ice is thinner than hi_max(jpl) hi_max(jpl) = 99._wp WRITE(*,*), "Ice categories boundaries [m] are", hi_max ! 3. Load variables ! ----------------- CALL load_variables() ! 4. Ice analysis sanity check ! ---------------------------- ! 1) Regularize ice concentration and other variables based on that CALL regularize() ! 2) Heat content update ! F. Massonnet - test CALL update_hc() ! CALL update_hc() ! 3) Rebin categories with thickness out of bounds ! The code here follows closely the contents of subroutine limi_itd_th_reb CALL itd_reb() ! 4) Several ice corrections CALL several_ice_corrections() ! 5) Shrink ice CALL shrink_ice() ! 6) Rebin categories with thickness out of bounds CALL itd_reb() ! 7) Final check. ! Stops if total conc > 1 or < 0, same for inidividual conc ! If exceeds 1 or is below zero for numerical reasons, reset. CALL conc_check() ! 5. Ocean analysis sanity check ! ------------------------------ ! 1) Salinity CALL salinity_check() ! 2) Temperature CALL temperature_check() ! 3) Velocity CALL velocity_check() ! 6. Record NetCDF ! ---------------- ! 1) Ice variables CALL record_ice() ! 2) Oce variables CALL record_oce() END PROGRAM sanity_check_LIM3 ! SUBROUTINES ! ¨¨¨¨¨¨¨¨¨¨¨ SUBROUTINE itd_reb() ! -------------------------------------------------------- ! This routine is strongly inspired from lim_itd_th_reb ! When in situ thickness exceeds bounds, it transfers ice ! to neighbouring categories ! This routine calls shiftice() defined below ! -------------------------------------------------------- USE my_variables IMPLICIT NONE ! Dummy variables INTEGER :: ji, jj, jl WRITE(*,*) 'Entering itd_reb()' ! 4.2.1 Compute in situ ice thickness in the categories (if there's ice) DO jl = 1, jpl DO ji = 1, jpi DO jj = 1 , jpj IF ( a_i(ji,jj,jl) > epsi10 ) THEN ht_i(ji,jj,jl) = v_i(ji,jj,jl) / a_i(ji,jj,jl) ELSE ht_i(ji,jj,jl) = rzero ENDIF END DO ! jj END DO ! ji END DO ! jl ! Print data at particular location before rebinning !WRITE(*,*) 'Before rebinning: ' !WRITE(*,*) 'a_i : ',a_i(jiindx,jjindx,:) !WRITE(*,*) 'ht_i : ',ht_i(jiindx,jjindx,:) !WRITE(*,*) 'Total concentration: ', SUM(a_i(jiindx,jjindx,:)) !WRITE(*,*) 'Total volume : ', SUM(v_i(jiindx,jjindx,:)) ! 4.2.2 Make sure thickness of first category is at least hi_max_typ ! Not sure to understand what to do here ! 4.2.3. If a category is not in bounds, shift the entire area, volume and ! energy to the neighbouring category ! Initialize shift arrays DO jl = 1, jpl idonor(:,:,jl) = 0 ! idonor is the index of the category that is giving ice with respect to ! the running category. ! Example: we are looping over categories. When jl = 3 (third category), ! we notice that this category has way too much ice. Then idonor(ji,jj,3) ! will be 3. Example 2 (jl=3): The fourth category has too much ice ! Then idonor(ji,jj,3) = 4 zdaice(:,:,jl) = 0 zdvice(:,:,jl) = 0 END DO ! 4.2.3.1 Move thin categories up ! Strategy: -loop over all categories up to the last-but-one ! -identify thicknesses in the current category that are too big ! -if applicable, transfer all volume, area and energy to the ! jpl + 1 category DO jl = 1, jpl - 1 zshiftflag = 0 DO ji = 1, jpi DO jj= 1, jpj IF ( a_i(ji,jj,jl) > epsi10 .AND. ht_i(ji,jj,jl) > hi_max(jl) ) THEN IF ( ldebug ) THEN WRITE(*,*), "Found excessive ice thickness in category",jl WRITE(*,*), ht_i(ji,jj,jl), "larger than" , hi_max(jl) WRITE(*,*), "At grid point" , ji, jj END IF zshiftflag = 1 ! There's at least one point in the domain ! with too thick a sea ice in this category idonor(ji,jj,jl) = jl ! The running category is then donor zdaice(ji,jj,jl) = a_i(ji,jj,jl) ! Amount of ice to be transferred zdvice(ji,jj,jl) = v_i(ji,jj,jl) ! (everything) END IF END DO! jj END DO ! ji IF (zshiftflag == 1) THEN ! this is the case if there was at least one ! excessive thickness in the running category ! somewhere on the domain CALL shiftice() idonor(:,:,jl) = 0 ! Reset before we move to next category zdaice(:,:,jl) = rzero zdvice(:,:,jl) = rzero END IF END DO ! jl, categories ! 4.2.3.2 Move thick categories down ! Strategy: -loop over all categories starting from last-but-one ! to the first one ! -identify if the thickness of the category above is smaller ! than the upper limit for the running category ! -if so, transfer everything to the running category DO jl = jpl-1, 1, -1 ! first, last, step zshiftflag = 0 DO ji = 1, jpi DO jj= 1, jpj IF ( a_i(ji,jj,jl+1) > epsi10 .AND. ht_i(ji,jj,jl+1) <= hi_max(jl) ) THEN IF ( ldebug ) THEN WRITE(*,*), "Found too small ice thickness in category",jl+1 WRITE(*,*), ht_i(ji,jj,jl+1), "smaller than" , hi_max(jl) END IF zshiftflag = 1 ! There's at least one point in the domain ! with too thin a sea ice in this category idonor(ji,jj,jl) = jl+1 ! The jl+1 category is then donor zdaice(ji,jj,jl) = a_i(ji,jj,jl+1) ! Amount of ice to be transferred zdvice(ji,jj,jl) = v_i(ji,jj,jl+1) ! (everything) END IF END DO! jj END DO ! ji IF (zshiftflag == 1) THEN ! this is the case if there was at least one ! too small thickness in the jl+1 category CALL shiftice() idonor(:,:,jl) = rzero ! Reset before we move to next category zdaice(:,:,jl) = rzero zdvice(:,:,jl) = rzero END IF END DO ! jl, categories ! 4.2.3.3 Conservation check ! Print data at particular location after rebinning !WRITE(*,*) 'After rebinning: ' !WRITE(*,*) 'a_i : ',a_i(jiindx,jjindx,:) !WRITE(*,*) 'ht_i : ',ht_i(jiindx,jjindx,:) !WRITE(*,*) 'Total concentration: ', SUM(a_i(jiindx,jjindx,:)) !WRITE(*,*) 'Total volume : ', SUM(v_i(jiindx,jjindx,:)) WRITE(*,*) 'Leaving itd_reb()' END SUBROUTINE itd_reb SUBROUTINE shiftice() !------------------------------------------------------------- ! This routine is (strongly) inspired by lim_itd_shiftice ! It re-arranges ice thickness in the categories ! ! This routine can potentially be called 2*4 times, ! for the first jpl-1 categories to upgrade their too thick ice ! and for the jpl-1 last categories to downgrade their too thin ice ! ! ! Francois Massonnet UCL 2012 francois.massonnet@uclouvain.be !------------------------------------------------------------- ! Variable declaration and importation USE my_variables IMPLICIT NONE LOGICAL :: zdaice_negative, zdvice_negative, zdaice_greater_aicen,& &zdvice_greater_vicen ! Number of cells with ice to transfer REAL(wp), DIMENSION(jpi,jpj) :: zworka REAL(wp), DIMENSION(jpi,jpj,jpl) :: zaTsfn REAL(wp) :: zdvsnow, zdesnow, zdo_aice ! Volume of snow transferred, ! Snow heat, ice age ! local variables REAL(wp) :: zdsm_vice, zdaTsf, zdeice ! ice age, surface temperature, ! ice heat content ! local variables ! Dummy variables INTEGER :: ji, jj, jl, jk, jl1, jl2 ! End definitions ! ---------------------------------------------------------------------- ! Welcome the user WRITE(*,*) 'Entering SUBROUTINE shiftice' ! Define surface temperature times concentration ! This has more dimensions like energy. It will be used ! when surface temperature will be "transferred" after rebinning DO jl = 1 , jpl zaTsfn(:,:,jl) = a_i(:,:,jl) * t_su(:,:,jl) END DO ! Subroutine DO jl = 1 , jpl - 1 zdaice_negative = .FALSE. zdvice_negative = .FALSE. zdaice_greater_aicen = .FALSE. zdvice_greater_vicen = .FALSE. DO ji = 1 , jpi DO jj = 1 , jpj !--------------------------------------------------------------------- ! 1) Check for daice or dvice out of bounds and reset if necessary !--------------------------------------------------------------------- IF ( idonor(ji,jj,jl) .GT. 0 ) THEN ! If the running category is giving jl1=idonor(ji,jj,jl) ! record the donor category !WRITE(*,*), "At grid point ", ji, jj, "category", jl1, "is giving ice" ! Tackle the cases with problems. Normally, zdaice and zdvice should ! be positive, but ... ! Check for negative transfers of ice given in input ! 1. Ice area IF (zdaice(ji,jj,jl) .LT. rzero ) THEN IF (zdaice(ji,jj,jl) .GT. -epsi10 ) THEN WRITE(*,*)," zdaice is negative but artificial" IF ( ( jl1 .EQ. jl .AND. ht_i(ji,jj,jl1) .GT. hi_max(jl) )& &.OR.& &( jl1 .EQ. jl+1 .AND. ht_i(ji,jj,jl1) .LE. hi_max(jl) )) THEN ! You are here if one of these 2 conditions are verified ! 1) The running category is the donor and has too thick a ! sea ice ! 2) The running category is one category below the donor, ! which has too thin a sea ice ! ! If you still with me, it means that a transfer needs to be ! done but the amount of ice is negative due to roundoff ! error or something. Let us reset zdaice and zdvice to the ! value of a_i and v_i in the category for the transfer zdaice(ji,jj,jl) = a_i(ji,jj,jl1) zdvice(ji,jj,jl) = v_i(ji,jj,jl1) ELSE ! None of the two conditions above is verified ! That is, ! 1) The running category is not the donor or has a correct ! ice ! AND ! 2) The category above the running one is not the donor, or ! the ice is correct in this above category ! ! Since everything is fine, nothing should be done zdaice(ji,jj,jl) = rzero zdvice(ji,jj,jl) = rzero END IF ELSE ! zdaice < - eps10 WRITE(*,*) "zdaice is really negative" zdaice_negative = .TRUE. END IF END IF ! zdaice < 0 ! 2. Repeat for ice volume IF (zdvice(ji,jj,jl) .LT. rzero ) THEN IF (zdvice(ji,jj,jl) .GT. -epsi10 ) THEN WRITE(*,*)," zdvice is negative but artificial" IF ( ( jl1 .EQ. jl .AND. ht_i(ji,jj,jl1) .GT. hi_max(jl) )& &.OR.& &( jl1 .EQ. jl+1 .AND. ht_i(ji,jj,jl1) .LE. hi_max(jl) )) THEN ! You are here if one of these 2 conditions are verified ! 1) The running category is the donor and has too thick a ! sea ice ! 2) The running category is one category below the donor, ! which has too thin a sea ice ! ! If you still with me, it means that a transfer needs to be ! done but the amount of ice is negative due to roundoff ! error or something. Let us reset zdaice and zdvice to the ! value of a_i and v_i in the category for the transfer zdaice(ji,jj,jl) = a_i(ji,jj,jl1) zdvice(ji,jj,jl) = v_i(ji,jj,jl1) ELSE ! None of the two conditions above is verified ! That is, ! 1) The running category is not the donor or has a correct ! ice ! AND ! 2) The category above the running one is not the donor, or ! the ice is correct in this above category ! ! Since everything is fine, nothing should be done zdaice(ji,jj,jl) = rzero zdvice(ji,jj,jl) = rzero END IF ELSE ! zdvice < - eps10 WRITE(*,*) "zdvice is really negative" zdvice_negative = .TRUE. END IF END IF ! zdvice < 0 ! 3.A. If the area to be transferred equals the area in the running ! category, then just update it to the exact value ! (i.e. round it ) IF ( zdaice(ji,jj,jl) .GT. a_i(ji,jj,jl1) - epsi10 .AND.& & zdaice(ji,jj,jl) .LT. a_i(ji,jj,jl1) + epsi10 ) THEN ! if concentration to be transferred is "equal" ! to the concentration of the donor. This is obviously the case ! if the running category is the donor zdaice(ji,jj,jl) = a_i(ji,jj,jl1) zdvice(ji,jj,jl) = v_i(ji,jj,jl1) ELSE ! Otherwise, set the switch to true zdaice_greater_aicen = .TRUE. END IF ! zdaice "=" a_i ! 3.B. (Repeat for volume) ! If the volume to be transferred equals the volume in the running ! category, then keep it as is IF ( zdvice(ji,jj,jl) .GT. v_i(ji,jj,jl1) - epsi10 .AND.& & zdvice(ji,jj,jl) .LT. v_i(ji,jj,jl1) + epsi10 ) THEN ! if volume to be transferred is "equal" ! to the volume of the donor. This is obviously the case ! if the running category is the donor zdaice(ji,jj,jl) = a_i(ji,jj,jl1) zdvice(ji,jj,jl) = v_i(ji,jj,jl1) ELSE zdvice_greater_vicen = .TRUE. END IF ! zdaice "=" a_i ELSE ! if idonor ! Nothing happens since this category is not giving ice END IF ! if idonor END DO ! jj END DO ! ji END DO ! jl, category !------------------------------------------------- ! 2) Transfer volume and energy between categories !------------------------------------------------- DO jl = 1, jpl - 1 DO ji = 1, jpi DO jj = 1, jpj IF (zdaice(ji,jj,jl) .GT. rzero ) THEN jl1 = idonor(ji,jj,jl) ! Define proportionality factor. ! zworka will be the ratio between transferred volume of ice and ! actual volume of ice in the category. Auxiliary quantities (snow volume, snow ! heat content, ice age, salinity, etc.) will be transferred following ! this ratio. Indeed, if out of 4 m of ice, 1 m is transferred, then ! 1/4 of the snow volume will be transferred, too. zindb = MAX( rzero , SIGN( rone , v_i(ji,jj,jl1) - epsi10 ) ) ! Thus zindb is equal to 1 if the donor has positive and not artificially ! close to zero ice volume to give, to zero otherwise zworka(ji,jj) = zdvice(ji,jj,jl) / MAX( v_i(ji,jj,jl1), epsi10 ) *& &zindb ! zworka is zero where the donor has no ice, otherwise ! equal to zdvice/vice of the current category ! We have a donor, but who is the benefiter? who will receive ice? IF ( jl1 == jl ) THEN ! If the donor is the current category, then jl2=jl1+1 ! the receiver is the one above ELSE ! Otherwise (the donor is the above categ) jl2=jl ! then it's the current category! END IF ! ----------------------- ! Go for the transfers!!! ! ----------------------- ! A. Ice areas ! ------------ a_i(ji,jj,jl1) = a_i(ji,jj,jl1) - zdaice(ji,jj,jl) a_i(ji,jj,jl2) = a_i(ji,jj,jl2) + zdaice(ji,jj,jl) ! The donor loses , the receiver gains ! B. Ice volumes ! -------------- v_i(ji,jj,jl1) = v_i(ji,jj,jl1) - zdvice(ji,jj,jl) v_i(ji,jj,jl2) = v_i(ji,jj,jl2) + zdvice(ji,jj,jl) ! C. Snow volumes ! --------------- zdvsnow = v_s(ji,jj,jl1) * zworka(ji,jj) v_s(ji,jj,jl1) = v_s(ji,jj,jl1) - zdvsnow v_s(ji,jj,jl2) = v_s(ji,jj,jl2) + zdvsnow ! D. Snow heat content ! -------------------- zdesnow = e_s(ji,jj,1,jl1) * zworka(ji,jj) e_s(ji,jj,1,jl1) = e_s(ji,jj,1,jl1) - zdesnow e_s(ji,jj,1,jl2) = e_s(ji,jj,1,jl2) + zdesnow ! E. Ice age ! ---------- ! Ice age is defined as areal. If 1/2 of the area ! is transferred then 1/2 of the age too ! Wat een cromme definitie! zdo_aice = oa_i(ji,jj,jl1) * zdaice(ji,jj,jl) oa_i(ji,jj,jl1) = oa_i(ji,jj,jl1) - zdo_aice oa_i(ji,jj,jl2) = oa_i(ji,jj,jl2) + zdo_aice ! F. Ice salinity ! --------------- zdsm_vice = smv_i(ji,jj,jl1) * zworka(ji,jj) smv_i(ji,jj,jl1) = smv_i(ji,jj,jl1) - zdsm_vice smv_i(ji,jj,jl2) = smv_i(ji,jj,jl2) + zdsm_vice ! G. Surface temperature ! ---------------------- zdaTsf = t_su(ji,jj,jl1) * zdaice(ji,jj,jl) zaTsfn(ji,jj,jl1) = zaTsfn(ji,jj,jl1) - zdaTsf zaTsfn(ji,jj,jl2) = zaTsfn(ji,jj,jl2) + zdaTsf ! H. Ice heat content ! ------------------- DO jk = 1 , nlay_i zdeice = e_i(ji,jj,jk,jl1) * zworka(ji,jj) e_i(ji,jj,jk,jl1) = e_i(ji,jj,jk,jl1) - zdeice e_i(ji,jj,jk,jl2) = e_i(ji,jj,jk,jl2) + zdeice END DO ELSE ! WRITE(*,*),"Nothing to transfer" END IF END DO ! jpj END DO ! jpi END DO ! jl, category ! --------------------------------------- ! 3) Update ice thickness and temperature ! --------------------------------------- DO jl = 1 , jpl DO ji = 1 , jpi DO jj = 1 , jpj IF ( a_i(ji,jj,jl) > epsi10 ) THEN ht_i(ji,jj,jl) = v_i(ji,jj,jl) / a_i(ji,jj,jl) t_su(ji,jj,jl) = zaTsfn(ji,jj,jl) / a_i(ji,jj,jl) ELSE ht_i(ji,jj,jl) = rzero t_su(ji,jj,jl) = rtt ! If no ice then set "ice" temperature to ! freezing point END IF END DO ! jj END DO !ji END DO ! jl WRITE(*,*) 'Leaving SUBROUTINE shiftice' END SUBROUTINE shiftice SUBROUTINE conc_check() USE my_variables IMPLICIT NONE ! Dummy variables INTEGER :: ji, jj, jl DO ji=1,jpi DO jj = 1,jpj DO jl=1,jpl ! 1. Individual concentrations greater than zero IF ( a_i(ji,jj,jl) .LT. rzero ) THEN IF (a_i(ji,jj,jl) .LT. -epsi10) THEN ! We have a "true" negative conc WRITE(*,*) 'ERROR: final check: a_i negative at ',ji,jj WRITE(*,*) 'for category ',jl WRITE(*,*) 'Value: ', a_i(ji,jj,jl) WRITE(*,*) 'ABORT' STOP ELSE ! We have a fake negative conc a_i(ji,jj,jl) = rzero ENDIF ENDIF ! 2. Individual concentrations less than one IF ( a_i(ji,jj,jl) .GT. zamax ) THEN IF (a_i(ji,jj,jl) - zamax .GT. epsi10) THEN ! We have a "true" positive conc WRITE(*,*) 'ERROR: final check: a_i more than zamax at ',ji,jj WRITE(*,*) 'for category ',jl WRITE(*,*) 'Value: ', a_i(ji,jj,jl) WRITE(*,*) 'ABORT' STOP ELSE ! We have a fake more than one conc a_i(ji,jj,jl) = zamax ENDIF ENDIF END DO ! jl ! 3. Total concentration greater than zero IF ( SUM(a_i(ji,jj,:)) .LT. rzero ) THEN IF (SUM(a_i(ji,jj,:)) .LT. -epsi10) THEN ! We have a "true" negative total conc WRITE(*,*) 'ERROR: final check: at_i negative at ',ji,jj WRITE(*,*) 'Value: ', SUM(a_i(ji,jj,:)) WRITE(*,*) 'ABORT' STOP ELSE ! We have a fake negative conc ! Let's reset all categories to zero DO jl = 1, jpl a_i(ji,jj,jl)=rzero END DO ENDIF ENDIF ! 4. Total concentration less than one IF ( SUM(a_i(ji,jj,:)) .GT. zamax ) THEN IF (SUM(a_i(ji,jj,:)) - zamax .GT. epsi10) THEN ! We have a "true"positive conc WRITE(*,*) 'ERROR: final check: at_i more than one at ',ji,jj WRITE(*,*) 'Value: ', SUM(a_i(ji,jj,:)) WRITE(*,*) 'Individual: ', a_i(ji,jj,:) WRITE(*,*) 'ABORT' STOP ELSE ! We have a fake more than one conc ! Define the excess, which is by definition negligible zda_i = SUM(a_i(ji,jj,:)) - zamax ! (positive) ! Give this excess to the first category than can accept it, i.e. that ! has less than zamax - zda_i l_adjust = .TRUE. DO jl = 1, jpl IF ( ( a_i(ji,jj,jl) .GT. zda_i ) .AND. l_adjust ) THEN a_i(ji,jj,jl) = a_i(ji,jj,jl) - zda_i l_adjust = .FALSE. END IF ENDDO ! jl IF ( l_adjust ) THEN WRITE(*,*) 'It was not possible to give the excess of ice.' WRITE(*,*) 'At grid point ',ji,jj WRITE(*,*) 'Category', jl WRITE(*,*) 'The excess is ', zda_i WRITE(*,*) 'The conc. in categories are ', a_i(ji,jj,:) WRITE(*,*) ' ABORT' STOP END IF ENDIF ENDIF END DO !jj END DO !ji END SUBROUTINE conc_check SUBROUTINE get_dimensions() USE my_variables USE netcdf IMPLICIT NONE ! Dummy variables WRITE(*,*), 'Entering get_dimensions()' ! 2.A Get dimensions ! ------------------ ! Open the oce restart file ierr = nf90_open(TRIM(cfile_analysis_oce),nf90_Write,incid_oce_an_in) IF (ierr .NE. nf90_noerr ) THEN WRITE(*,*), "(sanity_check_LIM3) Error NetCDF Opening. Abort" WRITE(*,*), "File :" // TRIM(cfile_analysis_oce) STOP END IF ! Inquire variable id (here the variable does not matter!) ierr = nf90_inq_varid(incid_oce_an_in, "sn", ivarid) IF (ierr .NE. nf90_noerr ) THEN WRITE(*,*), "(sanity_check_LIM3) Error NetCDF Variable ID Inquiring. Abort" WRITE(*,*), "File :"//TRIM(cfile_analysis_oce) WRITE(*,*), "Variable: sn" STOP END IF ! Inquire variable ierr = nf90_inquire_variable(incid_oce_an_in, ivarid, dimids=idimid) IF (ierr .NE. nf90_noerr ) THEN WRITE(*,*), "(sanity_check_LIM3) Error NetCDF Variable Inquiring. Abort" WRITE(*,*), "File :"//TRIM(cfile_analysis_oce) WRITE(*,*), "Variable: sn" STOP END IF ! Get dimensions ierr = nf90_inquire_dimension(incid_oce_an_in, idimid(1), len=jpi) ierr = nf90_inquire_dimension(incid_oce_an_in, idimid(2), len=jpj) ierr = nf90_inquire_dimension(incid_oce_an_in, idimid(3), len=jpk) IF (ierr .NE. nf90_noerr ) THEN WRITE(*,*), "(sanity_check_LIM3) Error NetCDF dimensions Inquiring. Abort" WRITE(*,*), "File :"//TRIM(cfile_analysis_oce) STOP END IF WRITE(*,*), "The model horizontal dimensions are" , jpi, "by",jpj WRITE(*,*), "The model vertical dimension is" , jpk WRITE(*,*), 'Leaving get_dimensions()' END SUBROUTINE get_dimensions SUBROUTINE get_nb_cat() USE my_variables USE netcdf IMPLICIT NONE ! Dummy variables INTEGER :: jl WRITE(*,*) 'Entering get_nb_cat()' ! Open file ierr = nf90_open(trim(cfile_analysis_ice),nf90_NoWrite,incid_ice_an_in) IF (ierr .NE. nf90_noerr ) THEN WRITE(*,*), "(sanity_check_LIM3) Error NetCDF Ice file. Abort" WRITE(*,*), TRIM(cfile_analysis_ice) STOP END IF ! Get number of ice categories cvarroot="a_i_htc" jl=1 WRITE(cinteger,"(i1)") jl ierr = nf90_inq_varid(incid_ice_an_in, TRIM(cvarroot)//TRIM(cinteger), ivarid) jl=jl+1 DO WHILE (ierr == nf90_noerr) WRITE(cinteger,"(i1)") jl !WRITE(*,*),"Checking existence of " // TRIM(cconcroot)//TRIM(cinteger) ierr = nf90_inq_varid(incid_ice_an_in, TRIM(cvarroot)//TRIM(cinteger), ivarid) jl=jl+1 ENDDO jpl=jl-2 WRITE(*,*), "There are ", jpl, "ice categories" ! Close ierr = nf90_close(incid_ice_an_in); IF (ierr .NE. nf90_noerr ) THEN WRITE(*,*), "(sanity_check_LIM3) Error NetCDF Ice close file. Abort" WRITE(*,*), TRIM(cfile_analysis_ice) STOP END IF WRITE(*,*) 'Leaving get_nb_cat()' END SUBROUTINE get_nb_cat SUBROUTINE get_nb_il() USE my_variables USE netcdf IMPLICIT NONE ! Dummy variables INTEGER :: jk WRITE(*,*) 'Entering get_nb_il()' ! Open file ierr = nf90_open(trim(cfile_analysis_ice),nf90_NoWrite,incid_ice_an_in) IF (ierr .NE. nf90_noerr ) THEN WRITE(*,*), "(sanity_check_LIM3) Error NetCDF Ice file. Abort" WRITE(*,*), TRIM(cfile_analysis_ice) STOP END IF cvarroot="tempt_il" jk=1 WRITE(cinteger,"(i1)") jk ierr = nf90_inq_varid(incid_ice_an_in, TRIM(cvarroot)//TRIM(cinteger)//'_htc1', ivarid) jk=jk+1 DO WHILE (ierr == nf90_noerr) WRITE(cinteger,"(i1)") jk ierr = nf90_inq_varid(incid_ice_an_in, TRIM(cvarroot)//TRIM(cinteger)//'_htc1' , ivarid) jk=jk+1 ENDDO nlay_i = jk-2 WRITE(*,*), "There are ", nlay_i, "ice layers" ! Close ierr = nf90_close(incid_ice_an_in); IF (ierr .NE. nf90_noerr ) THEN WRITE(*,*), "(sanity_check_LIM3) Error NetCDF Ice close file. Abort" WRITE(*,*), TRIM(cfile_analysis_ice) STOP END IF WRITE(*,*) 'Leaving get_nb_il()' END SUBROUTINE get_nb_il SUBROUTINE get_nb_sl() USE my_variables USE netcdf IMPLICIT NONE ! Dummy variables INTEGER :: jk WRITE(*,*) 'Entering get_nb_sl()' ! Open file ierr = nf90_open(trim(cfile_analysis_ice),nf90_NoWrite,incid_ice_an_in) IF (ierr .NE. nf90_noerr ) THEN WRITE(*,*), "(sanity_check_LIM3) Error NetCDF Ice file. Abort" WRITE(*,*), TRIM(cfile_analysis_ice) STOP END IF cvarroot="tempt_sl" jk=1 WRITE(cinteger,"(i1)") jk ierr = nf90_inq_varid(incid_ice_an_in, TRIM(cvarroot)//TRIM(cinteger)//'_htc1', ivarid) jk=jk+1 DO WHILE (ierr == nf90_noerr) WRITE(cinteger,"(i1)") jk ierr = nf90_inq_varid(incid_ice_an_in, TRIM(cvarroot)//TRIM(cinteger)//'_htc1' ,ivarid) jk=jk+1 ENDDO nlay_s = jk-2 WRITE(*,*), "There are ", nlay_s, "snow layers" ! Close ierr = nf90_close(incid_ice_an_in); IF (ierr .NE. nf90_noerr ) THEN WRITE(*,*), "(sanity_check_LIM3) Error NetCDF Ice close file. Abort" WRITE(*,*), TRIM(cfile_analysis_ice) STOP END IF WRITE(*,*) 'Leaving get_nb_sl()' END SUBROUTINE get_nb_sl SUBROUTINE allocate_variables() USE my_variables IMPLICIT NONE WRITE(*,*), 'Entering allocate_variables()' ALLOCATE( ilandmask(jpi, jpj) ) ALLOCATE( a_i( jpi, jpj, jpl) ,& &v_i( jpi, jpj, jpl) ,& &v_i_fc( jpi, jpj, jpl) ,& &v_s( jpi, jpj, jpl) ,& &v_s_fc( jpi, jpj, jpl) ,& &oa_i ( jpi, jpj, jpl) ,& &smv_i( jpi, jpj, jpl) ,& &t_su( jpi, jpj, jpl) ) ALLOCATE( ht_i( jpi, jpj, jpl) ) ALLOCATE( e_i( jpi, jpj,nlay_i,jpl) ) ALLOCATE( e_s( jpi, jpj,nlay_s,jpl) ) ALLOCATE( sss_m_an( jpi, jpj ) ,& sss_m_fc( jpi, jpj ) ,& sst_m_an( jpi, jpj ) ,& sst_m_fc( jpi, jpj ) ) ALLOCATE( sn_an ( jpi, jpj,jpk ) ,& sn_fc ( jpi, jpj,jpk ) ,& tn_an ( jpi, jpj,jpk ) ,& tn_fc ( jpi, jpj,jpk ) ,& un_an ( jpi, jpj,jpk ) ,& un_fc ( jpi, jpj,jpk ) ,& ub_an ( jpi, jpj,jpk ) ,& ub_fc ( jpi, jpj,jpk ) ,& vn_an ( jpi, jpj,jpk ) ,& vn_fc ( jpi, jpj,jpk ) ,& vb_an ( jpi, jpj,jpk ) ,& vb_fc ( jpi, jpj,jpk ) ) ALLOCATE( ssu_m_an( jpi, jpj ) ,& ssu_m_fc( jpi, jpj ) ,& ssv_m_an( jpi, jpj ) ,& ssv_m_fc( jpi, jpj ) ) ALLOCATE( hi_max( jpl) ) ALLOCATE( idonor( jpi, jpj, jpl) ) ALLOCATE( internal_melt(jpi, jpj, jpl) ) ALLOCATE( zdaice( jpi, jpj, jpl) ,& &zdvice( jpi, jpj, jpl) ) ALLOCATE( ze1t( jpi, jpj ) ,& &ze2t( jpi, jpj ) ,& &zcellarea(jpi, jpj ) ) ALLOCATE(zheat_res( jpi, jpj ) ,& &zdmsnif( jpi, jpj ) ) ALLOCATE(at_i( jpi, jpj ) ,& snwice_mass(jpi,jpj ) ,& snwice_mass_b(jpi,jpj ) ) WRITE(*,*) 'Leaving allocate_variables()' END SUBROUTINE allocate_variables SUBROUTINE load_variables() USE my_variables USE netcdf IMPLICIT NONE ! Dummy variables INTEGER :: jl, jk WRITE(*,*) 'Entering load_variables()' ! 3.A Mask ! -------- ! Open ierr = nf90_open(TRIM(cmaskfile),nf90_NoWrite,incid_mask) IF (ierr .NE. nf90_noerr ) THEN WRITE(*,*), "(sanity_check_LIM3) Error NetCDF Mask file. Abort" WRITE(*,*), TRIM(cmaskfile) STOP END IF ! Request variable ID ierr = nf90_inq_varid(incid_mask, cmaskvar, ivarid) IF (ierr .NE. nf90_noerr ) THEN WRITE(*,*), "(sanity_check_LIM3) Error NetCDF Inquire varID . Abort" WRITE(*,*), TRIM(cmaskfile) STOP END IF ! Get the NetCDF variable and put it in the Fortran variable ierr = nf90_get_var(incid_mask, ivarid, ilandmask) IF (ierr .NE. nf90_noerr ) THEN WRITE(*,*), "(sanity_check_LIM3) Error NetCDF getVar varID . Abort" WRITE(*,*), TRIM(cmaskfile) STOP END IF ! Close mask file ierr = nf90_close(incid_mask) IF (ierr .NE. nf90_noerr ) THEN WRITE(*,*), "(sanity_check_LIM3) Error Closing NetCDF . Abort" WRITE(*,*), TRIM(cmaskfile) STOP END IF ! 3.B Mesh file ! ------------- ! Open ierr = nf90_open(TRIM(cmeshfile),nf90_NoWrite,incid_mesh) IF (ierr .NE. nf90_noerr ) THEN WRITE(*,*), "(sanity_check_LIM3) Error NetCDF Mask file. Abort" WRITE(*,*), TRIM(cmeshfile) STOP END IF WRITE(*,*), "Mesh loaded" ! Request variable ID ierr = nf90_inq_varid(incid_mesh, ce1tvar, ivarid) IF (ierr .NE. nf90_noerr ) THEN WRITE(*,*), "(sanity_check_LIM3) Error NetCDF Inquire varID . Abort" WRITE(*,*), TRIM(cmeshfile) STOP END IF ! Get the NetCDF variable and put it in the Fortran variable ierr = nf90_get_var(incid_mesh, ivarid, ze1t) IF (ierr .NE. nf90_noerr ) THEN WRITE(*,*), "(sanity_check_LIM3) Error NetCDF getVar varID . Abort" WRITE(*,*), TRIM(cmeshfile) STOP END IF ! Repeat for e2t ! Request variable ID ierr = nf90_inq_varid(incid_mesh, ce2tvar, ivarid) IF (ierr .NE. nf90_noerr ) THEN WRITE(*,*), "(sanity_check_LIM3) Error NetCDF Inquire varID . Abort" WRITE(*,*), TRIM(cmeshfile) STOP END IF ! Get the NetCDF variable and put it in the Fortran variable ierr = nf90_get_var(incid_mesh, ivarid, ze2t) IF (ierr .NE. nf90_noerr ) THEN WRITE(*,*), "(sanity_check_LIM3) Error NetCDF getVar varID . Abort" WRITE(*,*), TRIM(cmeshfile) STOP END IF ! Close mesh file ierr = nf90_close(incid_mesh) IF (ierr .NE. nf90_noerr ) THEN WRITE(*,*), "(sanity_check_LIM3) Error Closing NetCDF . Abort" WRITE(*,*), TRIM(cmeshfile) STOP END IF ! Computing zcellarea zcellarea(:,:) = ze1t(:,:) * ze2t(:,:) ! 3.C Ocean variables ! ------------------- ! ! 3.C.1 Open analysis ierr = nf90_open(trim(cfile_analysis_oce),nf90_NoWrite,incid_oce_an_in) IF (ierr .NE. nf90_noerr ) THEN WRITE(*,*), "(sanity_check_LIM3) Error NetCDF Ocean file. Abort" WRITE(*,*), TRIM(cfile_analysis_oce) STOP END IF ! The following lines commented out because SSS does not appear anymore ! as restart variable (2020) !! 3.C.1.A Sea surface salinity !! Request variable ID !ierr = nf90_inq_varid(incid_oce_an_in, csss_m, ivarid) !IF (ierr .NE. nf90_noerr ) THEN ! WRITE(*,*), "(sanity_check_LIM3) Error NetCDF Ocean file var. inquiry. Abort" ! WRITE(*,*), TRIM(cfile_analysis_oce) ! STOP !END IF ! !! Get variable !ierr = nf90_get_var(incid_oce_an_in, ivarid, sss_m_an) !IF (ierr .NE. nf90_noerr ) THEN ! WRITE(*,*), "(sanity_check_LIM3) Error NetCDF Ocean file get. Abort" ! WRITE(*,*), TRIM(cfile_analysis_oce) ! STOP !END IF ! ! 3.C.1.B Salinity ! Request variable ID ierr = nf90_inq_varid(incid_oce_an_in, csn, ivarid) IF (ierr .NE. nf90_noerr ) THEN WRITE(*,*), "(sanity_check_LIM3) Error NetCDF Ocean file var. inquiry. Abort" WRITE(*,*), TRIM(cfile_analysis_oce) STOP END IF ! Get variable ierr = nf90_get_var(incid_oce_an_in, ivarid, sn_an) IF (ierr .NE. nf90_noerr ) THEN WRITE(*,*), "(sanity_check_LIM3) Error NetCDF Ocean file get. Abort" WRITE(*,*), TRIM(cfile_analysis_oce) STOP END IF ! The following lines commented out because SST does not longer ! appear as restart file (2020) !! 3.C.1.C Sea surface temperature !! Request variable ID !ierr = nf90_inq_varid(incid_oce_an_in, csst_m, ivarid) !IF (ierr .NE. nf90_noerr ) THEN ! WRITE(*,*), "(sanity_check_LIM3) Error NetCDF Ocean file var. inquiry. Abort" ! WRITE(*,*), TRIM(cfile_analysis_oce) ! STOP !END IF ! !! Get variable !ierr = nf90_get_var(incid_oce_an_in, ivarid, sst_m_an) !IF (ierr .NE. nf90_noerr ) THEN ! WRITE(*,*), "(sanity_check_LIM3) Error NetCDF Ocean file get. Abort" ! WRITE(*,*), TRIM(cfile_analysis_oce) ! STOP !END IF ! 3.C.1.D. Temperature ierr = nf90_inq_varid(incid_oce_an_in, ctn, ivarid) IF (ierr .NE. nf90_noerr ) THEN WRITE(*,*), "(sanity_check_LIM3) Error NetCDF Ocean file var. inquiry. Abort" WRITE(*,*), TRIM(cfile_analysis_oce) STOP END IF ! Get variable ierr = nf90_get_var(incid_oce_an_in, ivarid, tn_an) IF (ierr .NE. nf90_noerr ) THEN WRITE(*,*), "(sanity_check_LIM3) Error NetCDF Ocean file get. Abort" WRITE(*,*), TRIM(cfile_analysis_oce) STOP END IF ! 3.C.1.E U-velocity ("un") ierr = nf90_inq_varid(incid_oce_an_in, cun, ivarid) IF (ierr .NE. nf90_noerr ) THEN WRITE(*,*), "(sanity_check_LIM3) Error NetCDF Ocean file var. inquiry. Abort" WRITE(*,*), TRIM(cfile_analysis_oce) STOP END IF ! Get variable ierr = nf90_get_var(incid_oce_an_in, ivarid, un_an) IF (ierr .NE. nf90_noerr ) THEN WRITE(*,*), "(sanity_check_LIM3) Error NetCDF Ocean file get. Abort" WRITE(*,*), TRIM(cfile_analysis_oce) STOP END IF ! 3.C.1.F U-velocity ("ub") ierr = nf90_inq_varid(incid_oce_an_in, cub, ivarid) IF (ierr .NE. nf90_noerr ) THEN WRITE(*,*), "(sanity_check_LIM3) Error NetCDF Ocean file var. inquiry. Abort" WRITE(*,*), TRIM(cfile_analysis_oce) STOP END IF ! Get variable ierr = nf90_get_var(incid_oce_an_in, ivarid, ub_an) IF (ierr .NE. nf90_noerr ) THEN WRITE(*,*), "(sanity_check_LIM3) Error NetCDF Ocean file get. Abort" WRITE(*,*), TRIM(cfile_analysis_oce) STOP END IF ! 3.C.1.G V-velocity ("vn") ierr = nf90_inq_varid(incid_oce_an_in, cvn, ivarid) IF (ierr .NE. nf90_noerr ) THEN WRITE(*,*), "(sanity_check_LIM3) Error NetCDF Ocean file var. inquiry. Abort" WRITE(*,*), TRIM(cfile_analysis_oce) STOP END IF ! Get variable ierr = nf90_get_var(incid_oce_an_in, ivarid, vn_an) IF (ierr .NE. nf90_noerr ) THEN WRITE(*,*), "(sanity_check_LIM3) Error NetCDF Ocean file get. Abort" WRITE(*,*), TRIM(cfile_analysis_oce) STOP END IF ! 3.C.1.H V-velocity ("vb") ierr = nf90_inq_varid(incid_oce_an_in, cvb, ivarid) IF (ierr .NE. nf90_noerr ) THEN WRITE(*,*), "(sanity_check_LIM3) Error NetCDF Ocean file var. inquiry. Abort" WRITE(*,*), TRIM(cfile_analysis_oce) STOP END IF ! Get variable ierr = nf90_get_var(incid_oce_an_in, ivarid, vb_an) IF (ierr .NE. nf90_noerr ) THEN WRITE(*,*), "(sanity_check_LIM3) Error NetCDF Ocean file get. Abort" WRITE(*,*), TRIM(cfile_analysis_oce) STOP END IF ! The following lines commented out because surface velocities no longer ! appear as restart variables (2020) !! 3.C.1.I SSU-velocity ("ssu_m") !ierr = nf90_inq_varid(incid_oce_an_in, cssu_m, ivarid) !IF (ierr .NE. nf90_noerr ) THEN ! WRITE(*,*), "(sanity_check_LIM3) Error NetCDF Ocean file var. inquiry. Abort" ! WRITE(*,*), TRIM(cfile_analysis_oce) ! STOP !END IF ! !! Get variable !ierr = nf90_get_var(incid_oce_an_in, ivarid, ssu_m_an) !IF (ierr .NE. nf90_noerr ) THEN ! WRITE(*,*), "(sanity_check_LIM3) Error NetCDF Ocean file get. Abort" ! WRITE(*,*), TRIM(cfile_analysis_oce) ! STOP !END IF ! !! 3.C.1.J SSV-velocity ("ssv_m") !ierr = nf90_inq_varid(incid_oce_an_in, cssv_m, ivarid) !IF (ierr .NE. nf90_noerr ) THEN ! WRITE(*,*), "(sanity_check_LIM3) Error NetCDF Ocean file var. inquiry. Abort" ! WRITE(*,*), TRIM(cfile_analysis_oce) ! STOP !END IF ! !! Get variable !ierr = nf90_get_var(incid_oce_an_in, ivarid, ssv_m_an) !IF (ierr .NE. nf90_noerr ) THEN ! WRITE(*,*), "(sanity_check_LIM3) Error NetCDF Ocean file get. Abort" ! WRITE(*,*), TRIM(cfile_analysis_oce) ! STOP !END IF ! Close analysis ierr = nf90_close(incid_oce_an_in); IF (ierr .NE. nf90_noerr ) THEN WRITE(*,*), "(sanity_check_LIM3) Error NetCDF Ocean close file. Abort" WRITE(*,*), TRIM(cfile_analysis_oce) STOP END IF ! 3.C.2. Open forecast ierr = nf90_open(trim(cfile_forecast_oce),nf90_NoWrite,incid_oce_fc_in) IF (ierr .NE. nf90_noerr ) THEN WRITE(*,*), "(sanity_check_LIM3) Error NetCDF Ocean file. Abort" WRITE(*,*), TRIM(cfile_forecast_oce) STOP END IF ! Following lines commented out as the variable no longer appears ! in the restart files (2020 !! 3.C.2.A Sea surface salinity !! Request variable ID !ierr = nf90_inq_varid(incid_oce_fc_in, csss_m, ivarid) !IF (ierr .NE. nf90_noerr ) THEN ! WRITE(*,*), "(sanity_check_LIM3) Error NetCDF Ocean file var. inquiry. Abort" ! WRITE(*,*), TRIM(cfile_forecast_oce) ! STOP !END IF ! !! Get variable !ierr = nf90_get_var(incid_oce_fc_in, ivarid, sss_m_fc) !IF (ierr .NE. nf90_noerr ) THEN ! WRITE(*,*), "(sanity_check_LIM3) Error NetCDF Ocean file get. Abort" ! WRITE(*,*), TRIM(cfile_forecast_oce) ! STOP !END IF ! 3.C.2.B. Salinity ! Request variable ID ierr = nf90_inq_varid(incid_oce_fc_in, csn, ivarid) IF (ierr .NE. nf90_noerr ) THEN WRITE(*,*), "(sanity_check_LIM3) Error NetCDF Ocean file var. inquiry. Abort" WRITE(*,*), TRIM(cfile_forecast_oce) STOP END IF ! Get variable ierr = nf90_get_var(incid_oce_fc_in, ivarid, sn_fc) IF (ierr .NE. nf90_noerr ) THEN WRITE(*,*), "(sanity_check_LIM3) Error NetCDF Ocean file get. Abort" WRITE(*,*), TRIM(cfile_forecast_oce) STOP END IF ! Following lines commented out as the variable no longer appears ! in the restart files (2020) !! 3.C.2.C Sea surface temperature !! Request variable ID !ierr = nf90_inq_varid(incid_oce_fc_in, csst_m, ivarid) !IF (ierr .NE. nf90_noerr ) THEN ! WRITE(*,*), "(sanity_check_LIM3) Error NetCDF Ocean file var. inquiry. Abort" ! WRITE(*,*), TRIM(cfile_forecast_oce) ! STOP !END IF ! !! Get variable !ierr = nf90_get_var(incid_oce_fc_in, ivarid, sst_m_fc) !IF (ierr .NE. nf90_noerr ) THEN ! WRITE(*,*), "(sanity_check_LIM3) Error NetCDF Ocean file get. Abort" ! WRITE(*,*), TRIM(cfile_forecast_oce) ! STOP !END IF ! 3.C.2.D. Temperature ! Request variable ID ierr = nf90_inq_varid(incid_oce_fc_in, ctn, ivarid) IF (ierr .NE. nf90_noerr ) THEN WRITE(*,*), "(sanity_check_LIM3) Error NetCDF Ocean file var. inquiry. Abort" WRITE(*,*), TRIM(cfile_forecast_oce) STOP END IF ! Get variable ierr = nf90_get_var(incid_oce_fc_in, ivarid, tn_fc) IF (ierr .NE. nf90_noerr ) THEN WRITE(*,*), "(sanity_check_LIM3) Error NetCDF Ocean file get. Abort" WRITE(*,*), TRIM(cfile_forecast_oce) STOP END IF ! 3.C.2.E U-velocity ("un") ierr = nf90_inq_varid(incid_oce_fc_in, cun, ivarid) IF (ierr .NE. nf90_noerr ) THEN WRITE(*,*), "(sanity_check_LIM3) Error NetCDF Ocean file var. inquiry. Abort" WRITE(*,*), TRIM(cfile_forecast_oce) STOP END IF ! Get variable ierr = nf90_get_var(incid_oce_fc_in, ivarid, un_fc) IF (ierr .NE. nf90_noerr ) THEN WRITE(*,*), "(sanity_check_LIM3) Error NetCDF Ocean file get. Abort" WRITE(*,*), TRIM(cfile_forecast_oce) STOP END IF ! 3.C.2.F U-velocity ("ub") ierr = nf90_inq_varid(incid_oce_fc_in, cub, ivarid) IF (ierr .NE. nf90_noerr ) THEN WRITE(*,*), "(sanity_check_LIM3) Error NetCDF Ocean file var. inquiry. Abort" WRITE(*,*), TRIM(cfile_forecast_oce) STOP END IF ! Get variable ierr = nf90_get_var(incid_oce_fc_in, ivarid, ub_fc) IF (ierr .NE. nf90_noerr ) THEN WRITE(*,*), "(sanity_check_LIM3) Error NetCDF Ocean file get. Abort" WRITE(*,*), TRIM(cfile_forecast_oce) STOP END IF ! 3.C.2.G V-velocity ("vn") ierr = nf90_inq_varid(incid_oce_fc_in, cvn, ivarid) IF (ierr .NE. nf90_noerr ) THEN WRITE(*,*), "(sanity_check_LIM3) Error NetCDF Ocean file var. inquiry. Abort" WRITE(*,*), TRIM(cfile_forecast_oce) STOP END IF ! Get variable ierr = nf90_get_var(incid_oce_fc_in, ivarid, vn_fc) IF (ierr .NE. nf90_noerr ) THEN WRITE(*,*), "(sanity_check_LIM3) Error NetCDF Ocean file get. Abort" WRITE(*,*), TRIM(cfile_forecast_oce) STOP END IF ! 3.C.2.H V-velocity ("vb") ierr = nf90_inq_varid(incid_oce_fc_in, cvb, ivarid) IF (ierr .NE. nf90_noerr ) THEN WRITE(*,*), "(sanity_check_LIM3) Error NetCDF Ocean file var. inquiry. Abort" WRITE(*,*), TRIM(cfile_forecast_oce) STOP END IF ! Get variable ierr = nf90_get_var(incid_oce_fc_in, ivarid, vb_fc) IF (ierr .NE. nf90_noerr ) THEN WRITE(*,*), "(sanity_check_LIM3) Error NetCDF Ocean file get. Abort" WRITE(*,*), TRIM(cfile_forecast_oce) STOP END IF ! Following lines commented out as the variable no longer appears ! in the restart files (2020) !! 3.C.2.I SSU-velocity ("ssu_m") !ierr = nf90_inq_varid(incid_oce_fc_in, cssu_m, ivarid) !IF (ierr .NE. nf90_noerr ) THEN ! WRITE(*,*), "(sanity_check_LIM3) Error NetCDF Ocean file var. inquiry. Abort" ! WRITE(*,*), TRIM(cfile_forecast_oce) ! STOP !END IF ! !! Get variable !ierr = nf90_get_var(incid_oce_fc_in, ivarid, ssu_m_fc) !IF (ierr .NE. nf90_noerr ) THEN ! WRITE(*,*), "(sanity_check_LIM3) Error NetCDF Ocean file get. Abort" ! WRITE(*,*), TRIM(cfile_forecast_oce) ! STOP !END IF !! 3.C.2.J SSV-velocity ("ssv_m") !ierr = nf90_inq_varid(incid_oce_fc_in, cssv_m, ivarid) !IF (ierr .NE. nf90_noerr ) THEN ! WRITE(*,*), "(sanity_check_LIM3) Error NetCDF Ocean file var. inquiry. Abort" ! WRITE(*,*), TRIM(cfile_forecast_oce) ! STOP !END IF ! ! Get variable !ierr = nf90_get_var(incid_oce_fc_in, ivarid, ssv_m_fc) !IF (ierr .NE. nf90_noerr ) THEN ! WRITE(*,*), "(sanity_check_LIM3) Error NetCDF Ocean file get. Abort" ! WRITE(*,*), TRIM(cfile_forecast_oce) ! STOP !END IF ! Close forecast ierr = nf90_close(incid_oce_fc_in); IF (ierr .NE. nf90_noerr ) THEN WRITE(*,*), "(sanity_check_LIM3) Error NetCDF Ocean close file. Abort" WRITE(*,*), TRIM(cfile_forecast_oce) STOP END IF WRITE(*,*), "Ocean variables loaded" ! 3.D Ice variables ! ----------------- ! Open forecast ierr = nf90_open(trim(cfile_forecast_ice),nf90_NoWrite,incid_ice_fc_in) IF (ierr .NE. nf90_noerr ) THEN WRITE(*,*), "(sanity_check_LIM3) Error NetCDF Ice file. Abort" WRITE(*,*), TRIM(cfile_forecast_ice) STOP END IF DO jl = 1 , jpl ! Read ice volume forecast cvarroot='v_i_htc' WRITE(cinteger,'(i1)') jl cvarname = TRIM(cvarroot) // TRIM(cinteger) ! Inquire variable ID ierr = nf90_inq_varid(incid_ice_fc_in, cvarname, ivarid) IF (ierr .NE. nf90_noerr ) THEN WRITE(*,*), "(sanity_check_LIM3) Error NetCDF Ice file var. inquiry. Abort" WRITE(*,*), TRIM(cfile_forecast_ice) WRITE(*,*), cvarname STOP END IF ! Read variable ierr = nf90_get_var(incid_ice_fc_in, ivarid, v_i_fc(:,:,jl) ) IF (ierr .NE. nf90_noerr ) THEN WRITE(*,*), "(sanity_check_LIM3) Error NetCDF Ice file get. Abort" WRITE(*,*), TRIM(cfile_forecast_ice) WRITE(*,*), cvarname STOP END IF ! Read snow volume forecast cvarroot='v_s_htc' WRITE(cinteger,'(i1)') jl cvarname = TRIM(cvarroot) // TRIM(cinteger) ! Inquire variable ID ierr = nf90_inq_varid(incid_ice_fc_in, cvarname, ivarid) IF (ierr .NE. nf90_noerr ) THEN WRITE(*,*), "(sanity_check_LIM3) Error NetCDF Ice file var. inquiry. Abort" WRITE(*,*), TRIM(cfile_forecast_ice) WRITE(*,*), cvarname STOP END IF ! Read variable ierr = nf90_get_var(incid_ice_fc_in, ivarid, v_s_fc(:,:,jl) ) IF (ierr .NE. nf90_noerr ) THEN WRITE(*,*), "(sanity_check_LIM3) Error NetCDF Ice file get. Abort" WRITE(*,*), TRIM(cfile_forecast_ice) WRITE(*,*), cvarname STOP END IF END DO ! jl ! Close forecast ierr = nf90_close(incid_ice_fc_in); IF (ierr .NE. nf90_noerr ) THEN WRITE(*,*), "(sanity_check_LIM3) Error NetCDF Ice close file. Abort" WRITE(*,*), TRIM(cfile_forecast_ice) STOP END IF ! Open analysis ierr = nf90_open(trim(cfile_analysis_ice),nf90_NoWrite,incid_ice_an_in) IF (ierr .NE. nf90_noerr ) THEN WRITE(*,*), "(sanity_check_LIM3) Error NetCDF Ice file. Abort" WRITE(*,*), TRIM(cfile_analysis_ice) STOP END IF ! Time counter ! Request variable id cvarname = 'time_counter' ierr = nf90_inq_varid(incid_ice_an_in, cvarname, ivarid) IF (ierr .NE. nf90_noerr ) THEN WRITE(*,*), "(sanity_check_LIM3) Error NetCDF Ice file var. inquiry. Abort" WRITE(*,*), TRIM(cfile_analysis_ice) WRITE(*,*), cvarname STOP END IF ! Get variable ierr = nf90_get_var(incid_ice_an_in, ivarid, ztime_counter ) IF (ierr .NE. nf90_noerr ) THEN WRITE(*,*), "(sanity_check_LIM3) Error NetCDF Ice file get. Abort" WRITE(*,*), TRIM(cfile_analysis_ice) WRITE(*,*), cvarname STOP END IF ! Strategy: Loop over categories, and for specific variables over layers DO jl = 1 , jpl ! 3.D.1. Ice concentration cvarroot='a_i_htc' WRITE(cinteger,'(i1)') jl cvarname = TRIM(cvarroot) // TRIM(cinteger) ! WRITE(*,*), "Working with variable " // cvarname ! Request variable ID ierr = nf90_inq_varid(incid_ice_an_in, cvarname, ivarid) IF (ierr .NE. nf90_noerr ) THEN WRITE(*,*), "(sanity_check_LIM3) Error NetCDF Ice file var. inquiry. Abort" WRITE(*,*), TRIM(cfile_analysis_ice) WRITE(*,*), cvarname STOP END IF ! Get variable ierr = nf90_get_var(incid_ice_an_in, ivarid, a_i(:,:,jl) ) IF (ierr .NE. nf90_noerr ) THEN WRITE(*,*), "(sanity_check_LIM3) Error NetCDF Ice file get. Abort" WRITE(*,*), TRIM(cfile_analysis_ice) WRITE(*,*), cvarname STOP END IF ! 3.D.2. Ice volume per surface area cvarroot='v_i_htc' WRITE(cinteger,'(i1)') jl cvarname = TRIM(cvarroot) // TRIM(cinteger) ! WRITE(*,*), "Working with variable " // cvarname ! Request variable ID ierr = nf90_inq_varid(incid_ice_an_in, cvarname, ivarid) IF (ierr .NE. nf90_noerr ) THEN WRITE(*,*), "(sanity_check_LIM3) Error NetCDF Ice file var. inquiry. Abort" WRITE(*,*), TRIM(cfile_analysis_ice) WRITE(*,*), cvarname STOP END IF ! Get variable ierr = nf90_get_var(incid_ice_an_in, ivarid, v_i(:,:,jl) ) IF (ierr .NE. nf90_noerr ) THEN WRITE(*,*), "(sanity_check_LIM3) Error NetCDF Ice file get. Abort" WRITE(*,*), TRIM(cfile_analysis_ice) WRITE(*,*), cvarname STOP END IF ! 3.D.3. Snow volume per surface area cvarroot='v_s_htc' WRITE(cinteger,'(i1)') jl cvarname = TRIM(cvarroot) // TRIM(cinteger) ! WRITE(*,*), "Working with variable " // cvarname ! Request variable ID ierr = nf90_inq_varid(incid_ice_an_in, cvarname, ivarid) IF (ierr .NE. nf90_noerr ) THEN WRITE(*,*), "(sanity_check_LIM3) Error NetCDF Ice file var. inquiry. Abort" WRITE(*,*), TRIM(cfile_analysis_ice) WRITE(*,*), cvarname STOP END IF ! Get variable ierr = nf90_get_var(incid_ice_an_in, ivarid, v_s(:,:,jl) ) IF (ierr .NE. nf90_noerr ) THEN WRITE(*,*), "(sanity_check_LIM3) Error NetCDF Ice file get. Abort" WRITE(*,*), TRIM(cfile_analysis_ice) WRITE(*,*), cvarname STOP END IF ! 3.D.4. Ice age cvarroot='oa_i_htc' WRITE(cinteger,'(i1)') jl cvarname = TRIM(cvarroot) // TRIM(cinteger) ! WRITE(*,*), "Working with variable " // cvarname ! Request variable ID ierr = nf90_inq_varid(incid_ice_an_in, cvarname, ivarid) IF (ierr .NE. nf90_noerr ) THEN WRITE(*,*), "(sanity_check_LIM3) Error NetCDF Ice file var. inquiry. Abort" WRITE(*,*), TRIM(cfile_analysis_ice) WRITE(*,*), cvarname STOP END IF ! Get variable ierr = nf90_get_var(incid_ice_an_in, ivarid, oa_i(:,:,jl) ) IF (ierr .NE. nf90_noerr ) THEN WRITE(*,*), "(sanity_check_LIM3) Error NetCDF Ice file get. Abort" WRITE(*,*), TRIM(cfile_analysis_ice) WRITE(*,*), cvarname STOP END IF ! 3.D.5. Ice salinity cvarroot='smv_i_htc' WRITE(cinteger,'(i1)') jl cvarname = TRIM(cvarroot) // TRIM(cinteger) ! WRITE(*,*), "Working with variable " // cvarname ! Request variable ID ierr = nf90_inq_varid(incid_ice_an_in, cvarname, ivarid) IF (ierr .NE. nf90_noerr ) THEN WRITE(*,*), "(sanity_check_LIM3) Error NetCDF Ice file var. inquiry. Abort" WRITE(*,*), TRIM(cfile_analysis_ice) WRITE(*,*), cvarname STOP END IF ! Get variable ierr = nf90_get_var(incid_ice_an_in, ivarid, smv_i(:,:,jl) ) IF (ierr .NE. nf90_noerr ) THEN WRITE(*,*), "(sanity_check_LIM3) Error NetCDF Ice file get. Abort" WRITE(*,*), TRIM(cfile_analysis_ice) WRITE(*,*), cvarname STOP END IF ! 3.D.5. Ice surface temperature cvarroot='t_su_htc' WRITE(cinteger,'(i1)') jl cvarname = TRIM(cvarroot) // TRIM(cinteger) ! WRITE(*,*), "Working with variable " // cvarname ! Request variable ID ierr = nf90_inq_varid(incid_ice_an_in, cvarname, ivarid) IF (ierr .NE. nf90_noerr ) THEN WRITE(*,*), "(sanity_check_LIM3) Error NetCDF Ice file var. inquiry. Abort" WRITE(*,*), TRIM(cfile_analysis_ice) WRITE(*,*), cvarname STOP END IF ! Get variable ierr = nf90_get_var(incid_ice_an_in, ivarid, t_su(:,:,jl) ) IF (ierr .NE. nf90_noerr ) THEN WRITE(*,*), "(sanity_check_LIM3) Error NetCDF Ice file get. Abort" WRITE(*,*), TRIM(cfile_analysis_ice) WRITE(*,*), cvarname STOP END IF ! 3.D.X Variables that are defined on several ice layers DO jk = 1 , nlay_i ! Ice enthalpy in layers cvarroot='tempt_il' WRITE(cinteger2,'(i1)') jk cvarname = TRIM(cvarroot) // TRIM(cinteger2)// '_htc'//TRIM(cinteger) ! WRITE(*,*),"Working with variable " // cvarname ! Request variable ID ierr = nf90_inq_varid(incid_ice_an_in, cvarname, ivarid) IF (ierr .NE. nf90_noerr ) THEN WRITE(*,*), "(sanity_check_LIM3) Error NetCDF Ice file var. inquiry. Abort" WRITE(*,*), TRIM(cfile_analysis_ice) WRITE(*,*), cvarname STOP END IF ! Get variable ierr = nf90_get_var(incid_ice_an_in, ivarid, e_i(:,:,jk,jl) ) IF (ierr .NE. nf90_noerr ) THEN WRITE(*,*), "(sanity_check_LIM3) Error NetCDF Ice file get. Abort" WRITE(*,*), TRIM(cfile_analysis_ice) WRITE(*,*), cvarname STOP END IF ENDDO ! jk, layers DO jk = 1 , nlay_s ! Snow enthalpy in layers cvarroot='tempt_sl' WRITE(cinteger2,'(i1)') jk cvarname = TRIM(cvarroot) // TRIM(cinteger2)// '_htc'//TRIM(cinteger) ! WRITE(*,*),"Working with variable " // cvarname ! Request variable ID ierr = nf90_inq_varid(incid_ice_an_in, cvarname, ivarid) IF (ierr .NE. nf90_noerr ) THEN WRITE(*,*), "(sanity_check_LIM3) Error NetCDF Ice file var. inquiry. Abort" WRITE(*,*), TRIM(cfile_analysis_ice) WRITE(*,*), cvarname STOP END IF ! Get variable ierr = nf90_get_var(incid_ice_an_in, ivarid, e_s(:,:,jk,jl) ) IF (ierr .NE. nf90_noerr ) THEN WRITE(*,*), "(sanity_check_LIM3) Error NetCDF Ice file get. Abort" WRITE(*,*), TRIM(cfile_analysis_ice) WRITE(*,*), cvarname STOP END IF ENDDO ! jk, layers ENDDO ! jl, categories ! Load data that don't depend on categories ! Snow ice mass cvarname = 'snwice_mass' ! Request variable ID ierr = nf90_inq_varid(incid_ice_an_in, cvarname, ivarid) IF (ierr .NE. nf90_noerr ) THEN WRITE(*,*), "(sanity_check_LIM3) Error NetCDF Ice file var. inquiry. Abort" WRITE(*,*), TRIM(cfile_analysis_ice) WRITE(*,*), cvarname STOP END IF ! Get variable ierr = nf90_get_var(incid_ice_an_in, ivarid, snwice_mass(:,:) ) IF (ierr .NE. nf90_noerr ) THEN WRITE(*,*), "(sanity_check_LIM3) Error NetCDF Ice file get. Abort" WRITE(*,*), TRIM(cfile_analysis_ice) WRITE(*,*), cvarname STOP END IF ! Snow ice mass_b cvarname = 'snwice_mass_b' ! Request variable ID ierr = nf90_inq_varid(incid_ice_an_in, cvarname, ivarid) IF (ierr .NE. nf90_noerr ) THEN WRITE(*,*), "(sanity_check_LIM3) Error NetCDF Ice file var. inquiry. Abort" WRITE(*,*), TRIM(cfile_analysis_ice) WRITE(*,*), cvarname STOP END IF ! Get variable ierr = nf90_get_var(incid_ice_an_in, ivarid, snwice_mass_b(:,:) ) IF (ierr .NE. nf90_noerr ) THEN WRITE(*,*), "(sanity_check_LIM3) Error NetCDF Ice file get. Abort" WRITE(*,*), TRIM(cfile_analysis_ice) WRITE(*,*), cvarname STOP END IF ! Close ice analysis file ierr = nf90_close(incid_ice_an_in); IF (ierr .NE. nf90_noerr ) THEN WRITE(*,*), "(sanity_check_LIM3) Error NetCDF Ice close file. Abort" WRITE(*,*), TRIM(cfile_analysis_ice) STOP END IF WRITE(*,*) 'Leaving load_variables' END SUBROUTINE load_variables SUBROUTINE several_ice_corrections() USE my_variables IMPLICIT NONE ! Dummy variables INTEGER :: ji, jj, jl, jk WRITE(*,*), 'Entering several_ice_corrections()' DO ji = 1 , jpi DO jj = 1 , jpj DO jl = 1 , jpl ! Define switches (masks) zindb = MAX( rzero , SIGN( zamax , a_i(ji,jj,jl) - epsi06) ) zindsn = MAX(rzero , SIGN( zamax , v_s(ji,jj,jl) - epsi06) ) zindic = MAX(rzero , SIGN( zamax , v_i(ji,jj,jl) - epsi04) ) zindb = zindb * zindic ! Mask where there is conc. AND volume ! A. Corrections to ice age IF ( ( oa_i(ji,jj,jl) - rone ) * 86400.0 .GT. ( rdt_ice*ztime_counter*a_i(ji,jj,jl) ) )& &THEN !WRITE(*,*) 'Correction on ice age' oa_i(ji,jj,jl) = rdt_ice * ztime_counter / 86400.0 * a_i(ji,jj,jl) END IF oa_i(ji,jj,jl) = zindb*oa_i(ji,jj,jl) ! B. Corrections to snow thickness v_s(ji,jj,jl) = zindsn*zindb*v_s(ji,jj,jl) ! C. Corrections to ice thickness v_i(ji,jj,jl) = MAX( zindb * v_i(ji,jj,jl) , rzero ) v_i(ji,jj,jl) = zindic*v_i(ji,jj,jl) ! D. Snow transformed to ice if original ice cover disappears zindg = REAL(ilandmask(ji,jj) ) * MAX (rzero, SIGN( rone , - v_i (ji,jj,jl) ) ) ! (is it possible to have zindg not zero??) v_i(ji,jj,jl) = v_i(ji,jj,jl) + zindg * zrhosn * v_s(ji,jj,jl) / zrho0 ! The last term of the above equation is the water volume equivalent to the snow ! volume I guess v_s(ji,jj,jl) = (rone - zindg ) * v_s(ji,jj,jl) + & & zindg * v_i(ji,jj,jl) * (zrho0 - zrhoic ) / zrhosn ! E. Correction to ice concentration a_i(ji,jj,jl) = zindb * MAX(zindsn, zindic) * MAX(a_i(ji,jj,jl), epsi06) ! F. Correction to snow heat content IF ( ldebug .AND. ( ji == jiindx ) .AND. ( jj == jjindx ) ) THEN WRITE(*,*) 'Before:' , e_s(ji,jj,1,jl) END IF e_s(ji,jj,1,jl) = zindsn * & & ( MIN ( MAX ( rzero, e_s(ji,jj,1,jl) ), zbigvalue ) ) + & & ( rone - zindsn ) * rzero IF ( ldebug .AND. ( ji == jiindx ) .AND. ( jj == jjindx ) ) THEN WRITE(*,*) 'After:' , e_s(ji,jj,1,jl) END IF ! G. Correction to ice heat content DO jk = 1 , nlay_i e_i(ji,jj,jk,jl) = zindic * & &( MIN( MAX( rzero, e_i(ji,jj,jk,jl) ), zbigvalue ) ) & & + (rone - zindic ) * rzero END DO ! nlay_i ! H. Correction to ice salinity IF ( (num_sal .EQ. 2) .OR. (num_sal .EQ. 4) ) THEN ! If we are in salinity profile mode ! Salinity stays in bounds smv_i(ji,jj,jl) = MAX( MIN( (zrhoic-zrhosn)/zrhoic * & & sss_m_an(ji,jj) , smv_i(ji,jj,jl) ) , & 0.1 * v_i(ji,jj,jl) ) ENDIF ! I. Thickness of first category above 5cm IF ( jl == 1) THEN ht_i(ji,jj,jl) = zindb * v_i(ji,jj,jl) / MAX(a_i(ji,jj,jl),epsi06) zh = MAX(rone , zindb*zhiclim/ & & MAX(ht_i(ji,jj,jl),epsi20 ) ) ! This is a correction factor equal to zhiclim/h_insitu, i.e. the ratio ! between minimal and actual in situ thickness. ! Something to do for the snow ! The ice concentration is shrunk so that the ice thickness is at least ! zhiclim a_i(ji,jj,jl) = a_i(ji,jj,jl) / zh END IF END DO ! jl ! Reset snowice to zero if necessary IF ( (snwice_mass(ji,jj) .LT. rzero) .OR. (snwice_mass_b(ji,jj) .LT. rzero) ) THEN snwice_mass(ji,jj) = rzero snwice_mass_b(ji,jj)=rzero END IF END DO ! jj END DO ! ji WRITE(*,*), 'Leaving several_ice_corrections()' END SUBROUTINE several_ice_corrections SUBROUTINE shrink_ice() USE my_variables IMPLICIT NONE ! Dummy variables INTEGER :: ji, jj, jl WRITE(*,*), 'Entering shrink_ice()' ! Total concentration cannot exceed zamax at_i(:,:) = rzero DO jl = 1 , jpl at_i(:,:) = a_i(:,:,jl) + at_i(:,:) END DO DO ji = 1 , jpi DO jj = 1 , jpj ! Define the excessive concentration zda_ex = MAX( at_i(ji,jj) - zamax , rzero ) ! (i.e. rzero, except if excess) DO jl = 1 , jpl zindb = MAX( rzero , SIGN( rone , v_i(ji,jj,jl) ) ) ! (zindb is a mask) ! Spread the excess over the different categories with weight equal ! to concentration in each of them zda_i = a_i(ji,jj,jl) * zda_ex / MAX(at_i(ji,jj),epsi06) * zindb ! Correction of limupdate after Antoine Barthélemy ! Indeed the volume should not be changed. !! zdv_i = v_i(ji,jj,jl) * zda_i / MAX(at_i(ji,jj),epsi06) ! We take out the excess of ice and put it as volume a_i(ji,jj,jl) = a_i(ji,jj,jl) - zda_i !! v_i(ji,jj,jl) = v_i(ji,jj,jl) + zdv_i END DO END DO END DO WRITE(*,*), 'Leaving shrink_ice()' END SUBROUTINE shrink_ice SUBROUTINE record_ice() USE my_variables USE netcdf IMPLICIT NONE ! Dummy variables INTEGER :: jl, jk WRITE(*,*) 'Entering record_ice()' WRITE(*,*) 'Recording the NetCDF' ! 8.1 Record ice variables ! We copy the input file and store everything into this copy CALL system('cp -f '//trim(cfile_analysis_ice)//' '//trim(cfileout_ice) ) ierr = nf90_open(trim(cfileout_ice),nf90_Write,incid_ice_an_out) IF (ierr .NE. nf90_noerr ) THEN WRITE(*,*), "(sanity_check_LIM3) Error NetCDF Ocean file. Abort" WRITE(*,*), TRIM(cfileout_ice) STOP END IF ! Loop over categories DO jl = 1 , jpl ! 8.1.1 a_i !~~~~~~~~ cvarroot='a_i_htc' WRITE(cinteger,'(i1)') jl cvarname = TRIM(cvarroot) // TRIM(cinteger) ! WRITE(*,*), "Recording variable " // cvarname ! Request variable ID ierr = nf90_inq_varid(incid_ice_an_out, cvarname, ivarid) IF (ierr .NE. nf90_noerr ) THEN WRITE(*,*), "(sanity_check_LIM3) Error NetCDF ice file var. inquiry. Abort" WRITE(*,*), TRIM(cfileout_ice) STOP END IF ! Put variable ierr = nf90_put_var(incid_ice_an_out, ivarid, a_i(:,:,jl)) IF (ierr .NE. nf90_noerr ) THEN WRITE(*,*), "(sanity_check_LIM3) Error NetCDF ice file var. put. Abort" WRITE(*,*), TRIM(cfileout_ice) STOP END IF ! 8.1.2 v_i !~~~~~~~~ cvarroot='v_i_htc' WRITE(cinteger,'(i1)') jl cvarname = TRIM(cvarroot) // TRIM(cinteger) ! WRITE(*,*), "Recording variable " // cvarname ! Request variable ID ierr = nf90_inq_varid(incid_ice_an_out, cvarname, ivarid) IF (ierr .NE. nf90_noerr ) THEN WRITE(*,*), "(sanity_check_LIM3) Error NetCDF ice file var. inquiry. Abort" WRITE(*,*), TRIM(cfileout_ice) STOP END IF ! Put variable ierr = nf90_put_var(incid_ice_an_out, ivarid, v_i(:,:,jl)) IF (ierr .NE. nf90_noerr ) THEN WRITE(*,*), "(sanity_check_LIM3) Error NetCDF ice file var. put. Abort" WRITE(*,*), TRIM(cfileout_ice) STOP END IF ! 8.1.3 v_s !~~~~~~~~ cvarroot='v_s_htc' WRITE(cinteger,'(i1)') jl cvarname = TRIM(cvarroot) // TRIM(cinteger) ! WRITE(*,*), "Recording variable " // cvarname ! Request variable ID ierr = nf90_inq_varid(incid_ice_an_out, cvarname, ivarid) IF (ierr .NE. nf90_noerr ) THEN WRITE(*,*), "(sanity_check_LIM3) Error NetCDF ice file var. inquiry. Abort" WRITE(*,*), TRIM(cfileout_ice) STOP END IF ! Put variable ierr = nf90_put_var(incid_ice_an_out, ivarid, v_s(:,:,jl)) IF (ierr .NE. nf90_noerr ) THEN WRITE(*,*), "(sanity_check_LIM3) Error NetCDF ice file var. put. Abort" WRITE(*,*), TRIM(cfileout_ice) STOP END IF ! 8.1.4 smv_i !~~~~~~~~~~ cvarroot='smv_i_htc' WRITE(cinteger,'(i1)') jl cvarname = TRIM(cvarroot) // TRIM(cinteger) ! WRITE(*,*), "Recording variable " // cvarname ! Request variable ID ierr = nf90_inq_varid(incid_ice_an_out, cvarname, ivarid) IF (ierr .NE. nf90_noerr ) THEN WRITE(*,*), "(sanity_check_LIM3) Error NetCDF ice file var. inquiry. Abort" WRITE(*,*), TRIM(cfileout_ice) STOP END IF ! Put variable ierr = nf90_put_var(incid_ice_an_out, ivarid, smv_i(:,:,jl)) IF (ierr .NE. nf90_noerr ) THEN WRITE(*,*), "(sanity_check_LIM3) Error NetCDF ice file var. put. Abort" WRITE(*,*), TRIM(cfileout_ice) STOP END IF ! 8.1.5 oa_i ! ~~~~~~~~ cvarroot='oa_i_htc' WRITE(cinteger,'(i1)') jl cvarname = TRIM(cvarroot) // TRIM(cinteger) ! WRITE(*,*), "Recording variable " // cvarname ! Request variable ID ierr = nf90_inq_varid(incid_ice_an_out, cvarname, ivarid) IF (ierr .NE. nf90_noerr ) THEN WRITE(*,*), "(sanity_check_LIM3) Error NetCDF ice file var. inquiry. Abort" WRITE(*,*), TRIM(cfileout_ice) STOP END IF ! Put variable ierr = nf90_put_var(incid_ice_an_out, ivarid, oa_i(:,:,jl)) IF (ierr .NE. nf90_noerr ) THEN WRITE(*,*), "(sanity_check_LIM3) Error NetCDF ice file var. put. Abort" WRITE(*,*), TRIM(cfileout_ice) STOP END IF ! 8.1.6 t_su ! ~~~~~~~~ cvarroot='t_su_htc' WRITE(cinteger,'(i1)') jl cvarname = TRIM(cvarroot) // TRIM(cinteger) ! WRITE(*,*), "Recording variable " // cvarname ! Request variable ID ierr = nf90_inq_varid(incid_ice_an_out, cvarname, ivarid) IF (ierr .NE. nf90_noerr ) THEN WRITE(*,*), "(sanity_check_LIM3) Error NetCDF ice file var. inquiry. Abort" WRITE(*,*), TRIM(cfileout_ice) STOP END IF ! Put variable ierr = nf90_put_var(incid_ice_an_out, ivarid, t_su(:,:,jl)) IF (ierr .NE. nf90_noerr ) THEN WRITE(*,*), "(sanity_check_LIM3) Error NetCDF ice file var. put. Abort" WRITE(*,*), TRIM(cfileout_ice) STOP END IF ! 8.1.7 Ice enthalpy (layers) ! ~~~~~~~~~~~~~~~~~~~~~~~~~ cvarroot='tempt_il' DO jk = 1 , nlay_i WRITE(cinteger2,'(i1)') jk cvarname = TRIM(cvarroot) // TRIM(cinteger2)// '_htc'//TRIM(cinteger) ! WRITE(*,*),"Recording variable " // cvarname ! Request variable ID ierr = nf90_inq_varid(incid_ice_an_out, cvarname, ivarid) IF (ierr .NE. nf90_noerr ) THEN WRITE(*,*), "(sanity_check_LIM3) Error NetCDF ice file var. inquiry. Abort" WRITE(*,*), TRIM(cfileout_ice) STOP END IF ! Put variable ierr = nf90_put_var(incid_ice_an_out, ivarid, e_i(:,:,jk,jl)) IF (ierr .NE. nf90_noerr ) THEN WRITE(*,*), "(sanity_check_LIM3) Error NetCDF ice file var. put. Abort" WRITE(*,*), TRIM(cfileout_ice) STOP END IF END DO !jk ! 8.1.8 Snow enthalpy (layers) ! ~~~~~~~~~~~~~~~~~~~~~~~~~ cvarroot='tempt_sl' DO jk = 1 , nlay_s WRITE(cinteger2,'(i1)') jk cvarname = TRIM(cvarroot) // TRIM(cinteger2)// '_htc'//TRIM(cinteger) ! WRITE(*,*),"Recording variable " // cvarname ! Request variable ID ierr = nf90_inq_varid(incid_ice_an_out, cvarname, ivarid) IF (ierr .NE. nf90_noerr ) THEN WRITE(*,*), "(sanity_check_LIM3) Error NetCDF ice file var. inquiry. Abort" WRITE(*,*), TRIM(cfileout_ice) STOP END IF ! Put variable ierr = nf90_put_var(incid_ice_an_out, ivarid, e_s(:,:,jk,jl)) IF (ierr .NE. nf90_noerr ) THEN WRITE(*,*), "(sanity_check_LIM3) Error NetCDF ice file var. put. Abort" WRITE(*,*), TRIM(cfileout_ice) STOP END IF END DO !jk END DO ! jl, categories ! Put variables that don't depend on categories ! Snow ice mass ! Request variable id cvarname = 'snwice_mass' ierr = nf90_inq_varid(incid_ice_an_out, cvarname, ivarid) IF (ierr .NE. nf90_noerr ) THEN WRITE(*,*), "(sanity_check_LIM3) Error NetCDF ice file var. inquiry. Abort" WRITE(*,*), TRIM(cfileout_ice) STOP END IF ! Put variable ierr = nf90_put_var(incid_ice_an_out, ivarid, snwice_mass(:,:)) IF (ierr .NE. nf90_noerr ) THEN WRITE(*,*), "(sanity_check_LIM3) Error NetCDF ice file var. put. Abort" WRITE(*,*), TRIM(cfileout_ice) STOP END IF cvarname = 'snwice_mass_b' ierr = nf90_inq_varid(incid_ice_an_out, cvarname, ivarid) IF (ierr .NE. nf90_noerr ) THEN WRITE(*,*), "(sanity_check_LIM3) Error NetCDF ice file var. inquiry. Abort" WRITE(*,*), TRIM(cfileout_ice) STOP END IF ! Put variable ierr = nf90_put_var(incid_ice_an_out, ivarid, snwice_mass_b(:,:)) IF (ierr .NE. nf90_noerr ) THEN WRITE(*,*), "(sanity_check_LIM3) Error NetCDF ice file var. put. Abort" WRITE(*,*), TRIM(cfileout_ice) STOP END IF ! Close the NetCDF file ierr = nf90_close(incid_ice_an_out) IF (ierr .NE. nf90_noerr ) THEN WRITE(*,*), "(sanity_check_LIM3) Error Closing NetCDF . Abort" WRITE(*,*), TRIM(cfileout_ice) STOP END IF WRITE(*,*) 'Leaving record_ice()' END SUBROUTINE record_ice SUBROUTINE record_oce() USE my_variables USE netcdf IMPLICIT NONE WRITE(*,*) 'Entering record_oce()' ! Record oce variables ! We copy the input file and store everything into this copy CALL system('cp -f '//trim(cfile_analysis_oce)//' '//trim(cfileout_oce) ) ! Open the file ierr = nf90_open(trim(cfileout_oce),nf90_Write,incid_oce_an_out) IF (ierr .NE. nf90_noerr ) THEN WRITE(*,*), "(sanity_check_LIM3) Error NetCDF Ocean file. Abort" WRITE(*,*), TRIM(cfileout_oce) STOP END IF ! A. sn (salinity) cvarname= csn ierr = nf90_inq_varid(incid_oce_an_out, cvarname, ivarid) IF (ierr .NE. nf90_noerr ) THEN WRITE(*,*), "(sanity_check_LIM3) Error NetCDF oce file var. inquiry. Abort" WRITE(*,*), TRIM(cfileout_oce) STOP END IF ! Put variable ierr = nf90_put_var(incid_oce_an_out, ivarid, sn_an(:,:,:)) IF (ierr .NE. nf90_noerr ) THEN WRITE(*,*), "(sanity_check_LIM3) Error NetCDF oce file var. put. Abort" WRITE(*,*), TRIM(cfileout_oce) STOP END IF ! Following lines commented out as the variable no longer appears ! in the restart files (2020) !! B. sss_m (sea surface salinity) !cvarname= csss_m !ierr = nf90_inq_varid(incid_oce_an_out, cvarname, ivarid) !IF (ierr .NE. nf90_noerr ) THEN ! WRITE(*,*), "(sanity_check_LIM3) Error NetCDF oce file var. inquiry. Abort" ! WRITE(*,*), TRIM(cfileout_oce) ! STOP !END IF !! Put variable !ierr = nf90_put_var(incid_oce_an_out, ivarid, sss_m_an(:,:)) !IF (ierr .NE. nf90_noerr ) THEN ! WRITE(*,*), "(sanity_check_LIM3) Error NetCDF oce file var. put. Abort" ! WRITE(*,*), TRIM(cfileout_oce) ! STOP !END IF ! C. tn (temperature) cvarname= ctn ierr = nf90_inq_varid(incid_oce_an_out, cvarname, ivarid) IF (ierr .NE. nf90_noerr ) THEN WRITE(*,*), "(sanity_check_LIM3) Error NetCDF oce file var. inquiry. Abort" WRITE(*,*), TRIM(cfileout_oce) STOP END IF ! Put variable ierr = nf90_put_var(incid_oce_an_out, ivarid, tn_an(:,:,:)) IF (ierr .NE. nf90_noerr ) THEN WRITE(*,*), "(sanity_check_LIM3) Error NetCDF oce file var. put. Abort" WRITE(*,*), TRIM(cfileout_oce) STOP END IF ! Following lines commented out as the variable no longer appears ! in the restart files (2020) !! D. sst_m (sea surface temperature) !cvarname= csst_m !ierr = nf90_inq_varid(incid_oce_an_out, cvarname, ivarid) !IF (ierr .NE. nf90_noerr ) THEN ! WRITE(*,*), "(sanity_check_LIM3) Error NetCDF oce file var. inquiry. Abort" ! WRITE(*,*), TRIM(cfileout_oce) ! STOP !END IF !! Put variable !ierr = nf90_put_var(incid_oce_an_out, ivarid, sst_m_an(:,:)) !IF (ierr .NE. nf90_noerr ) THEN ! WRITE(*,*), "(sanity_check_LIM3) Error NetCDF oce file var. put. Abort" ! WRITE(*,*), TRIM(cfileout_oce) ! STOP !END IF ! E. un (sea velocity, "un") cvarname= cun ierr = nf90_inq_varid(incid_oce_an_out, cvarname, ivarid) IF (ierr .NE. nf90_noerr ) THEN WRITE(*,*), "(sanity_check_LIM3) Error NetCDF oce file var. inquiry. Abort" WRITE(*,*), TRIM(cfileout_oce) STOP END IF ! Put variable ierr = nf90_put_var(incid_oce_an_out, ivarid, un_an(:,:,:)) IF (ierr .NE. nf90_noerr ) THEN WRITE(*,*), "(sanity_check_LIM3) Error NetCDF oce file var. put. Abort" WRITE(*,*), TRIM(cfileout_oce) STOP END IF ! F. ub (sea velocity, "ub") cvarname= cub ierr = nf90_inq_varid(incid_oce_an_out, cvarname, ivarid) IF (ierr .NE. nf90_noerr ) THEN WRITE(*,*), "(sanity_check_LIM3) Error NetCDF oce file var. inquiry. Abort" WRITE(*,*), TRIM(cfileout_oce) STOP END IF ! Put variable ierr = nf90_put_var(incid_oce_an_out, ivarid, ub_an(:,:,:)) IF (ierr .NE. nf90_noerr ) THEN WRITE(*,*), "(sanity_check_LIM3) Error NetCDF oce file var. put. Abort" WRITE(*,*), TRIM(cfileout_oce) STOP END IF ! G. vn (sea velocity, "vn") cvarname= cvn ierr = nf90_inq_varid(incid_oce_an_out, cvarname, ivarid) IF (ierr .NE. nf90_noerr ) THEN WRITE(*,*), "(sanity_check_LIM3) Error NetCDF oce file var. inquiry. Abort" WRITE(*,*), TRIM(cfileout_oce) STOP END IF ! Put variable ierr = nf90_put_var(incid_oce_an_out, ivarid, vn_an(:,:,:)) IF (ierr .NE. nf90_noerr ) THEN WRITE(*,*), "(sanity_check_LIM3) Error NetCDF oce file var. put. Abort" WRITE(*,*), TRIM(cfileout_oce) STOP END IF ! H. vb (sea velocity, "vb") cvarname= cvb ierr = nf90_inq_varid(incid_oce_an_out, cvarname, ivarid) IF (ierr .NE. nf90_noerr ) THEN WRITE(*,*), "(sanity_check_LIM3) Error NetCDF oce file var. inquiry. Abort" WRITE(*,*), TRIM(cfileout_oce) STOP END IF ! Put variable ierr = nf90_put_var(incid_oce_an_out, ivarid, vb_an(:,:,:)) IF (ierr .NE. nf90_noerr ) THEN WRITE(*,*), "(sanity_check_LIM3) Error NetCDF oce file var. put. Abort" WRITE(*,*), TRIM(cfileout_oce) STOP END IF ! Following lines commented out as the variable no longer appears ! in the restart files (2020) !! I. ssu_m (sea surface velocity, "ssu_m") !cvarname= cssu_m !ierr = nf90_inq_varid(incid_oce_an_out, cvarname, ivarid) !IF (ierr .NE. nf90_noerr ) THEN ! WRITE(*,*), "(sanity_check_LIM3) Error NetCDF oce file var. inquiry. Abort" ! WRITE(*,*), TRIM(cfileout_oce) ! STOP !END IF ! !! Put variable !ierr = nf90_put_var(incid_oce_an_out, ivarid, ssu_m_an(:,:)) !IF (ierr .NE. nf90_noerr ) THEN ! WRITE(*,*), "(sanity_check_LIM3) Error NetCDF oce file var. put. Abort" ! WRITE(*,*), TRIM(cfileout_oce) ! STOP !END IF ! !! J. ssv_m (sea surface velocity, "ssv_m") !cvarname= cssv_m !ierr = nf90_inq_varid(incid_oce_an_out, cvarname, ivarid) !IF (ierr .NE. nf90_noerr ) THEN ! WRITE(*,*), "(sanity_check_LIM3) Error NetCDF oce file var. inquiry. Abort" ! WRITE(*,*), TRIM(cfileout_oce) ! STOP !END IF ! !! Put variable !ierr = nf90_put_var(incid_oce_an_out, ivarid, ssv_m_an(:,:)) !IF (ierr .NE. nf90_noerr ) THEN ! WRITE(*,*), "(sanity_check_LIM3) Error NetCDF oce file var. put. Abort" ! WRITE(*,*), TRIM(cfileout_oce) ! STOP !END IF ! Close file ierr = nf90_close(incid_oce_an_out) IF (ierr .NE. nf90_noerr ) THEN WRITE(*,*), "(sanity_check_LIM3) Error Closing NetCDF . Abort" WRITE(*,*), TRIM(cfileout_oce) STOP END IF WRITE(*,*) 'Leaving record_oce' END SUBROUTINE record_oce SUBROUTINE salinity_check() USE my_variables USE netcdf IMPLICIT NONE ! Dummy variables INTEGER :: ji, jj, jk REAL :: zmaxsaldiff=9999.0 WRITE(*,*) 'Entering salinity_check() ' DO ji=1,jpi DO jj=1,jpj DO jk=1,jpk zsaldiff=sn_an(ji,jj,jk) - sn_fc(ji,jj,jk) IF (ABS(zsaldiff) .GT. zmaxsaldiff ) THEN IF (ABS(zsaldiff) .GT. 100.0 * zmaxsaldiff ) THEN WRITE(*,*) "Very large salinity variation" WRITE(*,*) "at point ji,jj,jk = ", ji,jj,jk WRITE(*,*) "sn_an is ",sn_an(ji,jj,jk) WRITE(*,*) "sn_fc is ",sn_fc(ji,jj,jk) WRITE(*,*) "diff is", zsaldiff END IF !WRITE(*,*) "Salinity difference at ji,jj,jk=",ji,jj,jk !WRITE(*,*) "Difference is ",ABS(zsaldiff) !WRITE(*,*) "sn_an is ",sn_an(ji,jj,jk) !WRITE(*,*) "sn_fc is ",sn_fc(ji,jj,jk) sn_an(ji,jj,jk) = sn_fc(ji,jj,jk) + SIGN(zmaxsaldiff,zsaldiff) !WRITE(*,*) "sn_an is now ",sn_an(ji,jj,jk) END IF END DO ! jk zsaldiff=(sss_m_an(ji,jj) - sss_m_fc(ji,jj)) / REAL( nn_fsbc - 1 ) IF (ABS(zsaldiff) .GT. zmaxsaldiff ) THEN IF (ABS(zsaldiff) .GT. 100.0 * zmaxsaldiff ) THEN WRITE(*,*) "Very large SSS variation at ji,jj = ",ji,jj WRITE(*,*) "sss_m_an is ", sss_m_an(ji,jj) WRITE(*,*) "sss_m_fc is ", sss_m_fc(ji,jj) WRITE(*,*) "diff is" , zsaldiff END IF !WRITE(*,*) "Sea surface salinity difference at ji,jj=",ji,jj !WRITE(*,*) "Difference is ",ABS(zsaldiff) !WRITE(*,*) "sss_m_an is ",sss_m_an(ji,jj) / REAL( nn_fsbc - 1 ) !WRITE(*,*) "sss_m_fc is ",sss_m_fc(ji,jj) / REAL( nn_fsbc - 1 ) sss_m_an(ji,jj) = sss_m_fc(ji,jj) + SIGN(REAL( nn_fsbc - 1 ) * zmaxsaldiff , zsaldiff) !WRITE(*,*) "sss_m_an is now ",sss_m_an(ji,jj) / REAL( nn_fsbc - 1 ) ! Important note ! The reason for the nn_fsbc - 1 factor is the following. In the ! restarts, the variable sss_m is (nn_fsbc -1) times the first layer of ! the sea salinity, where nn_fsbc is the frequency call of the ice to ! the ocean. As I understood, this is intended to facilitate the switch ! from one frequency call to another. The point is that the variable ! sss_m is the cumulative sea surface salinity over the (nn_fsbc-1) time ! steps. Stated the other way around, it's nn_fsbc-1 times its mean ! value. Hence here we divide sss_m by (nn_fsbc-1), ensure that ! variations in salinity are not too strong, and give back the corrected ! quantity multiplied by (nn_fsbc-1) END IF END DO END DO WRITE(*,*) 'Leaving salinity_check()' END SUBROUTINE salinity_check SUBROUTINE temperature_check() USE my_variables USE netcdf IMPLICIT NONE ! Dummy variables INTEGER :: ji, jj, jk REAL :: zmaxtempdiff=100.0 REAL :: ztempmin ! Minimum temperature for surf sea water (fct of salinity) WRITE(*,*) 'Entering temperature_check() ' ! 3D- temperature DO ji=1,jpi DO jj=1,jpj DO jk=1,jpk ztempdiff=tn_an(ji,jj,jk) - tn_fc(ji,jj,jk) IF (ABS(ztempdiff) .GT. zmaxtempdiff ) THEN IF (ABS(ztempdiff) .GT. 100.0 * zmaxtempdiff ) THEN WRITE(*,*) "Very large temperature variation detected!!" WRITE(*,*) "at point ji,jj,jk = ", ji,jj,jk WRITE(*,*) "tn_an is ",tn_an(ji,jj,jk) WRITE(*,*) "tn_fc is ",tn_fc(ji,jj,jk) END IF !WRITE(*,*) "Temperature difference at ji,jj,jk=",ji,jj,jk !WRITE(*,*) "Difference is ",ABS(ztempdiff) !WRITE(*,*) "tn_an is ",tn_an(ji,jj,jk) !WRITE(*,*) "tn_fc is ",tn_fc(ji,jj,jk) tn_an(ji,jj,jk) = tn_fc(ji,jj,jk) + SIGN(zmaxtempdiff,ztempdiff) !WRITE(*,*) "tn_an is now ",tn_an(ji,jj,jk) END IF ! if diff is large END DO ! jk DO jk =1,1 ! Reset surface temperature to freezing point if below freezing point ! This depends on the formulation chosen in the namelist (nn_eos). ! The formula comes from the NEMO code (routine eosbn2.f90) ! ! In the case nn_eos = 0 (UNESCO formulation) ! ztempmin= ( - 0.0575_wp + 1.710523e-3_wp * SQRT( sn_an(ji,jj,jk) ) & ! & - 2.154996e-4_wp * sn_an(ji,jj,jk) ) * sn_an(ji,jj,jk) ! ! In the case nn_eos = -1 or 1 (TEOS-10 formulation) r1_S0 = 0.875_wp/35.16504_wp zs = sqrt( abs(sn_an(ji,jj,jk)) * r1_S0) ztempmin = ((((1.46873e-03_wp*zs-9.64972e-03_wp)*zs+2.28348e-02_wp)*zs & & - 3.12775e-02_wp)*zs+2.07679e-02_wp)*zs-5.87701e-02_wp ztempmin = ztempmin * sn_an(ji,jj,jk) ! This is the freezing point of sea water at the surface. IF ( tn_an(ji,jj,jk) .LT. ztempmin ) THEN WRITE(*,*) 'At grid point ', ji,jj,jk WRITE(*,*) 'tn_an reset from', tn_an(ji,jj,jk),'to ', ztempmin tn_an(ji,jj,jk) = ztempmin END IF END DO ! jk ! SST ztempdiff=(sst_m_an(ji,jj) - sst_m_fc(ji,jj)) / REAL( nn_fsbc - 1 ) IF (ABS(ztempdiff) .GT. zmaxtempdiff ) THEN IF (ABS(ztempdiff) .GT. 100.0 * zmaxtempdiff ) THEN WRITE(*,*) "Very large SST variation deteced at ji,jj = ",ji,jj WRITE(*,*) "sst_m_an is ", sst_m_an(ji,jj) WRITE(*,*) "sst_m_fc is ", sst_m_fc(ji,jj) WRITE(*,*) "diff is" , ztempdiff END IF !WRITE(*,*) "Sea surface temperature difference at ji,jj=",ji,jj !WRITE(*,*) "Difference is ",ABS(ztempdiff) !WRITE(*,*) "sst_m_an is ",sst_m_an(ji,jj) / REAL( nn_fsbc - 1 ) !WRITE(*,*) "sst_m_fc is ",sst_m_fc(ji,jj) / REAL( nn_fsbc - 1 ) sst_m_an(ji,jj) = sst_m_fc(ji,jj) + SIGN(REAL( nn_fsbc - 1 ) * zmaxtempdiff , ztempdiff) !WRITE(*,*) "sst_m_an is now ",sst_m_an(ji,jj) / REAL( nn_fsbc - 1 ) ! Important note ! The reason for the nn_fsbc - 1 factor is the following. In the ! restarts, the variable sst_m is (nn_fsbc -1) times the first layer of ! the sea temperature, where nn_fsbc is the frequency call of the ice to ! the ocean. As I understood, this is intended to facilitate the switch ! from one frequency call to another. The point is that the variable ! sst_m is the cumulative sea surface temperature over the (nn_fsbc-1) time ! steps. Stated the other way around, it's nn_fsbc-1 times its mean ! value. Hence here we divide sst_m by (nn_fsbc-1), ensure that ! variations in temperature are not too strong, and give back the corrected ! quantity multiplied by (nn_fsbc-1) END IF ! For nn_eos = 0 ! ztempmin= ( - 0.0575_wp + 1.710523e-3_wp * SQRT( sn_an(ji,jj,1) ) & ! & - 2.154996e-4_wp * sn_an(ji,jj,1) ) * sn_an(ji,jj,1) ! For nn_eos = -1 or 1 r1_S0 = 0.875_wp/35.16504_wp zs = sqrt( abs(sn_an(ji,jj,1)) * r1_S0) ztempmin = ((((1.46873e-03_wp*zs-9.64972e-03_wp)*zs+2.28348e-02_wp)*zs & & - 3.12775e-02_wp)*zs+2.07679e-02_wp)*zs-5.87701e-02_wp ztempmin = ztempmin * sn_an(ji,jj,1) IF ( ( sst_m_an(ji,jj) / REAL( nn_fsbc - 1 ) ) .LT. ztempmin ) THEN WRITE(*,*), 'At grid point ', ji,jj WRITE(*,*), 'SST reset from ',sst_m_an(ji,jj)/REAL( nn_fsbc - 1 ),' to ' ,ztempmin sst_m_an(ji,jj) = ztempmin * REAL (nn_fsbc - 1) END IF END DO ! jj END DO ! ji WRITE(*,*) 'Leaving temperature_check()' END SUBROUTINE temperature_check SUBROUTINE update_hc() USE my_variables IMPLICIT NONE INTEGER :: ji, jj, jk, jl REAL(wp) :: zratio, ztmelts REAL(wp), PARAMETER :: t_init = 270.0 ! Initial temperature of ice when it is forming. Inspired from limistate. REAL(wp) :: zhc ! Heat content (in ice or snow) REAL(wp), PARAMETER :: zunit_fac= 1.0e9! This 1.0e9 is because the e_i and e_s variables are saved in Giga Joules / m2 in ! the restart (but in Joules/m2 in the code), I guess because the restart cannot take ! numbers with large exponents. ! For info: the energy to melt 1 meter of ice at 0°C is ! 334 000 [J / kg] * 1 [m] * 1000 [kg/m3] = 0.334 x 10^9 J / m2 WRITE(*,*) 'Entering update_hc()' DO jl = 1, jpl DO ji = 1, jpi DO jj = 1, jpj ! Ice layers ! Case 1: there was ice in the forecast IF (v_i_fc(ji,jj,jl) .GT. epsi10) THEN zratio = v_i(ji,jj,jl) / v_i_fc(ji,jj,jl) ! Update the ice heat content by that amount in each layer DO jk = 1, nlay_i e_i(ji,jj,jk,jl) = zratio * e_i(ji,jj,jk,jl) END DO ! jk ! Case 2: there was no ice in the forecast ELSEIF (v_i(ji,jj,jl) .GT. epsi06 ) THEN ztmelts = - tmut * smv_i(ji,jj,jl) + rtt zhc = zrhoic * (& ! [kg/m3] &zcpic * (ztmelts - t_init)& ! [J/(kg.K) ] * [K ] = J/kg &+zlfus* (1- (ztmelts - rtt)/(t_init-rtt) )& ! [J/kg]*[] &-zrcp * (ztmelts - rtt)& ! [J/(kg.K)] * [K] &) ! zhc is in J/m3 ! The amount of heat in each layer is that divided by the number of ! layers, multiplied by the sea ice volume (v_i*cell area) and by ! unit_fac (see LIM3 code) to get the units in GigaJoule DO jk = 1,nlay_i e_i(ji,jj,jk,jl) = zhc * v_i(ji,jj,jl) * zcellarea(ji,jj) / nlay_i / zunit_fac END DO !WRITE(*,*) "Ice was created by the filter at point ", ji,jj !WRITE(*,*) "Ice volume in forecast was",v_i_fc(ji,jj,jl) !WRITE(*,*) "Ice volume in analysis is ",v_i(ji,jj,jl) !WRITE(*,*) "In category ", jl !WRITE(*,*) "New e_i is ",e_i(ji,jj,:,jl) END IF !IF ( ji ==127 .AND. jj == 124 .AND. jl==5 ) THEN ! WRITE(*,*) "ji = ",ji,"jj = ",jj, "jl =", jl ! WRITE(*,*) "e_i: ", e_i(ji,jj,:,jl) ! WRITE(*,*) "zratio: ",zratio ! WRITE(*,*) "v_i", v_i(ji,jj,jl) ! WRITE(*,*) "v_i_fc: ", v_i_fc(ji,jj,jl) ! WRITE(*,*) "Stop because requested so" ! !STOP !END IF ! Snow layer ! Case 1: there was snow in the forecast IF (v_s_fc(ji,jj,jl) .GT. epsi06) THEN zratio = v_s(ji,jj,jl) / v_s_fc(ji,jj,jl) ! Update the ice heat content by that amount in each layer DO jk = 1, nlay_s e_s(ji,jj,jk,jl) = zratio * e_s(ji,jj,jk,jl) END DO ! jk ! Case 2: there was no snow in the forecast ELSEIF (v_s(ji,jj,jl) .GT. epsi06) THEN zhc = zrhosn * (& ! [kg/m3] &zcpic * (rtt - t_init)& ! [J/(kg.K) ] * [K ] = J/kg &+zlfus& ! [J/(kg)] &) ! zhc is in J/m3 ! The amount of heat in each layer is that divided by the number of ! layers, multiplied by the snow volume (v_s*cell area) DO jk = 1,nlay_s e_s(ji,jj,jk,jl) = zhc * v_s(ji,jj,jl) * zcellarea(ji,jj) / nlay_s / zunit_fac END DO END IF END DO! jj END DO ! ji END DO !jl WRITE(*,*) 'Leaving update_hc()' END SUBROUTINE update_hc SUBROUTINE regularize() USE my_variables IMPLICIT NONE INTEGER :: ji, jj, jk, jl ! Resets < 0 concentrations to 0 ! Resets variables to zero where no ice concentration DO ji=1,jpi DO jj=1,jpj DO jl =1,jpl IF ( a_i(ji,jj,jl) .LT. rzero ) THEN a_i( ji,jj,jl) = rzero v_i( ji,jj,jl) = rzero smv_i(ji,jj,jl) = rzero v_s( ji,jj,jl) = rzero oa_i( ji,jj,jl) = rzero DO jk=1,nlay_i e_i(ji,jj,jk,jl) = rzero END DO DO jk=1,nlay_s e_s(ji,jj,jk,jl) = rzero END DO END IF END DO ! Variables that do not depend on categories IF ( SUM(a_i(ji,jj,:)) .LT. rzero ) THEN snwice_mass(ji,jj) = rzero snwice_mass_b(ji,jj) = rzero END IF END DO!jj END DO ! ji END SUBROUTINE regularize SUBROUTINE velocity_check() USE my_variables IMPLICIT NONE INTEGER :: ji,jj,jk REAL :: zmaxvel = 2.0 REAL :: zmaxsurfvel= 6.0 ! Resets ocean velocities to [-2,2] ms-1 ! Resets sea surface velocities to [-4,4] ms-1 DO ji=1,jpi DO jj = 1,jpj DO jk=1,jpk IF ( ABS(un_an(ji,jj,jk)) .GT. zmaxvel ) THEN WRITE(*,*) "Fast ocean found at ji,jj,jk=",ji,jj,jk WRITE(*,*) "un_an is ",un_an(ji,jj,jk) !un_an(ji,jj,jk) = un_fc(ji,jj,jk) !zmaxvel * SIGN( rone , un_an(ji,jj,jk) ) un_an(ji,jj,jk) = zmaxvel * SIGN( rone , un_an(ji,jj,jk) ) WRITE(*,*) "un_an reset to",un_an(ji,jj,jk) END IF IF ( ABS(ub_an(ji,jj,jk)) .GT. zmaxvel ) THEN WRITE(*,*) "Fast ocean found at ji,jj,jk=",ji,jj,jk WRITE(*,*) "ub_an is ",ub_an(ji,jj,jk) !ub_an(ji,jj,jk) = ub_fc(ji,jj,jk) !zmaxvel * SIGN( rone , ub_an(ji,jj,jk) ) ub_an(ji,jj,jk) = zmaxvel * SIGN( rone , ub_an(ji,jj,jk) ) WRITE(*,*) "ub_an reset to",ub_an(ji,jj,jk) END IF IF ( ABS(vn_an(ji,jj,jk)) .GT. zmaxvel ) THEN WRITE(*,*) "Fast ocean found at ji,jj,jk=",ji,jj,jk WRITE(*,*) "vn_an is ",vn_an(ji,jj,jk) !vn_an(ji,jj,jk) = vn_fc(ji,jj,jk) !zmaxvel * SIGN( rone , vn_an(ji,jj,jk) ) vn_an(ji,jj,jk) = zmaxvel * SIGN( rone , vn_an(ji,jj,jk) ) WRITE(*,*) "vn_an reset to",vn_an(ji,jj,jk) END IF IF ( ABS(vb_an(ji,jj,jk)) .GT. zmaxvel ) THEN WRITE(*,*) "Fast ocean found at ji,jj,jk=",ji,jj,jk WRITE(*,*) "vb_an is ",vb_an(ji,jj,jk) !vb_an(ji,jj,jk) = vb_fc(ji,jj,jk) !zmaxvel * SIGN( rone , vb_an(ji,jj,jk) ) vb_an(ji,jj,jk) = zmaxvel * SIGN( rone , vb_an(ji,jj,jk) ) WRITE(*,*) "vb_an reset to",vb_an(ji,jj,jk) END IF END DO !jk ! Surface velocity IF (ABS(ssu_m_an(ji,jj)) .GT. zmaxsurfvel ) THEN WRITE(*,*) "Fast sea surface velocity found at ji,jj=" , ji,jj WRITE(*,*) "ssu_m_an is ",ssu_m_an(ji,jj) !ssu_m_an(ji,jj) = ssu_m_fc(ji,jj) !zmaxsurfvel * SIGN(rone, ssu_m_an(ji,jj) ) ssu_m_an(ji,jj) = zmaxsurfvel * SIGN(rone, ssu_m_an(ji,jj) ) WRITE(*,*) "ssu_m_an reset to", ssu_m_an(ji,jj) END IF IF (ABS(ssv_m_an(ji,jj)) .GT. zmaxsurfvel ) THEN WRITE(*,*) "Fast sea surface velocity found at ji,jj=" , ji,jj WRITE(*,*) "ssv_m_an is ",ssv_m_an(ji,jj) !ssv_m_an(ji,jj) = ssv_m_fc(ji,jj) !zmaxsurfvel * SIGN(rone, ssv_m_an(ji,jj) ) ssv_m_an(ji,jj) = zmaxsurfvel * SIGN(rone, ssv_m_an(ji,jj) ) WRITE(*,*) "ssv_m_an reset to", ssv_m_an(ji,jj) END IF END DO END DO END SUBROUTINE velocity_check