123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700701702703704705706707708709710711712713714715716717718719720721722723724725726727728729730731732733734735736737738739740741742743744745746747748749750751752753754755756757758759760761762763764765766767768769770771772773774775776777778779780781782783784785786787788789790791792793794795796797798799800801802803804805806807808809810811812813814815816817818819820821822823824825826827828829830831832833834835836837838839840841842843844845846847848849850851852853854855856857858859860861862863864865866867868869870871872873874875876877878879880881882883884885886887888889890891892893894895896897898899900901902903904905906907908909910911912913914915916917918919920921922923924925926927928929930931932933934935936937938939940941942943944945946947948949950951952953954955956957958959960961962963964965966967968969970971972973974975976977978979980981982983984985986987988989990991992993994995996997998999100010011002100310041005100610071008100910101011101210131014101510161017101810191020102110221023102410251026102710281029103010311032103310341035103610371038103910401041104210431044104510461047104810491050105110521053105410551056105710581059106010611062106310641065106610671068106910701071107210731074107510761077107810791080108110821083108410851086108710881089109010911092109310941095109610971098109911001101110211031104110511061107110811091110111111121113111411151116111711181119112011211122112311241125112611271128112911301131113211331134113511361137113811391140114111421143114411451146114711481149115011511152115311541155115611571158115911601161116211631164116511661167116811691170117111721173117411751176117711781179118011811182118311841185118611871188118911901191119211931194119511961197119811991200120112021203120412051206120712081209121012111212121312141215121612171218121912201221122212231224122512261227122812291230123112321233123412351236123712381239124012411242124312441245124612471248124912501251125212531254125512561257125812591260126112621263126412651266126712681269127012711272127312741275127612771278127912801281128212831284128512861287128812891290129112921293129412951296129712981299130013011302130313041305130613071308130913101311131213131314131513161317131813191320132113221323132413251326132713281329133013311332133313341335133613371338133913401341134213431344134513461347134813491350135113521353135413551356135713581359136013611362136313641365136613671368136913701371137213731374137513761377137813791380138113821383138413851386138713881389139013911392139313941395139613971398139914001401140214031404140514061407140814091410141114121413141414151416141714181419142014211422142314241425142614271428142914301431143214331434143514361437143814391440144114421443144414451446144714481449145014511452145314541455145614571458145914601461146214631464146514661467146814691470147114721473147414751476147714781479148014811482148314841485148614871488148914901491149214931494149514961497149814991500150115021503150415051506150715081509151015111512151315141515151615171518151915201521152215231524152515261527152815291530153115321533153415351536153715381539154015411542154315441545154615471548154915501551155215531554155515561557155815591560156115621563156415651566156715681569157015711572157315741575157615771578157915801581158215831584158515861587158815891590159115921593159415951596159715981599160016011602160316041605160616071608160916101611161216131614161516161617161816191620162116221623162416251626162716281629163016311632163316341635163616371638163916401641164216431644164516461647164816491650165116521653165416551656165716581659166016611662166316641665166616671668166916701671167216731674167516761677167816791680168116821683168416851686168716881689169016911692169316941695169616971698169917001701170217031704170517061707170817091710171117121713171417151716171717181719172017211722172317241725172617271728172917301731173217331734173517361737173817391740174117421743174417451746174717481749175017511752175317541755175617571758175917601761176217631764176517661767176817691770177117721773177417751776177717781779178017811782178317841785178617871788178917901791179217931794179517961797179817991800180118021803180418051806180718081809181018111812181318141815181618171818181918201821182218231824182518261827182818291830183118321833183418351836183718381839184018411842184318441845184618471848184918501851185218531854185518561857185818591860186118621863186418651866186718681869187018711872187318741875187618771878187918801881188218831884188518861887188818891890189118921893189418951896189718981899190019011902190319041905190619071908190919101911191219131914191519161917191819191920192119221923192419251926192719281929193019311932193319341935193619371938193919401941194219431944194519461947194819491950195119521953195419551956195719581959196019611962196319641965196619671968196919701971197219731974197519761977197819791980198119821983198419851986198719881989199019911992199319941995199619971998199920002001200220032004200520062007200820092010201120122013201420152016201720182019202020212022202320242025202620272028202920302031203220332034203520362037203820392040204120422043204420452046204720482049205020512052205320542055205620572058205920602061206220632064206520662067206820692070207120722073207420752076207720782079208020812082208320842085208620872088208920902091209220932094209520962097209820992100210121022103210421052106210721082109211021112112211321142115211621172118211921202121212221232124212521262127212821292130213121322133213421352136213721382139214021412142214321442145214621472148214921502151215221532154215521562157215821592160216121622163216421652166216721682169217021712172217321742175217621772178217921802181218221832184218521862187218821892190219121922193219421952196219721982199220022012202220322042205220622072208220922102211221222132214221522162217221822192220222122222223222422252226222722282229223022312232223322342235223622372238223922402241224222432244224522462247224822492250225122522253225422552256225722582259226022612262226322642265226622672268226922702271227222732274227522762277227822792280228122822283228422852286228722882289229022912292229322942295229622972298229923002301230223032304230523062307230823092310231123122313231423152316231723182319232023212322232323242325232623272328232923302331233223332334233523362337233823392340234123422343234423452346234723482349235023512352235323542355235623572358235923602361236223632364236523662367236823692370237123722373237423752376237723782379238023812382238323842385238623872388238923902391239223932394239523962397239823992400240124022403240424052406240724082409241024112412241324142415241624172418241924202421242224232424242524262427242824292430243124322433243424352436243724382439244024412442244324442445244624472448244924502451245224532454245524562457245824592460246124622463246424652466246724682469247024712472247324742475247624772478247924802481248224832484248524862487248824892490249124922493249424952496249724982499250025012502250325042505250625072508250925102511251225132514251525162517251825192520252125222523252425252526252725282529253025312532253325342535253625372538253925402541254225432544254525462547254825492550255125522553255425552556255725582559256025612562256325642565256625672568256925702571257225732574257525762577257825792580258125822583258425852586258725882589259025912592259325942595259625972598259926002601260226032604260526062607260826092610261126122613261426152616261726182619262026212622262326242625262626272628262926302631263226332634263526362637263826392640264126422643264426452646264726482649265026512652265326542655265626572658265926602661266226632664266526662667266826692670267126722673267426752676267726782679268026812682268326842685268626872688268926902691269226932694269526962697269826992700270127022703270427052706270727082709271027112712271327142715271627172718271927202721272227232724272527262727272827292730273127322733273427352736273727382739274027412742274327442745274627472748274927502751275227532754275527562757275827592760276127622763276427652766276727682769277027712772277327742775277627772778277927802781278227832784278527862787278827892790279127922793279427952796279727982799280028012802280328042805280628072808280928102811281228132814 |
- 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
- ! 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
- ! 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
- ! 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
- ! 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
- ! 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
- ! 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
- ! 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
- ! 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
- ! 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
|