source: LMDZ6/trunk/libf/phylmd/phyaqua_mod.F90 @ 3568

Last change on this file since 3568 was 3540, checked in by Laurent Fairhead, 5 years ago

Modifications needed for "real" calendar in 1D model
MPL/EM
  • Property copyright set to
    Name of program: LMDZ
    Creation date: 1984
    Version: LMDZ5
    License: CeCILL version 2
    Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
    See the license file in the root directory
  • Property svn:keywords set to Id
File size: 31.2 KB
Line 
1!
2! $Id: phyaqua_mod.F90 3540 2019-06-25 14:50:13Z fairhead $
3!
4MODULE phyaqua_mod
5  ! Routines complementaires pour la physique planetaire.
6  IMPLICIT NONE
7
8CONTAINS
9
10  SUBROUTINE iniaqua(nlon, iflag_phys)
11
12    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
13    ! Creation d'un etat initial et de conditions aux limites
14    ! (resp startphy.nc et limit.nc) pour des configurations idealisees
15    ! du modele LMDZ dans sa version terrestre.
16    ! iflag_phys est un parametre qui controle
17    ! iflag_phys = N
18    ! de 100 a 199 : aqua planetes avec SST forcees
19    ! N-100 determine le type de SSTs
20    ! de 200 a 299 : terra planetes avec Ts calcule
21
22    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
23
24    USE dimphy, ONLY: klon
25    USE geometry_mod, ONLY : latitude
26    USE surface_data, ONLY: type_ocean, ok_veget
27    USE pbl_surface_mod, ONLY: pbl_surface_init
28    USE fonte_neige_mod, ONLY: fonte_neige_init
29    USE phys_state_var_mod
30    USE time_phylmdz_mod, ONLY: day_ref, ndays, pdtphys, &
31                                day_ini,day_end
32    USE indice_sol_mod
33    USE nrtype, ONLY: pi
34!    USE ioipsl
35    USE mod_phys_lmdz_para, ONLY: is_master
36    USE mod_phys_lmdz_transfert_para, ONLY: bcast
37    USE mod_grid_phy_lmdz
38    USE ioipsl_getin_p_mod, ONLY : getin_p
39    USE phys_cal_mod , ONLY: year_len
40    IMPLICIT NONE
41
42    include "YOMCST.h"
43    include "clesphys.h"
44    include "dimsoil.h"
45
46    INTEGER, INTENT (IN) :: nlon, iflag_phys
47    ! IM ajout latfi, lonfi
48!    REAL, INTENT (IN) :: lonfi(nlon), latfi(nlon)
49
50    INTEGER type_profil, type_aqua
51
52    ! Ajouts initialisation des surfaces
53    REAL :: run_off_lic_0(nlon)
54    REAL :: qsolsrf(nlon, nbsrf), snsrf(nlon, nbsrf)
55    REAL :: tsoil(nlon, nsoilmx, nbsrf)
56    REAL :: tslab(nlon), seaice(nlon)
57    REAL fder(nlon)
58
59
60
61    ! Arguments :
62    ! -----------
63
64    ! integer radpas
65    INTEGER it, unit, i, k, itap
66
67    REAL rugos, albedo
68    REAL tsurf
69    REAL time, timestep, day, day0
70    REAL qsol_f
71    REAL rugsrel(nlon)
72    LOGICAL alb_ocean
73
74    CHARACTER *80 ans, file_forctl, file_fordat, file_start
75    CHARACTER *100 file, var
76    CHARACTER *2 cnbl
77
78    REAL phy_nat(nlon, year_len)
79    REAL phy_alb(nlon, year_len)
80    REAL phy_sst(nlon, year_len)
81    REAL phy_bil(nlon, year_len)
82    REAL phy_rug(nlon, year_len)
83    REAL phy_ice(nlon, year_len)
84    REAL phy_fter(nlon, year_len)
85    REAL phy_foce(nlon, year_len)
86    REAL phy_fsic(nlon, year_len)
87    REAL phy_flic(nlon, year_len)
88
89    INTEGER, SAVE :: read_climoz = 0 ! read ozone climatology
90!$OMP THREADPRIVATE(read_climoz)
91
92    ! -------------------------------------------------------------------------
93    ! declaration pour l'appel a phyredem
94    ! -------------------------------------------------------------------------
95
96    ! real pctsrf(nlon,nbsrf),ftsol(nlon,nbsrf)
97    REAL falbe(nlon, nbsrf), falblw(nlon, nbsrf)
98    ! real pbl_tke(nlon,llm,nbsrf)
99    ! real rain_fall(nlon),snow_fall(nlon)
100    ! real solsw(nlon), sollw(nlon),radsol(nlon)
101    ! real t_ancien(nlon,llm),q_ancien(nlon,llm),rnebcon(nlon,llm)
102    ! real ratqs(nlon,llm)
103    ! real clwcon(nlon,llm)
104
105    INTEGER longcles
106    PARAMETER (longcles=20)
107    REAL clesphy0(longcles)
108
109
110    ! -----------------------------------------------------------------------
111    ! dynamial tendencies :
112    ! ---------------------
113
114    INTEGER l, ierr, aslun
115
116    REAL paire
117
118    ! Local
119    CHARACTER (LEN=20) :: modname='phyaqua'
120    CHARACTER (LEN=80) :: abort_message
121
122
123    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
124    ! INITIALISATIONS
125    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
126
127    ! -----------------------------------------------------------------------
128    ! Initialisations  des constantes
129    ! -------------------------------
130
131
132    if (year_len.ne.360) then
133      write (*,*) 'iniaqua: 360 day calendar is required !'
134      stop
135    endif
136
137    type_aqua = iflag_phys/100
138    type_profil = iflag_phys - type_aqua*100
139    PRINT *, 'iniaqua:type_aqua, type_profil', type_aqua, type_profil
140
141    IF (klon/=nlon) THEN
142      WRITE (*, *) 'iniaqua: klon=', klon, ' nlon=', nlon
143      abort_message= 'probleme de dimensions dans iniaqua'
144      CALL abort_physic(modname,abort_message,1)
145    END IF
146    CALL phys_state_var_init(read_climoz)
147
148
149    read_climoz = 0
150    day0 = 217.
151    day = day0
152    it = 0
153    time = 0.
154
155    ! -----------------------------------------------------------------------
156    ! initialisations de la physique
157    ! -----------------------------------------------------------------------
158
159    day_ini = day_ref
160    day_end = day_ini + ndays
161
162    nbapp_rad = 24
163    CALL getin_p('nbapp_rad', nbapp_rad)
164
165    ! ---------------------------------------------------------------------
166    ! Creation des conditions aux limites:
167    ! ------------------------------------
168    ! Initialisations des constantes
169    ! Ajouter les manquants dans planete.def... (albedo etc)
170    co2_ppm = 348.
171    CALL getin_p('co2_ppm', co2_ppm)
172
173    solaire = 1365.
174    CALL getin_p('solaire', solaire)
175 
176    ! CALL getin('albedo',albedo) ! albedo is set below, depending on
177    ! type_aqua
178    alb_ocean = .TRUE.
179    CALL getin_p('alb_ocean', alb_ocean)
180
181    WRITE (*, *) 'iniaqua: co2_ppm=', co2_ppm
182    WRITE (*, *) 'iniaqua: solaire=', solaire
183    WRITE (*, *) 'iniaqua: alb_ocean=', alb_ocean
184
185    radsol = 0.
186    qsol_f = 10.
187
188    ! Conditions aux limites:
189    ! -----------------------
190
191    qsol(:) = qsol_f
192    rugsrel = 0.0 ! (rugsrel = rugoro)
193    rugoro = 0.0
194    u_ancien = 0.0
195    v_ancien = 0.0
196    agesno = 50.0
197    ! Relief plat
198    zmea = 0.
199    zstd = 0.
200    zsig = 0.
201    zgam = 0.
202    zthe = 0.
203    zpic = 0.
204    zval = 0.
205
206    ! Une seule surface
207    pctsrf = 0.
208    IF (type_aqua==1) THEN
209      rugos = 1.E-4
210      albedo = 0.19
211      pctsrf(:, is_oce) = 1.
212    ELSE IF (type_aqua==2) THEN
213      rugos = 0.03
214      albedo = 0.1
215      pctsrf(:, is_ter) = 1.
216    END IF
217
218    CALL getin_p('rugos', rugos)
219
220    WRITE (*, *) 'iniaqua: rugos=', rugos
221    zmasq(:) = pctsrf(:, is_ter)
222
223    ! pctsrf_pot(:,is_oce) = 1. - zmasq(:)
224    ! pctsrf_pot(:,is_sic) = 1. - zmasq(:)
225
226    ! Si alb_ocean on calcule un albedo oceanique moyen
227    ! if (alb_ocean) then
228    ! Voir pourquoi on avait ca.
229    ! CALL ini_alb_oce(phy_alb)
230    ! else
231    phy_alb(:, :) = albedo ! albedo land only (old value condsurf_jyg=0.3)
232    ! endif !alb_ocean
233
234    DO i = 1, year_len
235      ! IM Terraplanete   phy_sst(:,i) = 260.+50.*cos(rlatd(:))**2
236      ! IM ajout calcul profil sst selon le cas considere (cf. FBr)
237
238      phy_nat(:, i) = 1.0 ! 0=ocean libre, 1=land, 2=glacier, 3=banquise
239      phy_bil(:, i) = 1.0 ! ne sert que pour les slab_ocean
240      phy_rug(:, i) = rugos ! longueur rugosite utilisee sur land only
241      phy_ice(:, i) = 0.0 ! fraction de glace (?)
242      phy_fter(:, i) = pctsrf(:, is_ter) ! fraction de glace (?)
243      phy_foce(:, i) = pctsrf(:, is_oce) ! fraction de glace (?)
244      phy_fsic(:, i) = pctsrf(:, is_sic) ! fraction de glace (?)
245      phy_flic(:, i) = pctsrf(:, is_lic) ! fraction de glace (?)
246    END DO
247    ! IM calcul profil sst
248    CALL profil_sst(nlon, latitude, type_profil, phy_sst)
249
250    IF (grid_type==unstructured) THEN
251      CALL writelim_unstruct(klon, phy_nat, phy_alb, phy_sst, phy_bil, phy_rug, phy_ice, &
252                             phy_fter, phy_foce, phy_flic, phy_fsic)
253    ELSE
254     
255       CALL writelim(klon, phy_nat, phy_alb, phy_sst, phy_bil, phy_rug, phy_ice, &
256                     phy_fter, phy_foce, phy_flic, phy_fsic)
257    ENDIF
258
259    ! ---------------------------------------------------------------------
260    ! Ecriture de l'etat initial:
261    ! ---------------------------
262
263
264    ! Ecriture etat initial physique
265
266    timestep = pdtphys
267    radpas = nint(rday/timestep/float(nbapp_rad))
268
269    DO i = 1, longcles
270      clesphy0(i) = 0.
271    END DO
272    clesphy0(1) = float(iflag_con)
273    clesphy0(2) = float(nbapp_rad)
274    ! IF( cycle_diurne  ) clesphy0(3) =  1.
275    clesphy0(3) = 1. ! cycle_diurne
276    clesphy0(4) = 1. ! soil_model
277    clesphy0(5) = 1. ! new_oliq
278    clesphy0(6) = 0. ! ok_orodr
279    clesphy0(7) = 0. ! ok_orolf
280    clesphy0(8) = 0. ! ok_limitvrai
281
282
283    ! =======================================================================
284    ! Profils initiaux
285    ! =======================================================================
286
287    ! On initialise les temperatures de surfaces comme les sst
288    DO i = 1, nlon
289      ftsol(i, :) = phy_sst(i, 1)
290      tsoil(i, :, :) = phy_sst(i, 1)
291      tslab(i) = phy_sst(i, 1)
292    END DO
293
294    falbe(:, :) = albedo
295    falblw(:, :) = albedo
296    rain_fall(:) = 0.
297    snow_fall(:) = 0.
298    solsw(:) = 0.
299    sollw(:) = 0.
300    radsol(:) = 0.
301
302    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
303    ! intialisation bidon mais pas grave
304    t_ancien(:, :) = 0.
305    q_ancien(:, :) = 0.
306    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
307    rnebcon = 0.
308    ratqs = 0.
309    clwcon = 0.
310    pbl_tke = 1.E-8
311
312    ! variables supplementaires pour appel a plb_surface_init
313    fder(:) = 0.
314    seaice(:) = 0.
315    run_off_lic_0 = 0.
316    fevap = 0.
317
318
319    ! Initialisations necessaires avant phyredem
320    type_ocean = 'force'
321    CALL fonte_neige_init(run_off_lic_0)
322    qsolsrf(:, :) = qsol(1) ! humidite du sol des sous surface
323    snsrf(:, :) = 0. ! couverture de neige des sous surface
324    z0m(:, :) = rugos ! couverture de neige des sous surface
325    z0h=z0m
326
327
328    CALL pbl_surface_init(fder, snsrf, qsolsrf, tsoil)
329
330    PRINT *, 'iniaqua: before phyredem'
331
332    pbl_tke(:,:,:) = 1.e-8
333    falb1 = albedo
334    falb2 = albedo
335    zmax0 = 0.
336    f0 = 0.
337    sig1 = 0.
338    w01 = 0.
339    wake_deltat = 0.
340    wake_deltaq = 0.
341    wake_s = 0.
342    wake_dens = 0.
343    wake_cstar = 0.
344    wake_pe = 0.
345    wake_fip = 0.
346    fm_therm = 0.
347    entr_therm = 0.
348    detr_therm = 0.
349    ale_bl = 0.
350    ale_bl_trig =0.
351    alp_bl =0.
352    treedrg(:,:,:)=0.
353
354!ym error : the sub surface dimension is the third not second : forgotten for iniaqua
355!    falb_dir(:,is_ter,:)=0.08; falb_dir(:,is_lic,:)=0.6
356!    falb_dir(:,is_oce,:)=0.5;  falb_dir(:,is_sic,:)=0.6
357    falb_dir(:,:,is_ter)=0.08; falb_dir(:,:,is_lic)=0.6
358    falb_dir(:,:,is_oce)=0.5;  falb_dir(:,:,is_sic)=0.6
359
360!ym falb_dif has been forgotten, initialize with defaukt value found in phyetat0 or 0 ?
361!ym probably the uninitialized value was 0 for standard (regular grid) case
362    falb_dif(:,:,:)=0
363
364
365    CALL phyredem('startphy.nc')
366
367    PRINT *, 'iniaqua: after phyredem'
368    CALL phys_state_var_end
369
370    RETURN
371  END SUBROUTINE iniaqua
372
373
374  ! ====================================================================
375  ! ====================================================================
376  SUBROUTINE zenang_an(cycle_diurne, gmtime, rlat, rlon, rmu0, fract)
377    USE dimphy
378    IMPLICIT NONE
379    ! ====================================================================
380    ! =============================================================
381    ! CALL zenang(cycle_diurne,gmtime,rlat,rlon,rmu0,fract)
382    ! Auteur : A. Campoy et F. Hourdin
383    ! Objet  : calculer les valeurs moyennes du cos de l'angle zenithal
384    ! et l'ensoleillement moyen entre gmtime1 et gmtime2
385    ! connaissant la declinaison, la latitude et la longitude.
386
387    ! Dans cette version particuliere, on calcule le rayonnement
388    ! moyen sur l'année à chaque latitude.
389    ! angle zenithal calculé pour obtenir un
390    ! Fit polynomial de  l'ensoleillement moyen au sommet de l'atmosphere
391    ! en moyenne annuelle.
392    ! Spécifique de la terre. Utilisé pour les aqua planetes.
393
394    ! Rque   : Different de la routine angle en ce sens que zenang
395    ! fournit des moyennes de pmu0 et non des valeurs
396    ! instantanees, du coup frac prend toutes les valeurs
397    ! entre 0 et 1.
398    ! Date   : premiere version le 13 decembre 1994
399    ! revu pour  GCM  le 30 septembre 1996
400    ! ===============================================================
401    ! longi----INPUT : la longitude vraie de la terre dans son plan
402    ! solaire a partir de l'equinoxe de printemps (degre)
403    ! gmtime---INPUT : temps universel en fraction de jour
404    ! pdtrad---INPUT : pas de temps du rayonnement (secondes)
405    ! lat------INPUT : latitude en degres
406    ! long-----INPUT : longitude en degres
407    ! pmu0-----OUTPUT: angle zenithal moyen entre gmtime et gmtime+pdtrad
408    ! frac-----OUTPUT: ensoleillement moyen entre gmtime et gmtime+pdtrad
409    ! ================================================================
410    include "YOMCST.h"
411    ! ================================================================
412    LOGICAL cycle_diurne
413    REAL gmtime
414    REAL rlat(klon), rlon(klon), rmu0(klon), fract(klon)
415    ! ================================================================
416    INTEGER i
417    REAL gmtime1, gmtime2
418    REAL pi_local
419
420
421    REAL rmu0m(klon), rmu0a(klon)
422
423
424    pi_local = 4.0*atan(1.0)
425
426    ! ================================================================
427    ! Calcul de l'angle zenithal moyen sur la journee
428    ! ================================================================
429
430    DO i = 1, klon
431      fract(i) = 1.
432      ! Calcule du flux moyen
433      IF (abs(rlat(i))<=28.75) THEN
434        rmu0m(i) = (210.1924+206.6059*cos(0.0174533*rlat(i))**2)/1365.
435      ELSE IF (abs(rlat(i))<=43.75) THEN
436        rmu0m(i) = (187.4562+236.1853*cos(0.0174533*rlat(i))**2)/1365.
437      ELSE IF (abs(rlat(i))<=71.25) THEN
438        rmu0m(i) = (162.4439+284.1192*cos(0.0174533*rlat(i))**2)/1365.
439      ELSE
440        rmu0m(i) = (172.8125+183.7673*cos(0.0174533*rlat(i))**2)/1365.
441      END IF
442    END DO
443
444    ! ================================================================
445    ! Avec ou sans cycle diurne
446    ! ================================================================
447
448    IF (cycle_diurne) THEN
449
450      ! On redecompose flux  au sommet suivant un cycle diurne idealise
451      ! identique a toutes les latitudes.
452
453      DO i = 1, klon
454        rmu0a(i) = 2.*rmu0m(i)*sqrt(2.)*pi_local/(4.-pi_local)
455        rmu0(i) = rmu0a(i)*abs(sin(pi_local*gmtime+pi_local*rlon(i)/360.)) - &
456          rmu0a(i)/sqrt(2.)
457      END DO
458
459      DO i = 1, klon
460        IF (rmu0(i)<=0.) THEN
461          rmu0(i) = 0.
462          fract(i) = 0.
463        ELSE
464          fract(i) = 1.
465        END IF
466      END DO
467
468      ! Affichage de l'angel zenitale
469      ! print*,'************************************'
470      ! print*,'************************************'
471      ! print*,'************************************'
472      ! print*,'latitude=',rlat(i),'longitude=',rlon(i)
473      ! print*,'rmu0m=',rmu0m(i)
474      ! print*,'rmu0a=',rmu0a(i)
475      ! print*,'rmu0=',rmu0(i)
476
477    ELSE
478
479      DO i = 1, klon
480        fract(i) = 0.5
481        rmu0(i) = rmu0m(i)*2.
482      END DO
483
484    END IF
485
486    RETURN
487  END SUBROUTINE zenang_an
488
489  ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
490
491  SUBROUTINE writelim_unstruct(klon, phy_nat, phy_alb, phy_sst, phy_bil, phy_rug, &
492      phy_ice, phy_fter, phy_foce, phy_flic, phy_fsic)
493
494    USE mod_phys_lmdz_para, ONLY: is_omp_master, klon_mpi
495    USE mod_phys_lmdz_transfert_para, ONLY: gather_omp
496#ifdef CPP_XIOS
497    USE xios
498#endif
499    IMPLICIT NONE
500
501    include "netcdf.inc"
502
503    INTEGER, INTENT (IN) :: klon
504    REAL, INTENT (IN) :: phy_nat(klon, 360)
505    REAL, INTENT (IN) :: phy_alb(klon, 360)
506    REAL, INTENT (IN) :: phy_sst(klon, 360)
507    REAL, INTENT (IN) :: phy_bil(klon, 360)
508    REAL, INTENT (IN) :: phy_rug(klon, 360)
509    REAL, INTENT (IN) :: phy_ice(klon, 360)
510    REAL, INTENT (IN) :: phy_fter(klon, 360)
511    REAL, INTENT (IN) :: phy_foce(klon, 360)
512    REAL, INTENT (IN) :: phy_flic(klon, 360)
513    REAL, INTENT (IN) :: phy_fsic(klon, 360)
514
515    REAL :: phy_mpi(klon_mpi, 360) ! temporary variable, to store phy_***(:)
516      ! on the whole physics grid
517 
518#ifdef CPP_XIOS
519    PRINT *, 'writelim: Ecriture du fichier limit'
520
521    CALL gather_omp(phy_foce, phy_mpi)
522    IF (is_omp_master) CALL xios_send_field('foce_limout',phy_mpi)
523
524    CALL gather_omp(phy_fsic, phy_mpi)
525    IF (is_omp_master) CALL xios_send_field('fsic_limout',phy_mpi)
526     
527    CALL gather_omp(phy_fter, phy_mpi)
528    IF (is_omp_master) CALL xios_send_field('fter_limout',phy_mpi)
529     
530    CALL gather_omp(phy_flic, phy_mpi)
531    IF (is_omp_master) CALL xios_send_field('flic_limout',phy_mpi)
532
533    CALL gather_omp(phy_sst, phy_mpi)
534    IF (is_omp_master) CALL xios_send_field('sst_limout',phy_mpi)
535
536    CALL gather_omp(phy_bil, phy_mpi)
537    IF (is_omp_master) CALL xios_send_field('bils_limout',phy_mpi)
538
539    CALL gather_omp(phy_alb, phy_mpi)
540    IF (is_omp_master) CALL xios_send_field('alb_limout',phy_mpi)
541
542    CALL gather_omp(phy_rug, phy_mpi)
543    IF (is_omp_master) CALL xios_send_field('rug_limout',phy_mpi)
544#endif
545  END SUBROUTINE writelim_unstruct
546
547
548
549  SUBROUTINE writelim(klon, phy_nat, phy_alb, phy_sst, phy_bil, phy_rug, &
550      phy_ice, phy_fter, phy_foce, phy_flic, phy_fsic)
551
552    USE mod_phys_lmdz_para, ONLY: is_master
553    USE mod_grid_phy_lmdz, ONLY: klon_glo
554    USE mod_phys_lmdz_transfert_para, ONLY: gather
555    USE phys_cal_mod, ONLY: year_len
556    IMPLICIT NONE
557    include "netcdf.inc"
558
559    INTEGER, INTENT (IN) :: klon
560    REAL, INTENT (IN) :: phy_nat(klon, year_len)
561    REAL, INTENT (IN) :: phy_alb(klon, year_len)
562    REAL, INTENT (IN) :: phy_sst(klon, year_len)
563    REAL, INTENT (IN) :: phy_bil(klon, year_len)
564    REAL, INTENT (IN) :: phy_rug(klon, year_len)
565    REAL, INTENT (IN) :: phy_ice(klon, year_len)
566    REAL, INTENT (IN) :: phy_fter(klon, year_len)
567    REAL, INTENT (IN) :: phy_foce(klon, year_len)
568    REAL, INTENT (IN) :: phy_flic(klon, year_len)
569    REAL, INTENT (IN) :: phy_fsic(klon, year_len)
570
571    REAL :: phy_glo(klon_glo, year_len) ! temporary variable, to store phy_***(:)
572      ! on the whole physics grid
573    INTEGER :: k
574    INTEGER ierr
575    INTEGER dimfirst(3)
576    INTEGER dimlast(3)
577
578    INTEGER nid, ndim, ntim
579    INTEGER dims(2), debut(2), epais(2)
580    INTEGER id_tim
581    INTEGER id_nat, id_sst, id_bils, id_rug, id_alb
582    INTEGER id_fter, id_foce, id_fsic, id_flic
583
584    IF (is_master) THEN
585
586      PRINT *, 'writelim: Ecriture du fichier limit'
587
588      ierr = nf_create('limit.nc', IOR(NF_CLOBBER,NF_64BIT_OFFSET), nid)
589
590      ierr = nf_put_att_text(nid, nf_global, 'title', 30, &
591        'Fichier conditions aux limites')
592      ! !        ierr = NF_DEF_DIM (nid, "points_physiques", klon, ndim)
593      ierr = nf_def_dim(nid, 'points_physiques', klon_glo, ndim)
594      ierr = nf_def_dim(nid, 'time', nf_unlimited, ntim)
595
596      dims(1) = ndim
597      dims(2) = ntim
598
599#ifdef NC_DOUBLE
600      ierr = nf_def_var(nid, 'TEMPS', nf_double, 1, ntim, id_tim)
601#else
602      ierr = nf_def_var(nid, 'TEMPS', nf_float, 1, ntim, id_tim)
603#endif
604      ierr = nf_put_att_text(nid, id_tim, 'title', 17, 'Jour dans l annee')
605
606#ifdef NC_DOUBLE
607      ierr = nf_def_var(nid, 'NAT', nf_double, 2, dims, id_nat)
608#else
609      ierr = nf_def_var(nid, 'NAT', nf_float, 2, dims, id_nat)
610#endif
611      ierr = nf_put_att_text(nid, id_nat, 'title', 23, &
612        'Nature du sol (0,1,2,3)')
613
614#ifdef NC_DOUBLE
615      ierr = nf_def_var(nid, 'SST', nf_double, 2, dims, id_sst)
616#else
617      ierr = nf_def_var(nid, 'SST', nf_float, 2, dims, id_sst)
618#endif
619      ierr = nf_put_att_text(nid, id_sst, 'title', 35, &
620        'Temperature superficielle de la mer')
621
622#ifdef NC_DOUBLE
623      ierr = nf_def_var(nid, 'BILS', nf_double, 2, dims, id_bils)
624#else
625      ierr = nf_def_var(nid, 'BILS', nf_float, 2, dims, id_bils)
626#endif
627      ierr = nf_put_att_text(nid, id_bils, 'title', 32, &
628        'Reference flux de chaleur au sol')
629
630#ifdef NC_DOUBLE
631      ierr = nf_def_var(nid, 'ALB', nf_double, 2, dims, id_alb)
632#else
633      ierr = nf_def_var(nid, 'ALB', nf_float, 2, dims, id_alb)
634#endif
635      ierr = nf_put_att_text(nid, id_alb, 'title', 19, 'Albedo a la surface')
636
637#ifdef NC_DOUBLE
638      ierr = nf_def_var(nid, 'RUG', nf_double, 2, dims, id_rug)
639#else
640      ierr = nf_def_var(nid, 'RUG', nf_float, 2, dims, id_rug)
641#endif
642      ierr = nf_put_att_text(nid, id_rug, 'title', 8, 'Rugosite')
643
644#ifdef NC_DOUBLE
645      ierr = nf_def_var(nid, 'FTER', nf_double, 2, dims, id_fter)
646#else
647      ierr = nf_def_var(nid, 'FTER', nf_float, 2, dims, id_fter)
648#endif
649      ierr = nf_put_att_text(nid, id_fter, 'title',10,'Frac. Land')
650#ifdef NC_DOUBLE
651      ierr = nf_def_var(nid, 'FOCE', nf_double, 2, dims, id_foce)
652#else
653      ierr = nf_def_var(nid, 'FOCE', nf_float, 2, dims, id_foce)
654#endif
655      ierr = nf_put_att_text(nid, id_foce, 'title',11,'Frac. Ocean')
656#ifdef NC_DOUBLE
657      ierr = nf_def_var(nid, 'FSIC', nf_double, 2, dims, id_fsic)
658#else
659      ierr = nf_def_var(nid, 'FSIC', nf_float, 2, dims, id_fsic)
660#endif
661      ierr = nf_put_att_text(nid, id_fsic, 'title',13,'Frac. Sea Ice')
662#ifdef NC_DOUBLE
663      ierr = nf_def_var(nid, 'FLIC', nf_double, 2, dims, id_flic)
664#else
665      ierr = nf_def_var(nid, 'FLIC', nf_float, 2, dims, id_flic)
666#endif
667      ierr = nf_put_att_text(nid, id_flic, 'title',14,'Frac. Land Ice')
668
669      ierr = nf_enddef(nid)
670      IF (ierr/=nf_noerr) THEN
671        WRITE (*, *) 'writelim error: failed to end define mode'
672        WRITE (*, *) nf_strerror(ierr)
673      END IF
674
675
676      ! write the 'times'
677      DO k = 1, year_len
678#ifdef NC_DOUBLE
679        ierr = nf_put_var1_double(nid, id_tim, k, dble(k))
680#else
681        ierr = nf_put_var1_real(nid, id_tim, k, float(k))
682#endif
683        IF (ierr/=nf_noerr) THEN
684          WRITE (*, *) 'writelim error with temps(k),k=', k
685          WRITE (*, *) nf_strerror(ierr)
686        END IF
687      END DO
688
689    END IF ! of if (is_master)
690
691    ! write the fields, after having collected them on master
692
693    CALL gather(phy_nat, phy_glo)
694    IF (is_master) THEN
695#ifdef NC_DOUBLE
696      ierr = nf_put_var_double(nid, id_nat, phy_glo)
697#else
698      ierr = nf_put_var_real(nid, id_nat, phy_glo)
699#endif
700      IF (ierr/=nf_noerr) THEN
701        WRITE (*, *) 'writelim error with phy_nat'
702        WRITE (*, *) nf_strerror(ierr)
703      END IF
704    END IF
705
706    CALL gather(phy_sst, phy_glo)
707    IF (is_master) THEN
708#ifdef NC_DOUBLE
709      ierr = nf_put_var_double(nid, id_sst, phy_glo)
710#else
711      ierr = nf_put_var_real(nid, id_sst, phy_glo)
712#endif
713      IF (ierr/=nf_noerr) THEN
714        WRITE (*, *) 'writelim error with phy_sst'
715        WRITE (*, *) nf_strerror(ierr)
716      END IF
717    END IF
718
719    CALL gather(phy_bil, phy_glo)
720    IF (is_master) THEN
721#ifdef NC_DOUBLE
722      ierr = nf_put_var_double(nid, id_bils, phy_glo)
723#else
724      ierr = nf_put_var_real(nid, id_bils, phy_glo)
725#endif
726      IF (ierr/=nf_noerr) THEN
727        WRITE (*, *) 'writelim error with phy_bil'
728        WRITE (*, *) nf_strerror(ierr)
729      END IF
730    END IF
731
732    CALL gather(phy_alb, phy_glo)
733    IF (is_master) THEN
734#ifdef NC_DOUBLE
735      ierr = nf_put_var_double(nid, id_alb, phy_glo)
736#else
737      ierr = nf_put_var_real(nid, id_alb, phy_glo)
738#endif
739      IF (ierr/=nf_noerr) THEN
740        WRITE (*, *) 'writelim error with phy_alb'
741        WRITE (*, *) nf_strerror(ierr)
742      END IF
743    END IF
744
745    CALL gather(phy_rug, phy_glo)
746    IF (is_master) THEN
747#ifdef NC_DOUBLE
748      ierr = nf_put_var_double(nid, id_rug, phy_glo)
749#else
750      ierr = nf_put_var_real(nid, id_rug, phy_glo)
751#endif
752      IF (ierr/=nf_noerr) THEN
753        WRITE (*, *) 'writelim error with phy_rug'
754        WRITE (*, *) nf_strerror(ierr)
755      END IF
756    END IF
757
758    CALL gather(phy_fter, phy_glo)
759    IF (is_master) THEN
760#ifdef NC_DOUBLE
761      ierr = nf_put_var_double(nid, id_fter, phy_glo)
762#else
763      ierr = nf_put_var_real(nid, id_fter, phy_glo)
764#endif
765      IF (ierr/=nf_noerr) THEN
766        WRITE (*, *) 'writelim error with phy_fter'
767        WRITE (*, *) nf_strerror(ierr)
768      END IF
769    END IF
770
771    CALL gather(phy_foce, phy_glo)
772    IF (is_master) THEN
773#ifdef NC_DOUBLE
774      ierr = nf_put_var_double(nid, id_foce, phy_glo)
775#else
776      ierr = nf_put_var_real(nid, id_foce, phy_glo)
777#endif
778      IF (ierr/=nf_noerr) THEN
779        WRITE (*, *) 'writelim error with phy_foce'
780        WRITE (*, *) nf_strerror(ierr)
781      END IF
782    END IF
783
784    CALL gather(phy_fsic, phy_glo)
785    IF (is_master) THEN
786#ifdef NC_DOUBLE
787      ierr = nf_put_var_double(nid, id_fsic, phy_glo)
788#else
789      ierr = nf_put_var_real(nid, id_fsic, phy_glo)
790#endif
791      IF (ierr/=nf_noerr) THEN
792        WRITE (*, *) 'writelim error with phy_fsic'
793        WRITE (*, *) nf_strerror(ierr)
794      END IF
795    END IF
796
797    CALL gather(phy_flic, phy_glo)
798    IF (is_master) THEN
799#ifdef NC_DOUBLE
800      ierr = nf_put_var_double(nid, id_flic, phy_glo)
801#else
802      ierr = nf_put_var_real(nid, id_flic, phy_glo)
803#endif
804      IF (ierr/=nf_noerr) THEN
805        WRITE (*, *) 'writelim error with phy_flic'
806        WRITE (*, *) nf_strerror(ierr)
807      END IF
808    END IF
809
810    ! close file:
811    IF (is_master) THEN
812      ierr = nf_close(nid)
813    END IF
814
815  END SUBROUTINE writelim
816
817  ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
818
819  SUBROUTINE profil_sst(nlon, rlatd, type_profil, phy_sst)
820    USE dimphy
821    USE phys_cal_mod , ONLY: year_len
822    IMPLICIT NONE
823
824    INTEGER nlon, type_profil, i, k, j
825    REAL :: rlatd(nlon), phy_sst(nlon, year_len)
826    INTEGER imn, imx, amn, amx, kmn, kmx
827    INTEGER p, pplus, nlat_max
828    PARAMETER (nlat_max=72)
829    REAL x_anom_sst(nlat_max)
830    CHARACTER (LEN=20) :: modname='profil_sst'
831    CHARACTER (LEN=80) :: abort_message
832
833    IF (klon/=nlon) THEN
834       abort_message='probleme de dimensions dans profil_sst'
835       CALL abort_physic(modname,abort_message,1)
836    ENDIF
837    WRITE (*, *) ' profil_sst: type_profil=', type_profil
838    DO i = 1, year_len
839      ! phy_sst(:,i) = 260.+50.*cos(rlatd(:))**2
840
841      ! Rajout fbrlmd
842
843      IF (type_profil==1) THEN
844        ! Méthode 1 "Control" faible plateau à l'Equateur
845        DO j = 1, klon
846          phy_sst(j, i) = 273. + 27.*(1-sin(1.5*rlatd(j))**2)
847          ! PI/3=1.047197551
848          IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
849            phy_sst(j, i) = 273.
850          END IF
851        END DO
852      END IF
853      IF (type_profil==2) THEN
854        ! Méthode 2 "Flat" fort plateau à l'Equateur
855        DO j = 1, klon
856          phy_sst(j, i) = 273. + 27.*(1-sin(1.5*rlatd(j))**4)
857          IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
858            phy_sst(j, i) = 273.
859          END IF
860        END DO
861      END IF
862
863
864      IF (type_profil==3) THEN
865        ! Méthode 3 "Qobs" plateau réel à l'Equateur
866        DO j = 1, klon
867          phy_sst(j, i) = 273. + 0.5*27.*(2-sin(1.5*rlatd(j))**2-sin(1.5* &
868            rlatd(j))**4)
869          IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
870            phy_sst(j, i) = 273.
871          END IF
872        END DO
873      END IF
874
875      IF (type_profil==4) THEN
876        ! Méthode 4 : Méthode 3 + SST+2 "Qobs" plateau réel à l'Equateur
877        DO j = 1, klon
878          phy_sst(j, i) = 273. + 0.5*29.*(2-sin(1.5*rlatd(j))**2-sin(1.5* &
879            rlatd(j))**4)
880          IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
881            phy_sst(j, i) = 273.
882          END IF
883        END DO
884      END IF
885
886      IF (type_profil==5) THEN
887        ! Méthode 5 : Méthode 3 + +2K "Qobs" plateau réel à l'Equateur
888        DO j = 1, klon
889          phy_sst(j, i) = 273. + 2. + 0.5*27.*(2-sin(1.5*rlatd(j))**2-sin(1.5 &
890            *rlatd(j))**4)
891          IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
892            phy_sst(j, i) = 273. + 2.
893          END IF
894
895        END DO
896      END IF
897
898      IF (type_profil==6) THEN
899        ! Méthode 6 "cst" valeur constante de SST
900        DO j = 1, klon
901          phy_sst(j, i) = 288.
902        END DO
903      END IF
904
905
906      IF (type_profil==7) THEN
907        ! Méthode 7 "cst" valeur constante de SST +2
908        DO j = 1, klon
909          phy_sst(j, i) = 288. + 2.
910        END DO
911      END IF
912
913      p = 0
914      IF (type_profil==8) THEN
915        ! Méthode 8 profil anomalies SST du modèle couplé AR4
916        DO j = 1, klon
917          IF (rlatd(j)==rlatd(j-1)) THEN
918            phy_sst(j, i) = 273. + x_anom_sst(pplus) + &
919              0.5*27.*(2-sin(1.5*rlatd(j))**2-sin(1.5*rlatd(j))**4)
920            IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
921              phy_sst(j, i) = 273. + x_anom_sst(pplus)
922            END IF
923          ELSE
924            p = p + 1
925            pplus = 73 - p
926            phy_sst(j, i) = 273. + x_anom_sst(pplus) + &
927              0.5*27.*(2-sin(1.5*rlatd(j))**2-sin(1.5*rlatd(j))**4)
928            IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
929              phy_sst(j, i) = 273. + x_anom_sst(pplus)
930            END IF
931            WRITE (*, *) rlatd(j), x_anom_sst(pplus), phy_sst(j, i)
932          END IF
933        END DO
934      END IF
935
936      IF (type_profil==9) THEN
937        ! Méthode 5 : Méthode 3 + -2K "Qobs" plateau réel à l'Equateur
938        DO j = 1, klon
939          phy_sst(j, i) = 273. - 2. + 0.5*27.*(2-sin(1.5*rlatd(j))**2-sin(1.5 &
940            *rlatd(j))**4)
941          IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
942            phy_sst(j, i) = 273. - 2.
943          END IF
944        END DO
945      END IF
946
947
948      IF (type_profil==10) THEN
949        ! Méthode 10 : Méthode 3 + +4K "Qobs" plateau réel à l'Equateur
950        DO j = 1, klon
951          phy_sst(j, i) = 273. + 4. + 0.5*27.*(2-sin(1.5*rlatd(j))**2-sin(1.5 &
952            *rlatd(j))**4)
953          IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
954            phy_sst(j, i) = 273. + 4.
955          END IF
956        END DO
957      END IF
958
959      IF (type_profil==11) THEN
960        ! Méthode 11 : Méthode 3 + 4CO2 "Qobs" plateau réel à l'Equateur
961        DO j = 1, klon
962          phy_sst(j, i) = 273. + 0.5*27.*(2-sin(1.5*rlatd(j))**2-sin(1.5* &
963            rlatd(j))**4)
964          IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
965            phy_sst(j, i) = 273.
966          END IF
967        END DO
968      END IF
969
970      IF (type_profil==12) THEN
971        ! Méthode 12 : Méthode 10 + 4CO2 "Qobs" plateau réel à l'Equateur
972        DO j = 1, klon
973          phy_sst(j, i) = 273. + 4. + 0.5*27.*(2-sin(1.5*rlatd(j))**2-sin(1.5 &
974            *rlatd(j))**4)
975          IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
976            phy_sst(j, i) = 273. + 4.
977          END IF
978        END DO
979      END IF
980
981      IF (type_profil==13) THEN
982        ! Méthode 13 "Qmax" plateau réel à l'Equateur augmenté !
983        DO j = 1, klon
984          phy_sst(j, i) = 273. + 0.5*29.*(2-sin(1.5*rlatd(j))**2-sin(1.5* &
985            rlatd(j))**4)
986          IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
987            phy_sst(j, i) = 273.
988          END IF
989        END DO
990      END IF
991
992      IF (type_profil==14) THEN
993        ! Méthode 13 "Qmax2K" plateau réel à l'Equateur augmenté +2K !
994        DO j = 1, klon
995          phy_sst(j, i) = 273. + 2. + 0.5*29.*(2-sin(1.5*rlatd(j))**2-sin(1.5 &
996            *rlatd(j))**4)
997          IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
998            phy_sst(j, i) = 273.
999          END IF
1000        END DO
1001      END IF
1002
1003      if (type_profil.EQ.20) then
1004      print*,'Profile SST 20'
1005!     Méthode 13 "Qmax2K" plateau réel �|  l'Equateur augmenté +2K
1006
1007      do j=1,klon
1008        phy_sst(j,i)=248.+55.*(1-sin(rlatd(j))**2)
1009      enddo
1010      endif
1011
1012      if (type_profil.EQ.21) then
1013      print*,'Profile SST 21'
1014!     Méthode 13 "Qmax2K" plateau réel �|  l'Equateur augmenté +2K
1015      do j=1,klon
1016        phy_sst(j,i)=252.+55.*(1-sin(rlatd(j))**2)
1017      enddo
1018      endif
1019
1020
1021
1022    END DO
1023
1024    ! IM beg : verif profil SST: phy_sst
1025    amn = min(phy_sst(1,1), 1000.)
1026    amx = max(phy_sst(1,1), -1000.)
1027    imn = 1
1028    kmn = 1
1029    imx = 1
1030    kmx = 1
1031    DO k = 1, year_len
1032      DO i = 2, nlon
1033        IF (phy_sst(i,k)<amn) THEN
1034          amn = phy_sst(i, k)
1035          imn = i
1036          kmn = k
1037        END IF
1038        IF (phy_sst(i,k)>amx) THEN
1039          amx = phy_sst(i, k)
1040          imx = i
1041          kmx = k
1042        END IF
1043      END DO
1044    END DO
1045
1046    PRINT *, 'profil_sst: imn, kmn, phy_sst(imn,kmn) ', imn, kmn, amn
1047    PRINT *, 'profil_sst: imx, kmx, phy_sst(imx,kmx) ', imx, kmx, amx
1048    ! IM end : verif profil SST: phy_sst
1049
1050    RETURN
1051  END SUBROUTINE profil_sst
1052
1053END MODULE phyaqua_mod
Note: See TracBrowser for help on using the repository browser.